71 using ExecutionSpace =
typename DeviceType::execution_space;
72 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
73 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
74 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
76 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
77 using TeamMember =
typename TeamPolicy::member_type;
81 OutputFieldType output_;
82 InputPointsType inputPoints_;
84 ordinal_type polyOrder_;
85 ordinal_type numFields_, numPoints_;
87 size_t fad_size_output_;
89 static constexpr ordinal_type numVertices = 4;
90 static constexpr ordinal_type numFaces = 4;
91 static constexpr ordinal_type numVerticesPerFace = 3;
92 static constexpr ordinal_type numInteriorFamilies = 3;
95 const ordinal_type face_vertices[numFaces*numVerticesPerFace] = {0,1,2,
101 const ordinal_type numFaceFunctionsPerFace_;
102 const ordinal_type numFaceFunctions_;
103 const ordinal_type numInteriorFunctionsPerFamily_;
104 const ordinal_type numInteriorFunctions_;
107 const ordinal_type faceOrdinalForInterior_[numInteriorFamilies] = {0,2,3};
108 const ordinal_type faceFamilyForInterior_[numInteriorFamilies] = {0,0,1};
109 const ordinal_type interiorCoordinateOrdinal_[numInteriorFamilies] = {3,0,1};
112 const ordinal_type interior_face_family_start_ [numInteriorFamilies] = {0,0,1};
113 const ordinal_type interior_face_family_middle_[numInteriorFamilies] = {1,1,2};
114 const ordinal_type interior_face_family_end_ [numInteriorFamilies] = {2,2,0};
116 KOKKOS_INLINE_FUNCTION
117 ordinal_type dofOrdinalForFace(
const ordinal_type &faceOrdinal,
118 const ordinal_type &zeroBasedFaceFamily,
119 const ordinal_type &i,
120 const ordinal_type &j)
const
123 const ordinal_type faceDofOffset = faceOrdinal * numFaceFunctionsPerFace_;
129 const ordinal_type max_ij_sum = polyOrder_ - 1;
131 ordinal_type fieldOrdinal = faceDofOffset + zeroBasedFaceFamily;
133 for (ordinal_type ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
135 for (ordinal_type ii=0; ii<ij_sum; ii++)
138 const ordinal_type jj = ij_sum - ii;
139 if ( (ii == i) && (jj == j))
151 : opType_(opType), output_(output), inputPoints_(inputPoints),
152 polyOrder_(polyOrder),
154 numFaceFunctionsPerFace_(polyOrder * (polyOrder+1)/2),
155 numFaceFunctions_(numFaceFunctionsPerFace_*numFaces),
156 numInteriorFunctionsPerFamily_((polyOrder-1)*polyOrder*(polyOrder+1)/6),
157 numInteriorFunctions_(numInteriorFunctionsPerFamily_ * numInteriorFamilies)
159 numFields_ = output.extent_int(0);
160 numPoints_ = output.extent_int(1);
162 const ordinal_type expectedCardinality = numFaceFunctions_ + numInteriorFunctions_;
168 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
169 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument,
"output field size does not match basis cardinality");
172 KOKKOS_INLINE_FUNCTION
173 void computeFaceJacobi(OutputScratchView &P_2ip1,
174 const ordinal_type &zeroBasedFaceOrdinal,
175 const ordinal_type &i,
176 const PointScalar* lambda)
const
179 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
180 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
181 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];
183 const auto & s0 = lambda[s0_index];
184 const auto & s1 = lambda[s1_index];
185 const auto & s2 = lambda[s2_index];
186 const PointScalar jacobiScaling = s0 + s1 + s2;
188 const double alpha = i*2.0 + 1;
189 Polynomials::shiftedScaledJacobiValues(P_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);
193 KOKKOS_INLINE_FUNCTION
195 const ordinal_type &zeroBasedInteriorFamilyOrdinal,
196 const ordinal_type &i,
197 const PointScalar* lambda)
const
199 const ordinal_type & relatedFaceOrdinal = faceOrdinalForInterior_[zeroBasedInteriorFamilyOrdinal];
200 const auto &s0_vertex_number = interior_face_family_start_ [zeroBasedInteriorFamilyOrdinal];
201 const auto &s1_vertex_number = interior_face_family_middle_[zeroBasedInteriorFamilyOrdinal];
202 const auto &s2_vertex_number = interior_face_family_end_ [zeroBasedInteriorFamilyOrdinal];
205 const auto &s0_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s0_vertex_number];
206 const auto &s1_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s1_vertex_number];
207 const auto &s2_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s2_vertex_number];
209 const auto & s0 = lambda[s0_index];
210 const auto & s1 = lambda[s1_index];
211 const auto & s2 = lambda[s2_index];
212 const PointScalar jacobiScaling = s0 + s1 + s2;
214 const double alpha = i*2.0 + 1;
215 Polynomials::shiftedScaledJacobiValues(P_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);
218 KOKKOS_INLINE_FUNCTION
219 void computeFaceLegendre(OutputScratchView &P,
220 const ordinal_type &zeroBasedFaceOrdinal,
221 const PointScalar* lambda)
const
224 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
225 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
227 const auto & s0 = lambda[s0_index];
228 const auto & s1 = lambda[s1_index];
229 const PointScalar legendreScaling = s0 + s1;
231 Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, s1, legendreScaling);
234 KOKKOS_INLINE_FUNCTION
235 void computeFaceLegendreForInterior(OutputScratchView &P,
236 const ordinal_type &zeroBasedInteriorFamilyOrdinal,
237 const PointScalar* lambda)
const
239 const ordinal_type & relatedFaceOrdinal = faceOrdinalForInterior_[zeroBasedInteriorFamilyOrdinal];
240 const auto &s0_vertex_number = interior_face_family_start_ [zeroBasedInteriorFamilyOrdinal];
241 const auto &s1_vertex_number = interior_face_family_middle_[zeroBasedInteriorFamilyOrdinal];
244 const auto &s0_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s0_vertex_number];
245 const auto &s1_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s1_vertex_number];
247 const auto & s0 = lambda[s0_index];
248 const auto & s1 = lambda[s1_index];
249 const PointScalar legendreScaling = s0 + s1;
251 Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, s1, legendreScaling);
254 KOKKOS_INLINE_FUNCTION
255 void computeFaceVectorWeight(OutputScalar &vectorWeight_x,
256 OutputScalar &vectorWeight_y,
257 OutputScalar &vectorWeight_z,
258 const ordinal_type &zeroBasedFaceOrdinal,
259 const PointScalar* lambda,
260 const PointScalar* lambda_dx,
261 const PointScalar* lambda_dy,
262 const PointScalar* lambda_dz)
const
266 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
267 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
268 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];
270 const auto & s0 = lambda [s0_index];
271 const auto & s0_dx = lambda_dx[s0_index];
272 const auto & s0_dy = lambda_dy[s0_index];
273 const auto & s0_dz = lambda_dz[s0_index];
275 const auto & s1 = lambda [s1_index];
276 const auto & s1_dx = lambda_dx[s1_index];
277 const auto & s1_dy = lambda_dy[s1_index];
278 const auto & s1_dz = lambda_dz[s1_index];
280 const auto & s2 = lambda [s2_index];
281 const auto & s2_dx = lambda_dx[s2_index];
282 const auto & s2_dy = lambda_dy[s2_index];
283 const auto & s2_dz = lambda_dz[s2_index];
285 vectorWeight_x = s0 * (s1_dy * s2_dz - s1_dz * s2_dy)
286 + s1 * (s2_dy * s0_dz - s2_dz * s0_dy)
287 + s2 * (s0_dy * s1_dz - s0_dz * s1_dy);
289 vectorWeight_y = s0 * (s1_dz * s2_dx - s1_dx * s2_dz)
290 + s1 * (s2_dz * s0_dx - s2_dx * s0_dz)
291 + s2 * (s0_dz * s1_dx - s0_dx * s1_dz);
293 vectorWeight_z = s0 * (s1_dx * s2_dy - s1_dy * s2_dx)
294 + s1 * (s2_dx * s0_dy - s2_dy * s0_dx)
295 + s2 * (s0_dx * s1_dy - s0_dy * s1_dx);
299 KOKKOS_INLINE_FUNCTION
300 void faceFunctionValue(OutputScalar &value_x,
301 OutputScalar &value_y,
302 OutputScalar &value_z,
303 const ordinal_type &i,
304 const ordinal_type &j,
305 const OutputScratchView &P,
306 const OutputScratchView &P_2ip1,
307 const OutputScalar &vectorWeight_x,
308 const OutputScalar &vectorWeight_y,
309 const OutputScalar &vectorWeight_z,
310 const PointScalar* lambda)
const
312 const auto &P_i = P(i);
313 const auto &P_2ip1_j = P_2ip1(j);
315 value_x = P_i * P_2ip1_j * vectorWeight_x;
316 value_y = P_i * P_2ip1_j * vectorWeight_y;
317 value_z = P_i * P_2ip1_j * vectorWeight_z;
320 KOKKOS_INLINE_FUNCTION
321 void computeFaceDivWeight(OutputScalar &divWeight,
322 const ordinal_type &zeroBasedFaceOrdinal,
323 const PointScalar* lambda_dx,
324 const PointScalar* lambda_dy,
325 const PointScalar* lambda_dz)
const
329 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
330 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
331 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];
333 const auto & s0_dx = lambda_dx[s0_index];
334 const auto & s0_dy = lambda_dy[s0_index];
335 const auto & s0_dz = lambda_dz[s0_index];
337 const auto & s1_dx = lambda_dx[s1_index];
338 const auto & s1_dy = lambda_dy[s1_index];
339 const auto & s1_dz = lambda_dz[s1_index];
341 const auto & s2_dx = lambda_dx[s2_index];
342 const auto & s2_dy = lambda_dy[s2_index];
343 const auto & s2_dz = lambda_dz[s2_index];
345 divWeight = s0_dx * (s1_dy * s2_dz - s1_dz * s2_dy)
346 + s0_dy * (s1_dz * s2_dx - s1_dx * s2_dz)
347 + s0_dz * (s1_dx * s2_dy - s1_dy * s2_dx);
350 KOKKOS_INLINE_FUNCTION
351 void computeInteriorIntegratedJacobi(OutputScratchView &L_2ipjp1,
352 const ordinal_type &i,
353 const ordinal_type &j,
354 const ordinal_type &zeroBasedFamilyOrdinal,
355 const PointScalar* lambda)
const
357 const auto &lambda_m = lambda[interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal]];
359 const double alpha = 2 * (i + j + 1);
361 const PointScalar jacobiScaling = 1.0;
362 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipjp1, alpha, polyOrder_-1, lambda_m, jacobiScaling);
365 KOKKOS_INLINE_FUNCTION
366 void computeInteriorJacobi(OutputScratchView &P_2ipjp1,
367 const ordinal_type &i,
368 const ordinal_type &j,
369 const ordinal_type &zeroBasedFamilyOrdinal,
370 const PointScalar* lambda)
const
372 const auto &lambda_m = lambda[interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal]];
374 const double alpha = 2 * (i + j + 1);
376 const PointScalar jacobiScaling = 1.0;
377 Polynomials::shiftedScaledJacobiValues(P_2ipjp1, alpha, polyOrder_-1, lambda_m, jacobiScaling);
381 KOKKOS_INLINE_FUNCTION
382 void faceFunctionDiv(OutputScalar &divValue,
383 const ordinal_type &i,
384 const ordinal_type &j,
385 const OutputScratchView &P,
386 const OutputScratchView &P_2ip1,
387 const OutputScalar &divWeight,
388 const PointScalar* lambda)
const
390 const auto &P_i = P(i);
391 const auto &P_2ip1_j = P_2ip1(j);
393 divValue = (i + j + 3.) * P_i * P_2ip1_j * divWeight;
397 KOKKOS_INLINE_FUNCTION
398 void gradInteriorIntegratedJacobi(OutputScalar &L_2ipjp1_dx,
399 OutputScalar &L_2ipjp1_dy,
400 OutputScalar &L_2ipjp1_dz,
401 const ordinal_type &zeroBasedFamilyOrdinal,
402 const ordinal_type &j,
403 const ordinal_type &k,
404 const OutputScratchView &P_2ipjp1,
405 const PointScalar* lambda,
406 const PointScalar* lambda_dx,
407 const PointScalar* lambda_dy,
408 const PointScalar* lambda_dz)
const
413 const ordinal_type &m = interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal];
415 L_2ipjp1_dx = P_2ipjp1(k-1) * lambda_dx[m];
416 L_2ipjp1_dy = P_2ipjp1(k-1) * lambda_dy[m];
417 L_2ipjp1_dz = P_2ipjp1(k-1) * lambda_dz[m];
420 KOKKOS_INLINE_FUNCTION
421 void interiorFunctionDiv(OutputScalar &outputDiv,
422 OutputScalar &L_2ipjp1_k,
423 OutputScalar &faceDiv,
424 OutputScalar &L_2ipjp1_k_dx,
425 OutputScalar &L_2ipjp1_k_dy,
426 OutputScalar &L_2ipjp1_k_dz,
427 OutputScalar &faceValue_x,
428 OutputScalar &faceValue_y,
429 OutputScalar &faceValue_z)
const
431 outputDiv = L_2ipjp1_k * faceDiv + L_2ipjp1_k_dx * faceValue_x + L_2ipjp1_k_dy * faceValue_y + L_2ipjp1_k_dz * faceValue_z;
434 KOKKOS_INLINE_FUNCTION
435 void operator()(
const TeamMember &teamMember)
const {
436 auto pointOrdinal = teamMember.league_rank();
437 OutputScratchView scratch0, scratch1, scratch2, scratch3;
438 if (fad_size_output_ > 0) {
439 scratch0 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
440 scratch1 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
441 scratch2 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
442 scratch3 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
445 scratch0 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
446 scratch1 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
447 scratch2 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
448 scratch3 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
451 const auto & x = inputPoints_(pointOrdinal,0);
452 const auto & y = inputPoints_(pointOrdinal,1);
453 const auto & z = inputPoints_(pointOrdinal,2);
456 const PointScalar lambda[4] = {1. - x - y - z, x, y, z};
457 const PointScalar lambda_dx[4] = {-1., 1., 0., 0.};
458 const PointScalar lambda_dy[4] = {-1., 0., 1., 0.};
459 const PointScalar lambda_dz[4] = {-1., 0., 0., 1.};
468 auto &scratchP = scratch0;
469 auto &scratchP_2ip1 = scratch1;
471 const ordinal_type max_ij_sum = polyOrder_ - 1;
473 for (ordinal_type faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
475 OutputScalar divWeight;
476 computeFaceDivWeight(divWeight, faceOrdinal, lambda_dx, lambda_dy, lambda_dz);
478 OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;
479 computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, faceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);
481 ordinal_type fieldOrdinal = faceOrdinal * numFaceFunctionsPerFace_;
482 computeFaceLegendre(scratchP, faceOrdinal, lambda);
484 for (
int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)
486 for (
int i=0; i<=ij_sum; i++)
488 computeFaceJacobi(scratchP_2ip1, faceOrdinal, i, lambda);
490 const int j = ij_sum - i;
492 auto & output_x = output_(fieldOrdinal,pointOrdinal,0);
493 auto & output_y = output_(fieldOrdinal,pointOrdinal,1);
494 auto & output_z = output_(fieldOrdinal,pointOrdinal,2);
496 faceFunctionValue(output_x, output_y, output_z, i, j,
497 scratchP, scratchP_2ip1, vectorWeight_x,
498 vectorWeight_y, vectorWeight_z, lambda);
509 auto &scratchP = scratch0;
510 auto &scratchP_2ip1 = scratch1;
511 auto &scratchL_2ipjp1 = scratch2;
513 const ordinal_type min_ijk_sum = 1;
514 const ordinal_type max_ijk_sum = polyOrder_-1;
515 const ordinal_type min_ij_sum = 0;
516 const ordinal_type min_k = 1;
517 const ordinal_type min_j = 0;
518 const ordinal_type min_i = 0;
520 OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;
522 for (
int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
526 ordinal_type interiorFamilyFieldOrdinal =
527 numFaceFunctions_ + interiorFamilyOrdinal - 1;
529 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
531 computeFaceLegendreForInterior(scratchP,
532 interiorFamilyOrdinal - 1, lambda);
533 computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, relatedFaceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);
535 for (
int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
537 for (
int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)
539 for (
int i=min_i; i<=ij_sum-min_j; i++)
541 const ordinal_type j = ij_sum-i;
542 const ordinal_type k = ijk_sum - ij_sum;
545 scratchP_2ip1, interiorFamilyOrdinal - 1, i, lambda);
546 computeInteriorIntegratedJacobi(scratchL_2ipjp1, i, j,
547 interiorFamilyOrdinal - 1,
550 OutputScalar V_x, V_y, V_z;
552 faceFunctionValue(V_x, V_y, V_z, i, j, scratchP,
553 scratchP_2ip1, vectorWeight_x,
554 vectorWeight_y, vectorWeight_z, lambda);
557 output_(interiorFamilyFieldOrdinal, pointOrdinal, 0);
559 output_(interiorFamilyFieldOrdinal, pointOrdinal, 1);
561 output_(interiorFamilyFieldOrdinal, pointOrdinal, 2);
563 output_x = V_x * scratchL_2ipjp1(k);
564 output_y = V_y * scratchL_2ipjp1(k);
565 output_z = V_z * scratchL_2ipjp1(k);
567 interiorFamilyFieldOrdinal +=
581 auto &scratchP = scratch0;
582 auto &scratchP_2ip1 = scratch1;
585 ordinal_type fieldOrdinal = 0;
586 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
588 const int max_ij_sum = polyOrder_ - 1;
589 computeFaceLegendre(scratchP, faceOrdinal, lambda);
590 OutputScalar divWeight;
591 computeFaceDivWeight(divWeight, faceOrdinal, lambda_dx, lambda_dy, lambda_dz);
592 for (
int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)
594 for (
int i=0; i<=ij_sum; i++)
596 const int j = ij_sum - i;
598 computeFaceJacobi(scratchP_2ip1, faceOrdinal, i, lambda);
599 auto &outputValue = output_(fieldOrdinal,pointOrdinal);
600 faceFunctionDiv(outputValue, i, j, scratchP, scratchP_2ip1,
611 auto &scratchL_2ipjp1 = scratch2;
612 auto &scratchP_2ipjp1 = scratch3;
614 const int interiorFieldOrdinalOffset = numFaceFunctions_;
615 for (
int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
619 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
621 computeFaceLegendreForInterior(scratchP,
622 interiorFamilyOrdinal - 1, lambda);
623 OutputScalar divWeight;
624 computeFaceDivWeight(divWeight, relatedFaceOrdinal, lambda_dx, lambda_dy, lambda_dz);
626 OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;
627 computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, relatedFaceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);
629 ordinal_type interiorFieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
631 const ordinal_type min_ijk_sum = 1;
632 const ordinal_type max_ijk_sum = polyOrder_-1;
633 const ordinal_type min_ij_sum = 0;
634 const ordinal_type min_k = 1;
635 const ordinal_type min_j = 0;
636 const ordinal_type min_i = 0;
637 for (
int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
639 for (
int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)
641 for (
int i=min_i; i<=ij_sum-min_j; i++)
643 const ordinal_type j = ij_sum-i;
644 const ordinal_type k = ijk_sum - ij_sum;
646 scratchP_2ip1, interiorFamilyOrdinal - 1, i, lambda);
648 OutputScalar faceDiv;
649 faceFunctionDiv(faceDiv, i, j, scratchP, scratchP_2ip1,
652 OutputScalar faceValue_x, faceValue_y, faceValue_z;
654 faceFunctionValue(faceValue_x, faceValue_y, faceValue_z, i,
655 j, scratchP, scratchP_2ip1,
656 vectorWeight_x, vectorWeight_y,
657 vectorWeight_z, lambda);
658 computeInteriorJacobi(scratchP_2ipjp1, i, j,
659 interiorFamilyOrdinal - 1, lambda);
661 computeInteriorIntegratedJacobi(scratchL_2ipjp1, i, j,
662 interiorFamilyOrdinal - 1,
665 OutputScalar L_2ipjp1_k_dx, L_2ipjp1_k_dy, L_2ipjp1_k_dz;
666 gradInteriorIntegratedJacobi(
667 L_2ipjp1_k_dx, L_2ipjp1_k_dy, L_2ipjp1_k_dz,
668 interiorFamilyOrdinal - 1, j, k, scratchP_2ipjp1,
669 lambda, lambda_dx, lambda_dy, lambda_dz);
671 auto &outputDiv = output_(interiorFieldOrdinal, pointOrdinal);
672 interiorFunctionDiv(outputDiv, scratchL_2ipjp1(k), faceDiv,
673 L_2ipjp1_k_dx, L_2ipjp1_k_dy,
674 L_2ipjp1_k_dz, faceValue_x, faceValue_y,
677 interiorFieldOrdinal += numInteriorFamilies;
696 INTREPID2_TEST_FOR_ABORT(
true,
697 ">>> ERROR: (Intrepid2::Hierarchical_HDIV_TET_Functor) Unsupported differential operator");
707 size_t team_shmem_size (
int team_size)
const
710 size_t shmem_size = 0;
711 if (fad_size_output_ > 0)
712 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
714 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);