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