Intrepid2
Intrepid2_HCURL_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_HCURL_HEX_IN_FEM_DEF_HPP__
17#define __INTREPID2_HCURL_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_HCURL_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() + 3*cardLine*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
59 ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
60 ViewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
61 ViewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
62 ViewType outputBubble(Kokkos::view_wrap(ptr3, vcprop), cardBubble, npts);
63
64 // tensor product
65 ordinal_type idx = 0;
66 {
67 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
68 getValues(outputBubble, input_x, workLine, vinvBubble);
69
70 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
71 getValues(outputLine_A, input_y, workLine, vinvLine);
72
73 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
74 getValues(outputLine_B, input_z, workLine, vinvLine);
75
76 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
77 const auto output_x = outputBubble;
78 const auto output_y = outputLine_A;
79 const auto output_z = outputLine_B;
80
81 for (ordinal_type k=0;k<cardLine;++k) // z
82 for (ordinal_type j=0;j<cardLine;++j) // y
83 for (ordinal_type i=0;i<cardBubble;++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(outputLine_A, input_x, workLine, vinvLine);
93
94 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
95 getValues(outputBubble, input_y, workLine, vinvBubble);
96
97 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
98 // getValues(outputLine_B, input_z, workLine, vinvLine);
99
100 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
101 const auto output_x = outputLine_A;
102 const auto output_y = outputBubble;
103 const auto output_z = outputLine_B;
104
105 for (ordinal_type k=0;k<cardLine;++k) // z
106 for (ordinal_type j=0;j<cardBubble;++j) // y
107 for (ordinal_type i=0;i<cardLine;++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(outputLine_A, input_x, workLine, vinvLine);
117
118 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
119 getValues(outputLine_B, input_y, workLine, vinvLine);
120
121 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
122 getValues(outputBubble, input_z, workLine, vinvBubble);
123
124 // z component (bubbleBasis(z) bubbleBasis(y) lineBasis(x))
125 const auto output_x = outputLine_A;
126 const auto output_y = outputLine_B;
127 const auto output_z = outputBubble;
128
129 for (ordinal_type k=0;k<cardBubble;++k) // z
130 for (ordinal_type j=0;j<cardLine;++j) // y
131 for (ordinal_type i=0;i<cardLine;++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_CURL: {
141
142 auto ptr4 = work.data() + 4*cardLine*npts*dim_s;
143 auto ptr5 = work.data() + 5*cardLine*npts*dim_s;
144
145 ViewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
146 ViewType outputLine_A(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
147 ViewType outputLine_B(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
148 ViewType outputLine_DA(Kokkos::view_wrap(ptr3, vcprop), cardLine, npts, 1);
149 ViewType outputLine_DB(Kokkos::view_wrap(ptr4, vcprop), cardLine, npts, 1);
150 ViewType outputBubble(Kokkos::view_wrap(ptr5, vcprop), cardBubble, npts);
151
152 // tensor product
153 ordinal_type idx = 0;
154
155 { // x - component
156 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
157 getValues(outputBubble, input_x, workLine, vinvBubble);
158
159 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
160 getValues(outputLine_A, input_y, workLine, vinvLine);
161
162 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
163 getValues(outputLine_DA, input_y, workLine, vinvLine, 1);
164
165 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
166 getValues(outputLine_B, input_z, workLine, vinvLine);
167
168 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
169 getValues(outputLine_DB, input_z, workLine, vinvLine, 1);
170
171 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
172 const auto output_x = outputBubble;
173 const auto output_y = outputLine_A;
174 const auto output_dy = outputLine_DA;
175 const auto output_z = outputLine_B;
176 const auto output_dz = outputLine_DB;
177
178 for (ordinal_type k=0;k<cardLine;++k) // z
179 for (ordinal_type j=0;j<cardLine;++j) // y
180 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
181 for (ordinal_type l=0;l<npts;++l) {
182 output.access(idx,l,0) = 0.0;
183 output.access(idx,l,1) = output_x.access(i,l)*output_y.access (j,l) *output_dz.access(k,l,0);
184 output.access(idx,l,2) = -output_x.access(i,l)*output_dy.access(j,l,0)*output_z.access (k,l);
185 }
186 }
187 { // y - component
188 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
189 getValues(outputLine_A, input_x, workLine, vinvLine);
190
191 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
192 getValues(outputLine_DA, input_x, workLine, vinvLine, 1);
193
194 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
195 getValues(outputBubble, input_y, workLine, vinvBubble);
196
197 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
198 // getValues(outputLine_B, input_z, workLine, vinvLine);
199
200 //Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
201 // getValues(outputLine_DB, input_z, workLine, vinvLine, 1);
202
203 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
204 const auto output_x = outputLine_A;
205 const auto output_dx = outputLine_DA;
206 const auto output_y = outputBubble;
207 const auto output_z = outputLine_B;
208 const auto output_dz = outputLine_DB;
209
210 for (ordinal_type k=0;k<cardLine;++k) // z
211 for (ordinal_type j=0;j<cardBubble;++j) // y
212 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
213 for (ordinal_type l=0;l<npts;++l) {
214 output.access(idx,l,0) = -output_x.access (i,l) *output_y.access(j,l)*output_dz.access(k,l,0);
215 output.access(idx,l,1) = 0.0;
216 output.access(idx,l,2) = output_dx(i,l,0)*output_y.access(j,l)*output_z.access (k,l);
217 }
218 }
219 { // z - component
220 // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
221 // getValues(outputLine_A, input_x, workLine, vinvLine);
222
223 // Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
224 // getValues(outputLine_DA, input_x, workLine, vinvLine, 1);
225
226 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
227 getValues(outputLine_B, input_y, workLine, vinvLine);
228
229 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
230 getValues(outputLine_DB, input_y, workLine, vinvLine, 1);
231
232 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
233 getValues(outputBubble, input_z, workLine, vinvBubble);
234
235 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
236 const auto output_x = outputLine_A;
237 const auto output_dx = outputLine_DA;
238 const auto output_y = outputLine_B;
239 const auto output_dy = outputLine_DB;
240 const auto output_z = outputBubble;
241
242 for (ordinal_type k=0;k<cardBubble;++k) // z
243 for (ordinal_type j=0;j<cardLine;++j) // y
244 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
245 for (ordinal_type l=0;l<npts;++l) {
246 output.access(idx,l,0) = output_x.access (i,l) *output_dy.access(j,l,0)*output_z.access(k,l);
247 output.access(idx,l,1) = -output_dx(i,l,0)*output_y.access (j,l) *output_z.access(k,l);
248 output.access(idx,l,2) = 0.0;
249 }
250 }
251 break;
252 }
253 default: {
254 INTREPID2_TEST_FOR_ABORT( true,
255 ">>> ERROR: (Intrepid2::Basis_HCURL_HEX_In_FEM::Serial::getValues) operator is not supported" );
256 }
257 }
258 }
259 template<typename DT, ordinal_type numPtsPerEval,
260 typename outputValueValueType, class ...outputValueProperties,
261 typename inputPointValueType, class ...inputPointProperties,
262 typename vinvValueType, class ...vinvProperties>
263 void
264 Basis_HCURL_HEX_In_FEM::
265 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
266 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
267 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
268 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
269 const EOperator operatorType ) {
270 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
271 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
272 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
273 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
274
275 // loopSize corresponds to cardinality
276 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
277 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
278 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
279 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
280
281 typedef typename inputPointViewType::value_type inputPointType;
282
283 const ordinal_type cardinality = outputValues.extent(0);
284 //get basis order based on basis cardinality.
285 ordinal_type order = 0;
286 ordinal_type cardBubble; // = std::cbrt(cardinality/3);
287 ordinal_type cardLine; // = cardBubble+1;
288 do {
289 cardBubble = Intrepid2::getPnCardinality<1>(order);
290 cardLine = Intrepid2::getPnCardinality<1>(++order);
291 } while((3*cardBubble*cardLine*cardLine != cardinality) && (order != Parameters::MaxOrder));
292 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
293 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
294
295 switch (operatorType) {
296 case OPERATOR_VALUE: {
297 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
298 workViewType work(Kokkos::view_alloc("Basis_CURL_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
299 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
300 OPERATOR_VALUE,numPtsPerEval> FunctorType;
301 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
302 break;
303 }
304 case OPERATOR_CURL: {
305 auto workSize = Serial<OPERATOR_CURL>::getWorkSizePerPoint(order);
306 workViewType work(Kokkos::view_alloc("Basis_CURL_HEX_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
307 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
308 OPERATOR_CURL,numPtsPerEval> FunctorType;
309 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
310 break;
311 }
312 default: {
313 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
314 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Operator type not implemented" );
315 // break;commented out since exception is thrown
316 }
317 }
318 }
319 }
320
321 // -------------------------------------------------------------------------------------
322
323 template<typename DT, typename OT, typename PT>
324 Basis_HCURL_HEX_In_FEM<DT,OT,PT>::
325 Basis_HCURL_HEX_In_FEM( const ordinal_type order,
326 const EPointType pointType ) {
327
328 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
329 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
330 ">>> ERROR (Basis_HCURL_HEX_In_FEM): pointType must be either equispaced or warpblend.");
331
332 // this should be created in host and vinv should be deep copied into device space
333 Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
334 Basis_HVOL_LINE_Cn_FEM<DT,OT,PT> bubbleBasis( order - 1, POINTTYPE_GAUSS );
335
336 const ordinal_type
337 cardLine = lineBasis.getCardinality(),
338 cardBubble = bubbleBasis.getCardinality();
339
340 this->vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvLine", cardLine, cardLine);
341 this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hcurl::Hex::In::vinvBubble", cardBubble, cardBubble);
342
343 lineBasis.getVandermondeInverse(this->vinvLine_);
344 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
345
346 const ordinal_type spaceDim = 3;
347 this->basisCardinality_ = 3*cardLine*cardLine*cardBubble;
348 this->basisDegree_ = order;
349 this->basisCellTopologyKey_ = shards::Hexahedron<8>::key;
350 this->basisType_ = BASIS_FEM_LAGRANGIAN;
351 this->basisCoordinates_ = COORDINATES_CARTESIAN;
352 this->functionSpace_ = FUNCTION_SPACE_HCURL;
353 pointType_ = pointType;
354
355 // initialize tags
356 {
357 // Basis-dependent initializations
358 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
359 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
360 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
361 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
362
363 // Note: the only reason why equispaced can't support higher order than Parameters::MaxOrder appears to be the fact that the tags below get stored into a fixed-length array.
364 // TODO: relax the maximum order requirement by setting up tags in a different container, perhaps directly into an OrdinalTypeArray1DHost (tagView, below). (As of this writing (1/25/22), looks like other nodal bases do this in a similar way -- those should be fixed at the same time; maybe search for Parameters::MaxOrder.)
365 INTREPID2_TEST_FOR_EXCEPTION( order > Parameters::MaxOrder, std::invalid_argument, "polynomial order exceeds the max supported by this class");
366
367 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
368 constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
369 ordinal_type tags[3*maxCardLine*maxCardLine*maxCardLine][4];
370
371 const ordinal_type edge_x[2][2] = { {0, 4}, {2, 6} };
372 const ordinal_type edge_y[2][2] = { {3, 7}, {1, 5} };
373 const ordinal_type edge_z[2][2] = { {8,11}, {9,10} };
374
375 const ordinal_type face_yz[2] = {3, 1};
376 const ordinal_type face_xz[2] = {0, 2};
377 const ordinal_type face_xy[2] = {4, 5};
378
379 {
380 ordinal_type idx = 0;
381
385
386 // since there are x/y components in the interior
387 // dof sum should be computed before the information
388 // is assigned to tags
389 const ordinal_type
390 face_ndofs_per_direction = (cardLine-2)*cardBubble,
391 face_ndofs = 2*face_ndofs_per_direction,
392 intr_ndofs_per_direction = (cardLine-2)*(cardLine-2)*cardBubble,
393 intr_ndofs = 3*intr_ndofs_per_direction;
394
395 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
396 for (ordinal_type k=0;k<cardLine;++k) { // z
397 const auto tag_z = lineBasis.getDofTag(k);
398 for (ordinal_type j=0;j<cardLine;++j) { // y
399 const auto tag_y = lineBasis.getDofTag(j);
400 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
401 const auto tag_x = bubbleBasis.getDofTag(i);
402
403 if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 0) {
404 // edge: x edge, y vert, z vert
405 tags[idx][0] = 1; // edge dof
406 tags[idx][1] = edge_x[tag_y(1)][tag_z(1)]; // edge id
407 tags[idx][2] = tag_x(2); // local dof id
408 tags[idx][3] = tag_x(3); // total number of dofs in this edge
409 } else if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
410 // face, x edge, y vert, z edge
411 tags[idx][0] = 2; // face dof
412 tags[idx][1] = face_xz[tag_y(1)]; // face id
413 tags[idx][2] = tag_x(2) + tag_x(3)*tag_z(2); // local dof id
414 tags[idx][3] = face_ndofs; // total number of dofs in this face
415 } else if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
416 // face, x edge, y edge, z vert
417 tags[idx][0] = 2; // face dof
418 tags[idx][1] = face_xy[tag_z(1)]; // face id
419 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
420 tags[idx][3] = face_ndofs; // total number of dofs in this face
421 } else {
422 // interior
423 tags[idx][0] = 3; // interior dof
424 tags[idx][1] = 0;
425 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2) + tag_x(3)*tag_y(3)*tag_z(2); // local dof id
426 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
427 }
428 }
429 }
430 }
431
432 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
433 for (ordinal_type k=0;k<cardLine;++k) { // z
434 const auto tag_z = lineBasis.getDofTag(k);
435 for (ordinal_type j=0;j<cardBubble;++j) { // y
436 const auto tag_y = bubbleBasis.getDofTag(j);
437 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
438 const auto tag_x = lineBasis.getDofTag(i);
439
440 if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 0) {
441 // edge: x vert, y edge, z vert
442 tags[idx][0] = 1; // edge dof
443 tags[idx][1] = edge_y[tag_x(1)][tag_z(1)]; // edge id
444 tags[idx][2] = tag_y(2); // local dof id
445 tags[idx][3] = tag_y(3); // total number of dofs in this vertex
446 } else if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
447 // face, x vert, y edge, z edge
448 tags[idx][0] = 2; // face dof
449 tags[idx][1] = face_yz[tag_x(1)]; // face id
450 tags[idx][2] = tag_y(2) + tag_y(3)*tag_z(2); // local dof id
451 tags[idx][3] = face_ndofs; // total number of dofs in this vertex
452 } else if (tag_x(0) == 1 && tag_y(0) == 1 && tag_z(0) == 0) {
453 // face, x edge, y edge, z vert
454 tags[idx][0] = 2; // face dof
455 tags[idx][1] = face_xy[tag_z(1)]; // face id
456 tags[idx][2] = face_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2); // local dof id
457 tags[idx][3] = face_ndofs; // total number of dofs in this vertex
458 } else {
459 // interior
460 tags[idx][0] = 3; // interior dof
461 tags[idx][1] = 0;
462 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
463 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
464 }
465 }
466 }
467 }
468
469 // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
470 for (ordinal_type k=0;k<cardBubble;++k) { // y
471 const auto tag_z = bubbleBasis.getDofTag(k);
472 for (ordinal_type j=0;j<cardLine;++j) { // z
473 const auto tag_y = lineBasis.getDofTag(j);
474 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
475 const auto tag_x = lineBasis.getDofTag(i);
476
477 if (tag_x(0) == 0 && tag_y(0) == 0 && tag_z(0) == 1) {
478 // edge: x vert, y vert, z edge
479 tags[idx][0] = 1; // edge dof
480 tags[idx][1] = edge_z[tag_x(1)][tag_y(1)]; // edge id
481 tags[idx][2] = tag_z(2); // local dof id
482 tags[idx][3] = tag_z(3); // total number of dofs in this vertex
483 } else if (tag_x(0) == 0 && tag_y(0) == 1 && tag_z(0) == 1) {
484 // face, x vert, y edge, z edge
485 tags[idx][0] = 2; // face dof
486 tags[idx][1] = face_yz[tag_x(1)]; // face id
487 tags[idx][2] = face_ndofs_per_direction + tag_y(2) + tag_y(3)*tag_z(2); // local dof id
488 tags[idx][3] = face_ndofs; // total number of dofs in this vertex
489 } else if (tag_x(0) == 1 && tag_y(0) == 0 && tag_z(0) == 1) {
490 // face, x edge, y vert, z edge
491 tags[idx][0] = 2; // face dof
492 tags[idx][1] = face_xz[tag_y(1)]; // face id
493 tags[idx][2] = face_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_z(2); // local dof id
494 tags[idx][3] = face_ndofs; // total number of dofs in this vertex
495 } else {
496 // interior
497 tags[idx][0] = 3; // interior dof
498 tags[idx][1] = 0;
499 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
500 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
501 }
502 }
503 }
504 }
505
506 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
507 ">>> ERROR (Basis_HCURL_HEX_In_FEM): " \
508 "counted tag index is not same as cardinality." );
509 }
510
511 OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
512
513 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
514 // tags are constructed on host
515 this->setOrdinalTagData(this->tagToOrdinal_,
516 this->ordinalToTag_,
517 tagView,
518 this->basisCardinality_,
519 tagSize,
520 posScDim,
521 posScOrd,
522 posDfOrd);
523 }
524
525 // dofCoords on host and create its mirror view to device
526 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
527 dofCoordsHost("dofCoordsHost", this->basisCardinality_, spaceDim);
528
529 // dofCoords on host and create its mirror view to device
530 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
531 dofCoeffsHost("dofCoeffsHost", this->basisCardinality_, spaceDim);
532
533 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
534 dofCoordsLine("dofCoordsLine", cardLine, 1),
535 dofCoordsBubble("dofCoordsBubble", cardBubble, 1);
536
537 lineBasis.getDofCoords(dofCoordsLine);
538 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
539 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
540
541 bubbleBasis.getDofCoords(dofCoordsBubble);
542 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
543 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
544
545 {
546 ordinal_type idx = 0;
547
548 // x component (lineBasis(z) lineBasis(y) bubbleBasis(x))
549 for (ordinal_type k=0;k<cardLine;++k) { // z
550 for (ordinal_type j=0;j<cardLine;++j) { // y
551 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
552 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
553 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
554 dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
555 dofCoeffsHost(idx,0) = 1.0;
556 }
557 }
558 }
559
560 // y component (lineBasis(z) bubbleBasis(y) lineBasis(x))
561 for (ordinal_type k=0;k<cardLine;++k) { // z
562 for (ordinal_type j=0;j<cardBubble;++j) { // y
563 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
564 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
565 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
566 dofCoordsHost(idx,2) = dofCoordsLineHost(k,0);
567 dofCoeffsHost(idx,1) = 1.0;
568 }
569 }
570 }
571
572 // z component (bubbleBasis(z) lineBasis(y) lineBasis(x))
573 for (ordinal_type k=0;k<cardBubble;++k) { // z
574 for (ordinal_type j=0;j<cardLine;++j) { // y
575 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
576 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
577 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
578 dofCoordsHost(idx,2) = dofCoordsBubbleHost(k,0);
579 dofCoeffsHost(idx,2) = 1.0;
580 }
581 }
582 }
583 }
584
585 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
586 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
587
588 this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
589 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
590 }
591
592 template<typename DT, typename OT, typename PT>
593 void
594 Basis_HCURL_HEX_In_FEM<DT,OT,PT>::getScratchSpaceSize(
595 ordinal_type& perTeamSpaceSize,
596 ordinal_type& perThreadSpaceSize,
597 const PointViewType inputPoints,
598 const EOperator operatorType) const {
599 perTeamSpaceSize = 0;
600 ordinal_type scalarWorkViewExtent = (operatorType == OPERATOR_VALUE) ?
601 3*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0):
602 5*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0);
603 perThreadSpaceSize = scalarWorkViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
604 }
605
606 template<typename DT, typename OT, typename PT>
607 KOKKOS_INLINE_FUNCTION
608 void
609 Basis_HCURL_HEX_In_FEM<DT,OT,PT>::getValues(
610 OutputViewType outputValues,
611 const PointViewType inputPoints,
612 const EOperator operatorType,
613 const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
614 const typename DT::execution_space::scratch_memory_space & scratchStorage,
615 const ordinal_type subcellDim,
616 const ordinal_type subcellOrdinal) const {
617
618 INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
619 ">>> ERROR: (Intrepid2::Basis_HCURL_HEX_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
620
621 const int numPoints = inputPoints.extent(0);
622 using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
623 using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
624 ordinal_type scalarSizePerPoint = (operatorType == OPERATOR_VALUE) ?
625 3*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0):
626 5*this->vinvLine_.extent(0)+this->vinvBubble_.extent(0);
627 ordinal_type sizePerPoint = scalarSizePerPoint*get_dimension_scalar(inputPoints);
628 WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
629 using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
630
631 switch(operatorType) {
632 case OPERATOR_VALUE:
633 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinvLine_ = this->vinvLine_, &vinvBubble_ = this->vinvBubble_] (ordinal_type& pt) {
634 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
635 const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
636 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
637 Impl::Basis_HCURL_HEX_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, vinvLine_, vinvBubble_ );
638 });
639 break;
640 case OPERATOR_CURL:
641 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &vinvLine_ = this->vinvLine_, &vinvBubble_ = this->vinvBubble_] (ordinal_type& pt) {
642 auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
643 const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
644 WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
645 Impl::Basis_HCURL_HEX_In_FEM::Serial<OPERATOR_CURL>::getValues( output, input, work, vinvLine_, vinvBubble_ );
646 });
647 break;
648 default: {
649 INTREPID2_TEST_FOR_ABORT( true,
650 ">>> ERROR (Basis_HCURL_HEX_In_FEM): getValues not implemented for this operator");
651 }
652 }
653 }
654
655} // namespace Intrepid2
656
657#endif
Intrepid2::Basis_HCURL_HEX_In_FEM
Implementation of the default H(curl)-compatible FEM basis on Hexahedron cell.
Definition Intrepid2_HCURL_HEX_In_FEM.hpp:149
Intrepid2::Basis_HCURL_HEX_In_FEM::Basis_HCURL_HEX_In_FEM
Basis_HCURL_HEX_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition Intrepid2_HCURL_HEX_In_FEMDef.hpp:325
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