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