Intrepid2
Intrepid2_HDIV_HEX_In_FEMDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_HDIV_HEX_IN_FEM_DEF_HPP__
17#define __INTREPID2_HDIV_HEX_IN_FEM_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // -------------------------------------------------------------------------------------
22 namespace Impl {
23
24 template<EOperator OpType>
25 template<typename OutputViewType,
26 typename InputViewType,
27 typename WorkViewType,
28 typename VinvViewType>
29 KOKKOS_INLINE_FUNCTION
30 void
31 Basis_HDIV_HEX_In_FEM::Serial<OpType>::
32 getValues( OutputViewType output,
33 const InputViewType input,
34 WorkViewType work,
35 const VinvViewType vinvLine,
36 const VinvViewType vinvBubble) {
37 const ordinal_type cardLine = vinvLine.extent(0);
38 const ordinal_type cardBubble = vinvBubble.extent(0);
39
40 const ordinal_type npts = input.extent(0);
41
42 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
43 const auto input_x = Kokkos::subview(input, Kokkos::ALL(), range_type(0,1));
44 const auto input_y = Kokkos::subview(input, Kokkos::ALL(), range_type(1,2));
45 const auto input_z = Kokkos::subview(input, Kokkos::ALL(), range_type(2,3));
46
47 const ordinal_type dim_s = get_dimension_scalar(input);
48 auto ptr0 = work.data();
49 auto ptr1 = work.data()+cardLine*npts*dim_s;
50 auto ptr2 = work.data()+2*cardLine*npts*dim_s;
51 auto ptr3 = work.data()+(2*cardLine+cardBubble)*npts*dim_s;
52
53 typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
54 auto vcprop = Kokkos::common_view_alloc_prop(input);
55
56 switch (OpType) {
57 case OPERATOR_VALUE: {
58 ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
59 ViewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
60 ViewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
61 ViewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
62
63 // tensor product
64 ordinal_type idx = 0;
65 {
66 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
67 getValues(outputLine, input_x, workLine, vinvLine);
68
69 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
70 getValues(outputBubble_A, input_y, workLine, vinvBubble);
71
72 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
73 getValues(outputBubble_B, input_z, workLine, vinvBubble);
74
75
76 // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
77 const auto output_x = outputLine;
78 const auto output_y = outputBubble_A;
79 const auto output_z = outputBubble_B;
80
81 for (ordinal_type k=0;k<cardBubble;++k) // z
82 for (ordinal_type j=0;j<cardBubble;++j) // y
83 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
84 for (ordinal_type l=0;l<npts;++l) {
85 output.access(idx,l,0) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
86 output.access(idx,l,1) = 0.0;
87 output.access(idx,l,2) = 0.0;
88 }
89 }
90 {
91 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
92 getValues(outputBubble_A, input_x, workLine, vinvBubble);
93
94 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
95 getValues(outputLine, input_y, workLine, vinvLine);
96
97 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
98 // getValues(outputBubble_B, input_z, workLine, vinvBubble);
99
100 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
101 const auto output_x = outputBubble_A;
102 const auto output_y = outputLine;
103 const auto output_z = outputBubble_B;
104
105 for (ordinal_type k=0;k<cardBubble;++k) // z
106 for (ordinal_type j=0;j<cardLine;++j) // y
107 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
108 for (ordinal_type l=0;l<npts;++l) {
109 output.access(idx,l,0) = 0.0;
110 output.access(idx,l,1) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
111 output.access(idx,l,2) = 0.0;
112 }
113 }
114 {
115 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
116 // getValues(outputBubble_A, input_x, workLine, vinvBubble);
117
118 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
119 getValues(outputBubble_B, input_y, workLine, vinvBubble);
120
121 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
122 getValues(outputLine, input_z, workLine, vinvLine);
123
124 // z component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
125 const auto output_x = outputBubble_A;
126 const auto output_y = outputBubble_B;
127 const auto output_z = outputLine;
128
129 for (ordinal_type k=0;k<cardLine;++k) // z
130 for (ordinal_type j=0;j<cardBubble;++j) // y
131 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
132 for (ordinal_type l=0;l<npts;++l) {
133 output.access(idx,l,0) = 0.0;
134 output.access(idx,l,1) = 0.0;
135 output.access(idx,l,2) = output_x.access(i,l)*output_y.access(j,l)*output_z.access(k,l);
136 }
137 }
138 break;
139 }
140 case OPERATOR_DIV: {
141 ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
142 // A line value
143 ViewType outputBubble_A(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
144 // B line value
145 ViewType outputBubble_B(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
146 // Line grad
147 ViewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
148
149 // tensor product
150 ordinal_type idx = 0;
151
152 { // x - component
153 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
154 getValues(outputLine, input_x, workLine, vinvLine, 1);
155
156 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
157 getValues(outputBubble_A, input_y, workLine, vinvBubble);
158
159 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
160 getValues(outputBubble_B, input_z, workLine, vinvBubble);
161
162 // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
163 const auto output_dx = outputLine;
164 const auto output_y = outputBubble_A;
165 const auto output_z = outputBubble_B;
166
167 for (ordinal_type k=0;k<cardBubble;++k) // z
168 for (ordinal_type j=0;j<cardBubble;++j) // y
169 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
170 for (ordinal_type l=0;l<npts;++l)
171 output.access(idx,l) = output_dx.access(i,l,0)*output_y.access (j,l) *output_z.access(k,l);
172 }
173 { // y - component
174 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
175 getValues(outputBubble_A, input_x, workLine, vinvBubble);
176
177 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
178 getValues(outputLine, input_y, workLine, vinvLine, 1);
179
180 // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
181 // getValues(outputBubble_B, input_z, workLine, vinvBubble);
182
183 //(bubbleBasis(z) lineBasis(y) bubbleBasis(x))
184 const auto output_x = outputBubble_A;
185 const auto output_dy = outputLine;
186 const auto output_z = outputBubble_B;
187
188 for (ordinal_type k=0;k<cardBubble;++k) // z
189 for (ordinal_type j=0;j<cardLine;++j) // y
190 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
191 for (ordinal_type l=0;l<npts;++l)
192 output.access(idx,l) = output_x.access(i,l)*output_dy.access(j,l,0)*output_z.access(k,l);
193 }
194 { // z - component
195 // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
196 // getValues(outputBubble_A, input_x, workLine, vinvBubble);
197
198 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
199 getValues(outputBubble_B, input_y, workLine, vinvBubble);
200
201 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
202 getValues(outputLine, input_z, workLine, vinvLine, 1);
203
204 // (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
205 const auto output_x = outputBubble_A;
206 const auto output_y = outputBubble_B;
207 const auto output_dz = outputLine;
208
209 for (ordinal_type k=0;k<cardLine;++k) // z
210 for (ordinal_type j=0;j<cardBubble;++j) // y
211 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
212 for (ordinal_type l=0;l<npts;++l)
213 output.access(idx,l) = output_x.access(i,l)*output_y.access(j,l)*output_dz.access(k,l,0);
214 }
215 break;
216 }
217 default: {
218 INTREPID2_TEST_FOR_ABORT( true,
219 ">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::Serial::getValues) operator is not supported" );
220 }
221 }
222 }
223
224 template<typename DT, ordinal_type numPtsPerEval,
225 typename outputValueValueType, class ...outputValueProperties,
226 typename inputPointValueType, class ...inputPointProperties,
227 typename vinvValueType, class ...vinvProperties>
228 void
229 Basis_HDIV_HEX_In_FEM::
230 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
231 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
232 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
233 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
234 const EOperator operatorType ) {
235 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
236 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
237 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
238 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
239
240 // loopSize corresponds to cardinality
241 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
242 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
243 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
244 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
245
246 typedef typename inputPointViewType::value_type inputPointType;
247
248 const ordinal_type cardinality = outputValues.extent(0);
249 //get basis order based on basis cardinality.
250 ordinal_type order = 0;
251 ordinal_type cardBubble; // = std::cbrt(cardinality/3);
252 ordinal_type cardLine; // = cardBubble+1;
253 do {
254 cardBubble = Intrepid2::getPnCardinality<1>(order);
255 cardLine = Intrepid2::getPnCardinality<1>(++order);
256 } while((3*cardBubble*cardBubble*cardLine != cardinality) && (order != Parameters::MaxOrder));
257
258 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
259 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
260
261 switch (operatorType) {
262 case OPERATOR_VALUE: {
263 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
264 workViewType work(Kokkos::view_alloc("Basis_HDIV_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
265 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
266 OPERATOR_VALUE,numPtsPerEval> FunctorType;
267 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
268 break;
269 }
270 case OPERATOR_DIV: {
271 auto workSize = Serial<OPERATOR_DIV>::getWorkSizePerPoint(order);
272 workViewType work(Kokkos::view_alloc("Basis_HDIV_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
273 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
274 OPERATOR_DIV,numPtsPerEval> FunctorType;
275 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
276 break;
277 }
278 default: {
279 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
280 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Operator type not implemented" );
281 // break; commented out since exception is thrown
282 }
283 }
284 }
285 }
286
287 // -------------------------------------------------------------------------------------
288
289 template<typename DT, typename OT, typename PT>
290 Basis_HDIV_HEX_In_FEM<DT,OT,PT>::
291 Basis_HDIV_HEX_In_FEM( const ordinal_type order,
292 const EPointType pointType ) {
293
294 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
295 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
296 ">>> ERROR (Basis_HDIV_HEX_In_FEM): pointType must be either equispaced or warpblend.");
297
298 // this should be created in host and vinv should be deep copied into device space
299 Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
300 Basis_HVOL_LINE_Cn_FEM<DT,OT,PT> bubbleBasis( order - 1, POINTTYPE_GAUSS );
301
302 const ordinal_type
303 cardLine = lineBasis.getCardinality(),
304 cardBubble = bubbleBasis.getCardinality();
305
306 this->vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvLine", cardLine, cardLine);
307 this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvBubble", cardBubble, cardBubble);
308
309 lineBasis.getVandermondeInverse(this->vinvLine_);
310 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
311
312 const ordinal_type spaceDim = 3;
313 this->basisCardinality_ = 3*cardLine*cardBubble*cardBubble;
314 this->basisDegree_ = order;
315 this->basisCellTopologyKey_ = shards::Hexahedron<8>::key;
316 this->basisType_ = BASIS_FEM_LAGRANGIAN;
317 this->basisCoordinates_ = COORDINATES_CARTESIAN;
318 this->functionSpace_ = FUNCTION_SPACE_HDIV;
319 pointType_ = pointType;
320
321 // initialize tags
322 {
323 // Basis-dependent initializations
324 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
325 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
326 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
327 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
328
329 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
330 constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
331 ordinal_type tags[3*maxCardLine*maxCardLine*maxCardLine][4];
332
333 const ordinal_type face_yz[2] = {3, 1};
334 const ordinal_type face_xz[2] = {0, 2};
335 const ordinal_type face_xy[2] = {4, 5};
336
337 {
338 ordinal_type idx = 0;
339
343
344 // since there are x/y components in the interior
345 // dof sum should be computed before the information
346 // is assigned to tags
347 const ordinal_type
348 intr_ndofs_per_direction = (cardLine-2)*cardBubble*cardBubble,
349 intr_ndofs = 3*intr_ndofs_per_direction;
350
351 // x component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
352 for (ordinal_type k=0;k<cardBubble;++k) { // z
353 const auto tag_z = bubbleBasis.getDofTag(k);
354 for (ordinal_type j=0;j<cardBubble;++j) { // y
355 const auto tag_y = bubbleBasis.getDofTag(j);
356 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
357 const auto tag_x = lineBasis.getDofTag(i);
358
359 if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
360 // face, x vert, y edge, z edge
361 tags[idx][0] = 2; // face dof
362 tags[idx][1] = face_yz[tag_x(1)]; // face id
363 tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2); // local dof id
364 tags[idx][3] = tag_y(3)*tag_z(3); // total number of dofs in this vertex
365 } else {
366 // interior
367 tags[idx][0] = 3; // interior dof
368 tags[idx][1] = 0;
369 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
370 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
371 }
372 }
373 }
374 }
375
376 // y component (bubbleBasis(z) lineBasis(y) bubbleBasis(x))
377 for (ordinal_type k=0;k<cardBubble;++k) { // z
378 const auto tag_z = bubbleBasis.getDofTag(k);
379 for (ordinal_type j=0;j<cardLine;++j) { // y
380 const auto tag_y = lineBasis.getDofTag(j);
381 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
382 const auto tag_x = bubbleBasis.getDofTag(i);
383
384 if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
385 // face, x edge, y vert, z edge
386 tags[idx][0] = 2; // face dof
387 tags[idx][1] = face_xz[tag_y(1)]; // face id
388 tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2); // local dof id
389 tags[idx][3] = tag_x(3)*tag_z(3); // total number of dofs in this vertex
390 } else {
391 // interior
392 tags[idx][0] = 3; // interior dof
393 tags[idx][1] = 0;
394 tags[idx][2] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
395 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
396 }
397 }
398 }
399 }
400
401 // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
402 for (ordinal_type k=0;k<cardLine;++k) { // y
403 const auto tag_z = lineBasis.getDofTag(k);
404 for (ordinal_type j=0;j<cardBubble;++j) { // z
405 const auto tag_y = bubbleBasis.getDofTag(j);
406 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
407 const auto tag_x = bubbleBasis.getDofTag(i);
408
409 if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
410 // face, x edge, y edge, z vert
411 tags[idx][0] = 2; // face dof
412 tags[idx][1] = face_xy[tag_z(1)]; // face id
413 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
414 tags[idx][3] = tag_x(3)*tag_y(3); // total number of dofs in this vertex
415 } else {
416 // interior
417 tags[idx][0] = 3; // interior dof
418 tags[idx][1] = 0;
419 tags[idx][2] = 2*intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
420 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
421 }
422 }
423 }
424 }
425
426 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
427 ">>> ERROR (Basis_HDIV_HEX_In_FEM): " \
428 "counted tag index is not same as cardinality." );
429 }
430
431 OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
432
433 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
434 // tags are constructed on host
435 this->setOrdinalTagData(this->tagToOrdinal_,
436 this->ordinalToTag_,
437 tagView,
438 this->basisCardinality_,
439 tagSize,
440 posScDim,
441 posScOrd,
442 posDfOrd);
443 }
444
445 // dofCoords on host and create its mirror view to device
446 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
447 dofCoordsHost("dofCoordsHost", this->basisCardinality_, spaceDim);
448
449 // dofCoeffs on host and create its mirror view to device
450 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
451 dofCoeffsHost("dofCoeffsHost", this->basisCardinality_, spaceDim);
452
453 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
454 dofCoordsLine("dofCoordsLine", cardLine, 1),
455 dofCoordsBubble("dofCoordsBubble", cardBubble, 1);
456
457 lineBasis.getDofCoords(dofCoordsLine);
458 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
459 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
460
461 bubbleBasis.getDofCoords(dofCoordsBubble);
462 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
463 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
464
465 {
466 ordinal_type idx = 0;
467
468 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
469 for (ordinal_type k=0;k<cardBubble;++k) { // z
470 for (ordinal_type j=0;j<cardBubble;++j) { // y
471 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
472 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
473 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
474 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
475 dofCoeffsHost(idx,0) = 1.0;
476 }
477 }
478 }
479
480 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
481 for (ordinal_type k=0;k<cardBubble;++k) { // z
482 for (ordinal_type j=0;j<cardLine;++j) { // y
483 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
484 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
485 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
486 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
487 dofCoeffsHost(idx,1) = 1.0;
488 }
489 }
490 }
491
492 // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
493 for (ordinal_type k=0;k<cardLine;++k) { // z
494 for (ordinal_type j=0;j<cardBubble;++j) { // y
495 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
496 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
497 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
498 dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
499 dofCoeffsHost(idx,2) = 1.0;
500 }
501 }
502 }
503 }
504
505 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
506 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
507
508 this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
509 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
510 }
511
512 template<typename DT, typename OT, typename PT>
513 void
514 Basis_HDIV_HEX_In_FEM<DT,OT,PT>::getScratchSpaceSize(
515 ordinal_type& perTeamSpaceSize,
516 ordinal_type& perThreadSpaceSize,
517 const PointViewType inputPoints,
518 const EOperator operatorType) const {
519 perTeamSpaceSize = 0;
520 perThreadSpaceSize = (2*this->vinvLine_.extent(0)+2*this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
521 }
522
523 template<typename DT, typename OT, typename PT>
524 KOKKOS_INLINE_FUNCTION
525 void
526 Basis_HDIV_HEX_In_FEM<DT,OT,PT>::getValues(
527 OutputViewType outputValues,
528 const PointViewType inputPoints,
529 const EOperator operatorType,
530 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
531 const typename DT::execution_space::scratch_memory_space & scratchStorage,
532 const ordinal_type subcellDim,
533 const ordinal_type subcellOrdinal) const {
534
535 INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
536 ">>> ERROR: (Intrepid2::Basis_HDIV_HEX_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
537
538 const int numPoints = inputPoints.extent(0);
539 using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
540 using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
541 ordinal_type sizePerPoint = (2*this->vinvLine_.extent(0)+2*this->vinvBubble_.extent(0))*get_dimension_scalar(inputPoints);
542 WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
543 using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
544
545 switch(operatorType) {
546 case OPERATOR_VALUE:
547 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinvLine_ = this->vinvLine_, &vinvBubble_ = this->vinvBubble_] (ordinal_type& pt) {
548 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
549 const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
550 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
551 Impl::Basis_HDIV_HEX_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, vinvLine_, vinvBubble_ );
552 });
553 break;
554 case OPERATOR_DIV:
555 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinvLine_ = this->vinvLine_, &vinvBubble_ = this->vinvBubble_] (ordinal_type& pt) {
556 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
557 const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
558 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
559 Impl::Basis_HDIV_HEX_In_FEM::Serial<OPERATOR_DIV>::getValues( output, input, work, vinvLine_, vinvBubble_ );
560 });
561 break;
562 default: {
563 INTREPID2_TEST_FOR_ABORT( true,
564 ">>> ERROR (Basis_HDIV_HEX_In_FEM): getValues not implemented for this operator");
565 }
566 }
567 }
568
569} // namespace Intrepid2
570
571#endif
Intrepid2::Basis_HDIV_HEX_In_FEM
Implementation of the default H(div)-compatible FEM basis on Hexahedron cell.
Definition Intrepid2_HDIV_HEX_In_FEM.hpp:139
Intrepid2::Basis_HDIV_HEX_In_FEM::Basis_HDIV_HEX_In_FEM
Basis_HDIV_HEX_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition Intrepid2_HDIV_HEX_In_FEMDef.hpp:291
Intrepid2::Basis_HGRAD_LINE_Cn_FEM
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
Definition Intrepid2_HGRAD_LINE_Cn_FEM.hpp:144
Intrepid2::Basis_HGRAD_LINE_Cn_FEM::getDofCoords
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
Definition Intrepid2_HGRAD_LINE_Cn_FEM.hpp:214
Intrepid2::Basis_HVOL_LINE_Cn_FEM
Implementation of the locally HVOL-compatible FEM basis of variable order on the [-1,...
Definition Intrepid2_HVOL_LINE_Cn_FEM.hpp:145
Intrepid2::Basis_HVOL_LINE_Cn_FEM::getDofCoords
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
Definition Intrepid2_HVOL_LINE_Cn_FEM.hpp:205
Intrepid2::Basis::getDofTag
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
Definition Intrepid2_Basis.hpp:941
Intrepid2::Basis::getCardinality
ordinal_type getCardinality() const
Returns cardinality of the basis.
Definition Intrepid2_Basis.hpp:805
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition Intrepid2_Types.hpp:102
Generated by doxygen 1.9.8