54 #ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_PYR_h
55 #define Intrepid2_IntegratedLegendreBasis_HGRAD_PYR_h
57 #include <Kokkos_DynRankView.hpp>
59 #include <Intrepid2_config.h>
69 #include "Teuchos_RCP.hpp"
78 template<
class DeviceType,
class OutputScalar,
class PointScalar,
79 class OutputFieldType,
class InputPointsType>
82 using ExecutionSpace =
typename DeviceType::execution_space;
83 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
84 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
85 using OutputScratchView2D = Kokkos::View<OutputScalar**,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
86 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
88 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
89 using TeamMember =
typename TeamPolicy::member_type;
93 OutputFieldType output_;
94 InputPointsType inputPoints_;
97 bool defineVertexFunctions_;
98 int numFields_, numPoints_;
100 size_t fad_size_output_;
102 static const int numVertices = 5;
103 static const int numMixedEdges = 4;
104 static const int numTriEdges = 4;
105 static const int numEdges = 8;
109 const int edge_start_[numEdges] = {0,1,2,3,0,1,2,3};
110 const int edge_end_[numEdges] = {1,2,3,0,4,4,4,4};
113 static const int numQuadFaces = 1;
114 static const int numTriFaces = 4;
117 const int tri_face_vertex_0[numTriFaces] = {0,1,3,0};
118 const int tri_face_vertex_1[numTriFaces] = {1,2,2,3};
119 const int tri_face_vertex_2[numTriFaces] = {4,4,4,4};
122 int polyOrder,
bool defineVertexFunctions)
123 : opType_(opType), output_(output), inputPoints_(inputPoints),
124 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
125 fad_size_output_(getScalarDimensionForView(output))
127 numFields_ = output.extent_int(0);
128 numPoints_ = output.extent_int(1);
129 const auto & p = polyOrder;
130 const auto p3 = p * p * p;
131 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
132 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != p3 + 3 * p + 1, std::invalid_argument,
"output field size does not match basis cardinality");
135 KOKKOS_INLINE_FUNCTION
136 void operator()(
const TeamMember & teamMember )
const
138 auto pointOrdinal = teamMember.league_rank();
139 OutputScratchView scratch1D_1, scratch1D_2, scratch1D_3;
140 OutputScratchView scratch1D_4, scratch1D_5, scratch1D_6;
141 OutputScratchView scratch1D_7, scratch1D_8, scratch1D_9;
142 OutputScratchView2D scratch2D_1, scratch2D_2, scratch2D_3;
143 const int numAlphaValues = (polyOrder_-1 > 1) ? (polyOrder_-1) : 1;
144 if (fad_size_output_ > 0) {
145 scratch1D_1 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
146 scratch1D_2 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
147 scratch1D_3 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
148 scratch1D_4 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
149 scratch1D_5 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
150 scratch1D_6 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
151 scratch1D_7 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
152 scratch1D_8 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
153 scratch1D_9 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
154 scratch2D_1 = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
155 scratch2D_2 = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
156 scratch2D_3 = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
159 scratch1D_1 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
160 scratch1D_2 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
161 scratch1D_3 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
162 scratch1D_4 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
163 scratch1D_5 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
164 scratch1D_6 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
165 scratch1D_7 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
166 scratch1D_8 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
167 scratch1D_9 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
168 scratch2D_1 = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
169 scratch2D_2 = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
170 scratch2D_3 = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
173 const auto & x = inputPoints_(pointOrdinal,0);
174 const auto & y = inputPoints_(pointOrdinal,1);
175 const auto & z = inputPoints_(pointOrdinal,2);
181 Kokkos::Array<PointScalar,3> coords;
182 transformToESEASPyramid<>(coords[0], coords[1], coords[2], x, y, z);
186 Array<PointScalar,5> lambda;
187 Array<Kokkos::Array<PointScalar,3>,5> lambdaGrad;
189 Array<Array<PointScalar,3>,2> mu;
190 Array<Array<Array<PointScalar,3>,3>,2> muGrad;
192 Array<Array<PointScalar,2>,3> nu;
193 Array<Array<Array<PointScalar,3>,2>,3> nuGrad;
195 affinePyramid(lambda, lambdaGrad, mu, muGrad, nu, nuGrad, coords);
202 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
204 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
206 if (!defineVertexFunctions_)
210 output_(0,pointOrdinal) = 1.0;
214 auto & Li_a1 = scratch1D_1;
215 auto & Li_a2 = scratch1D_2;
217 Polynomials::shiftedScaledIntegratedLegendreValues(Li_a1, polyOrder_, nu[1][0], nu[0][0] + nu[1][0]);
218 Polynomials::shiftedScaledIntegratedLegendreValues(Li_a2, polyOrder_, nu[1][1], nu[0][1] + nu[1][1]);
222 int fieldOrdinalOffset = numVertices;
223 for (
int edgeOrdinal=0; edgeOrdinal<numMixedEdges; edgeOrdinal++)
227 int a = (edgeOrdinal % 2 == 0) ? 1 : 2;
229 auto & Li = (a == 1) ? Li_a1 : Li_a2;
232 int c = ((edgeOrdinal == 0) || (edgeOrdinal == 3)) ? 0 : 1;
233 for (
int i=2; i<=polyOrder_; i++)
235 output_(fieldOrdinalOffset,pointOrdinal) = mu[c][b-1] * Li(i);
236 fieldOrdinalOffset++;
242 auto & Li = scratch1D_1;
243 for (
int edgeOrdinal=0; edgeOrdinal<numMixedEdges; edgeOrdinal++)
245 const auto & lambda_a = lambda[edgeOrdinal];
246 Polynomials::shiftedScaledIntegratedLegendreValues(Li, polyOrder_, lambda[4], lambda_a + lambda[4]);
247 for (
int i=2; i<=polyOrder_; i++)
249 output_(fieldOrdinalOffset,pointOrdinal) = Li(i);
250 fieldOrdinalOffset++;
259 auto & Lj = scratch1D_2;
260 Polynomials::shiftedScaledIntegratedLegendreValues(Li, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
261 Polynomials::shiftedScaledIntegratedLegendreValues(Lj, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
263 for (
int j=2; j<=polyOrder_; j++)
265 for (
int i=2; i<=polyOrder_; i++)
267 output_(fieldOrdinalOffset,pointOrdinal) = mu[0][2] * Li(i) * Lj(j);
268 fieldOrdinalOffset++;
282 for (
int faceOrdinal=0; faceOrdinal<numTriFaces; faceOrdinal++)
286 int a = (faceOrdinal % 2 == 0) ? 1 : 2;
290 int c = ((faceOrdinal == 0) || (faceOrdinal == 3)) ? 0 : 1;
292 const auto & s0 = nu[0][a-1];
293 const auto & s1 = nu[1][a-1];
294 const auto & s2 = nu[2][a-1];
295 const PointScalar jacobiScaling = s0 + s1 + s2;
299 auto & jacobi = scratch2D_1;
300 for (
int n=2; n<=polyOrder_; n++)
302 const double alpha = n*2;
303 const int alphaOrdinal = n-2;
304 using Kokkos::subview;
306 auto jacobi_alpha = subview(jacobi, alphaOrdinal, ALL);
307 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
310 auto & edge_s0s1 = scratch1D_1;
311 Polynomials::shiftedScaledIntegratedLegendreValues(edge_s0s1, polyOrder_, s1, s0 + s1);
313 for (
int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
315 for (
int i=2; i<totalPolyOrder; i++)
317 const auto & edgeValue = edge_s0s1(i);
318 const int alphaOrdinal = i-2;
320 const int j = totalPolyOrder - i;
321 const auto & jacobiValue = jacobi(alphaOrdinal,j);
322 output_(fieldOrdinalOffset,pointOrdinal) = mu[c][b-1] * edgeValue * jacobiValue;
324 fieldOrdinalOffset++;
336 auto & Lk = scratch1D_3;
337 Polynomials::shiftedScaledIntegratedLegendreValues(Li, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
338 Polynomials::shiftedScaledIntegratedLegendreValues(Lj, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
339 Polynomials::shiftedScaledIntegratedLegendreValues(Lk, polyOrder_, mu[1][2], mu[0][2] + mu[1][2]);
341 for (
int k=2; k<=polyOrder_; k++)
343 for (
int j=2; j<=polyOrder_; j++)
345 for (
int i=2; i<=polyOrder_; i++)
347 output_(fieldOrdinalOffset,pointOrdinal) = Lk(k) * Li(i) * Lj(j);
348 fieldOrdinalOffset++;
359 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
361 for (
int d=0; d<3; d++)
363 output_(vertexOrdinal,pointOrdinal,d) = lambdaGrad[vertexOrdinal][d];
367 if (!defineVertexFunctions_)
371 output_(0,pointOrdinal,0) = 0.0;
372 output_(0,pointOrdinal,1) = 0.0;
373 output_(0,pointOrdinal,2) = 0.0;
377 int fieldOrdinalOffset = numVertices;
380 auto & P_i_minus_1 = scratch1D_1;
381 auto & L_i_dt = scratch1D_2;
382 auto & L_i = scratch1D_3;
384 for (
int edgeOrdinal=0; edgeOrdinal<numMixedEdges; edgeOrdinal++)
388 int a = (edgeOrdinal % 2 == 0) ? 1 : 2;
391 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, nu[1][a-1], nu[0][a-1] + nu[1][a-1]);
392 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, nu[1][a-1], nu[0][a-1] + nu[1][a-1]);
393 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, nu[1][a-1], nu[0][a-1] + nu[1][a-1]);
397 int c = ((edgeOrdinal == 0) || (edgeOrdinal == 3)) ? 0 : 1;
398 for (
int i=2; i<=polyOrder_; i++)
401 const auto & R_i_minus_1 = L_i_dt(i);
403 for (
int d=0; d<3; d++)
407 OutputScalar grad_Li_d = P_i_minus_1(i-1) * nuGrad[1][a-1][d] + R_i_minus_1 * (nuGrad[0][a-1][d] + nuGrad[1][a-1][d]);
408 output_(fieldOrdinalOffset,pointOrdinal,d) = muGrad[c][b-1][d] * L_i(i) + mu[c][b-1] * grad_Li_d;
410 fieldOrdinalOffset++;
415 P_i_minus_1 = scratch1D_1;
416 L_i_dt = scratch1D_2;
418 for (
int edgeOrdinal=0; edgeOrdinal<numMixedEdges; edgeOrdinal++)
420 const auto & lambda_a = lambda [edgeOrdinal];
421 const auto & lambdaGrad_a = lambdaGrad[edgeOrdinal];
422 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, lambda[4], lambda_a + lambda[4]);
423 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, lambda[4], lambda_a + lambda[4]);
424 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, lambda[4], lambda_a + lambda[4]);
426 for (
int i=2; i<=polyOrder_; i++)
428 const auto & R_i_minus_1 = L_i_dt(i);
429 for (
int d=0; d<3; d++)
433 OutputScalar grad_Li_d = P_i_minus_1(i-1) * lambdaGrad[4][d] + R_i_minus_1 * (lambdaGrad_a[d] + lambdaGrad[4][d]);
434 output_(fieldOrdinalOffset,pointOrdinal,d) = grad_Li_d;
436 fieldOrdinalOffset++;
442 P_i_minus_1 = scratch1D_1;
443 L_i_dt = scratch1D_2;
445 auto & P_j_minus_1 = scratch1D_4;
446 auto & L_j_dt = scratch1D_5;
447 auto & L_j = scratch1D_6;
448 Polynomials::shiftedScaledIntegratedLegendreValues(L_i, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
449 Polynomials::shiftedScaledIntegratedLegendreValues(L_j, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
451 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, mu[1][0], mu[0][0] + mu[1][0]);
452 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
453 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
455 Polynomials::shiftedScaledLegendreValues (P_j_minus_1, polyOrder_-1, mu[1][1], mu[0][1] + mu[1][1]);
456 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_j_dt, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
457 Polynomials::shiftedScaledIntegratedLegendreValues (L_j, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
460 for (
int j=2; j<=polyOrder_; j++)
462 const auto & R_j_minus_1 = L_j_dt(j);
464 for (
int i=2; i<=polyOrder_; i++)
466 const auto & R_i_minus_1 = L_i_dt(i);
468 OutputScalar phi_quad = L_i(i) * L_j(j);
470 for (
int d=0; d<3; d++)
474 OutputScalar grad_Lj_d = P_j_minus_1(j-1) * muGrad[1][1][d] + R_j_minus_1 * (muGrad[0][1][d] + muGrad[1][1][d]);
476 OutputScalar grad_Li_d = P_i_minus_1(i-1) * muGrad[1][0][d] + R_i_minus_1 * (muGrad[0][0][d] + muGrad[1][0][d]);
478 OutputScalar grad_phi_quad_d = L_i(i) * grad_Lj_d + L_j(j) * grad_Li_d;
480 output_(fieldOrdinalOffset,pointOrdinal,d) = mu[0][2] * grad_phi_quad_d + phi_quad * muGrad[0][2][d];
483 fieldOrdinalOffset++;
488 for (
int faceOrdinal=0; faceOrdinal<numTriFaces; faceOrdinal++)
492 int a = (faceOrdinal % 2 == 0) ? 1 : 2;
496 int c = ((faceOrdinal == 0) || (faceOrdinal == 3)) ? 0 : 1;
498 const auto & s0 = nu[0][a-1];
499 const auto & s1 = nu[1][a-1];
500 const auto & s2 = nu[2][a-1];
502 const auto & s0Grad = nuGrad[0][a-1];
503 const auto & s1Grad = nuGrad[1][a-1];
504 const auto & s2Grad = nuGrad[2][a-1];
506 const PointScalar jacobiScaling = s0 + s1 + s2;
511 P_i_minus_1 = scratch1D_1;
512 L_i_dt = scratch1D_2;
515 auto & L_2i_j_dt = scratch2D_1;
516 auto & L_2i_j = scratch2D_2;
517 auto & P_2i_j_minus_1 = scratch2D_3;
518 for (
int n=2; n<=polyOrder_; n++)
520 const double alpha = n*2;
521 const int alphaOrdinal = n-2;
522 using Kokkos::subview;
524 auto L_2i_j_dt_alpha = subview(L_2i_j_dt, alphaOrdinal, ALL);
525 auto L_2i_j_alpha = subview(L_2i_j, alphaOrdinal, ALL);
526 auto P_2i_j_minus_1_alpha = subview(P_2i_j_minus_1, alphaOrdinal, ALL);
527 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2i_j_dt_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
528 Polynomials::shiftedScaledIntegratedJacobiValues ( L_2i_j_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
529 Polynomials::shiftedScaledJacobiValues (P_2i_j_minus_1_alpha, alpha, polyOrder_-1, s2, jacobiScaling);
531 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, s1, s0 + s1);
532 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, s1, s0 + s1);
533 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, s1, s0 + s1);
535 for (
int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
537 for (
int i=2; i<totalPolyOrder; i++)
539 const int alphaOrdinal = i-2;
540 const int j = totalPolyOrder - i;
542 const auto & R_i_minus_1 = L_i_dt(i);
543 OutputScalar phi_tri = L_2i_j(alphaOrdinal,j) * L_i(i);
545 for (
int d=0; d<3; d++)
548 OutputScalar grad_Li_d = P_i_minus_1(i-1) * s1Grad[d] + R_i_minus_1 * (s0Grad[d] + s1Grad[d]);
549 OutputScalar grad_L2i_j_d = P_2i_j_minus_1(alphaOrdinal,j-1) * s2Grad[d] + L_2i_j_dt(alphaOrdinal,j) * (s0Grad[d] + s1Grad[d] + s2Grad[d]);
550 OutputScalar grad_phi_tri_d = L_i(i) * grad_L2i_j_d + L_2i_j(alphaOrdinal,j) * grad_Li_d;
552 output_(fieldOrdinalOffset,pointOrdinal,d) = mu[c][b-1] * grad_phi_tri_d + phi_tri * muGrad[c][b-1][d];
554 fieldOrdinalOffset++;
560 P_i_minus_1 = scratch1D_1;
561 L_i_dt = scratch1D_2;
563 P_j_minus_1 = scratch1D_4;
564 L_j_dt = scratch1D_5;
566 auto & P_k_minus_1 = scratch1D_7;
567 auto & L_k_dt = scratch1D_8;
568 auto & L_k = scratch1D_9;
570 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, mu[1][0], mu[0][0] + mu[1][0]);
571 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
572 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, mu[1][0], mu[0][0] + mu[1][0]);
574 Polynomials::shiftedScaledLegendreValues (P_j_minus_1, polyOrder_-1, mu[1][1], mu[0][1] + mu[1][1]);
575 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_j_dt, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
576 Polynomials::shiftedScaledIntegratedLegendreValues (L_j, polyOrder_, mu[1][1], mu[0][1] + mu[1][1]);
578 Polynomials::shiftedScaledLegendreValues (P_k_minus_1, polyOrder_-1, mu[1][2], mu[0][2] + mu[1][2]);
579 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_k_dt, polyOrder_, mu[1][2], mu[0][2] + mu[1][2]);
580 Polynomials::shiftedScaledIntegratedLegendreValues (L_k, polyOrder_, mu[1][2], mu[0][2] + mu[1][2]);
583 for (
int k=2; k<=polyOrder_; k++)
585 const auto & R_k_minus_1 = L_k_dt(k);
587 for (
int j=2; j<=polyOrder_; j++)
589 const auto & R_j_minus_1 = L_j_dt(j);
591 for (
int i=2; i<=polyOrder_; i++)
593 const auto & R_i_minus_1 = L_i_dt(i);
595 OutputScalar phi_quad = L_i(i) * L_j(j);
597 for (
int d=0; d<3; d++)
601 OutputScalar grad_Li_d = P_i_minus_1(i-1) * muGrad[1][0][d] + R_i_minus_1 * (muGrad[0][0][d] + muGrad[1][0][d]);
603 OutputScalar grad_Lj_d = P_j_minus_1(j-1) * muGrad[1][1][d] + R_j_minus_1 * (muGrad[0][1][d] + muGrad[1][1][d]);
605 OutputScalar grad_Lk_d = P_k_minus_1(k-1) * muGrad[1][2][d] + R_k_minus_1 * (muGrad[0][2][d] + muGrad[1][2][d]);
607 OutputScalar grad_phi_quad_d = L_i(i) * grad_Lj_d + L_j(j) * grad_Li_d;
609 output_(fieldOrdinalOffset,pointOrdinal,d) = L_k(k) * grad_phi_quad_d + phi_quad * grad_Lk_d;
612 fieldOrdinalOffset++;
617 for (
int basisOrdinal=0; basisOrdinal<numFields_; basisOrdinal++)
620 const auto dx_eseas = output_(basisOrdinal,pointOrdinal,0);
621 const auto dy_eseas = output_(basisOrdinal,pointOrdinal,1);
622 const auto dz_eseas = output_(basisOrdinal,pointOrdinal,2);
624 auto &dx_int2 = output_(basisOrdinal,pointOrdinal,0);
625 auto &dy_int2 = output_(basisOrdinal,pointOrdinal,1);
626 auto &dz_int2 = output_(basisOrdinal,pointOrdinal,2);
628 transformFromESEASPyramidGradient(dx_int2, dy_int2, dz_int2, dx_eseas, dy_eseas, dz_eseas);
642 INTREPID2_TEST_FOR_ABORT(
true,
643 ">>> ERROR: (Intrepid2::Hierarchical_HGRAD_PYR_Functor) Computing of second and higher-order derivatives is not currently supported");
646 device_assert(
false);
653 size_t team_shmem_size (
int team_size)
const
660 const int numAlphaValues = std::max(polyOrder_-1, 1);
661 size_t shmem_size = 0;
662 if (fad_size_output_ > 0)
665 shmem_size += 9 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
667 shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1, fad_size_output_);
672 shmem_size += 9 * OutputScratchView::shmem_size(polyOrder_ + 1);
674 shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1);
698 template<
typename DeviceType,
699 typename OutputScalar = double,
700 typename PointScalar = double,
701 bool defineVertexFunctions =
true>
703 :
public Basis<DeviceType,OutputScalar,PointScalar>
719 EPointType pointType_;
733 polyOrder_(polyOrder),
734 pointType_(pointType)
736 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,
"PointType not supported");
739 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<> >() );
744 const int degreeLength = 1;
748 int fieldOrdinalOffset = 0;
751 for (
int i=0; i<numVertices; i++)
758 if (!defineVertexFunctions)
763 fieldOrdinalOffset += numVertices;
766 const int numFunctionsPerEdge = polyOrder - 1;
768 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
770 for (
int i=0; i<numFunctionsPerEdge; i++)
775 fieldOrdinalOffset += numFunctionsPerEdge;
780 const int numFunctionsPerQuadFace = (polyOrder-1)*(polyOrder-1);
783 for (
int j=2; j<=polyOrder_; j++)
785 for (
int i=2; i<=polyOrder_; i++)
789 fieldOrdinalOffset++;
793 const int numFunctionsPerTriFace = ((polyOrder-1)*(polyOrder-2))/2;
794 const int numTriFaces = 4;
795 for (
int faceOrdinal=0; faceOrdinal<numTriFaces; faceOrdinal++)
797 for (
int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
799 const int totalFaceDofs = (totalPolyOrder-2)*(totalPolyOrder-1)/2;
800 const int totalFaceDofsPrevious = (totalPolyOrder-3)*(totalPolyOrder-2)/2;
801 const int faceDofsForPolyOrder = totalFaceDofs - totalFaceDofsPrevious;
802 for (
int i=0; i<faceDofsForPolyOrder; i++)
806 fieldOrdinalOffset++;
812 const int numFunctionsPerVolume = (polyOrder-1)*(polyOrder-1)*(polyOrder-1);
813 const int numVolumes = 1;
814 for (
int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
817 for (
int k=2; k<=polyOrder_; k++)
819 for (
int j=2; j<=polyOrder_; j++)
821 for (
int i=2; i<=polyOrder_; i++)
823 const int max_ij = std::max(i,j);
824 const int max_ijk = std::max(max_ij,k);
827 fieldOrdinalOffset++;
833 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
859 const int intrepid2FaceOrdinals[5] {4,0,1,2,3};
864 const ordinal_type tagSize = 4;
865 const ordinal_type posScDim = 0;
866 const ordinal_type posScOrd = 1;
867 const ordinal_type posDfOrd = 2;
870 const int vertexDim = 0, edgeDim = 1, faceDim = 2, volumeDim = 3;
872 if (defineVertexFunctions) {
875 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
877 tagView(tagNumber*tagSize+0) = vertexDim;
878 tagView(tagNumber*tagSize+1) = vertexOrdinal;
879 tagView(tagNumber*tagSize+2) = 0;
880 tagView(tagNumber*tagSize+3) = 1;
883 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
885 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
887 tagView(tagNumber*tagSize+0) = edgeDim;
888 tagView(tagNumber*tagSize+1) = edgeOrdinal;
889 tagView(tagNumber*tagSize+2) = functionOrdinal;
890 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge;
896 const int faceOrdinalESEAS = 0;
897 const int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
899 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerQuadFace; functionOrdinal++)
901 tagView(tagNumber*tagSize+0) = faceDim;
902 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2;
903 tagView(tagNumber*tagSize+2) = functionOrdinal;
904 tagView(tagNumber*tagSize+3) = numFunctionsPerQuadFace;
908 for (
int triFaceOrdinalESEAS=0; triFaceOrdinalESEAS<numTriFaces; triFaceOrdinalESEAS++)
910 const int faceOrdinalESEAS = triFaceOrdinalESEAS + 1;
911 const int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
912 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerTriFace; functionOrdinal++)
914 tagView(tagNumber*tagSize+0) = faceDim;
915 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2;
916 tagView(tagNumber*tagSize+2) = functionOrdinal;
917 tagView(tagNumber*tagSize+3) = numFunctionsPerTriFace;
921 for (
int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
923 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerVolume; functionOrdinal++)
925 tagView(tagNumber*tagSize+0) = volumeDim;
926 tagView(tagNumber*tagSize+1) = volumeOrdinal;
927 tagView(tagNumber*tagSize+2) = functionOrdinal;
928 tagView(tagNumber*tagSize+3) = numFunctionsPerVolume;
932 INTREPID2_TEST_FOR_EXCEPTION(tagNumber != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis tag enumeration is incorrect");
935 for (ordinal_type i=0;i<cardinality;++i) {
936 tagView(i*tagSize+0) = volumeDim;
937 tagView(i*tagSize+1) = 0;
938 tagView(i*tagSize+2) = i;
939 tagView(i*tagSize+3) = cardinality;
961 return "Intrepid2_IntegratedLegendreBasis_HGRAD_PYR";
994 const EOperator operatorType = OPERATOR_VALUE )
const override
996 auto numPoints = inputPoints.extent_int(0);
1000 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
1002 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
1003 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
1004 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
1005 const int teamSize = 1;
1007 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
1008 Kokkos::parallel_for(
"Hierarchical_HGRAD_PYR_Functor", policy, functor);
1019 BasisPtr<DeviceType,OutputScalar,PointScalar>
1027 return Teuchos::rcp(
new HGRAD_LINE(p));
1029 else if (subCellDim == 2)
1031 if (subCellOrd == 4)
1033 return Teuchos::rcp(
new HGRAD_QUAD(p));
1037 return Teuchos::rcp(
new HGRAD_TRI(p));
1040 INTREPID2_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Input parameters out of bounds");
1047 virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>
1049 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
1051 return Teuchos::rcp(
new HostBasisType(polyOrder_, pointType_) );
H(grad) basis on the line based on integrated Legendre polynomials.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
const char * getName() const override
Returns basis name.
virtual void getValues(const ExecutionSpace &, OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
EBasis basisType_
Type of the basis.
Defines several coordinates and their gradients on the pyramid; maps from Intrepid2 (shards) pyramid ...
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
Free functions, callable from device code, that implement various polynomials useful in basis definit...
EFunctionSpace functionSpace_
The function space in which the basis is defined.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
Header function for Intrepid2::Util class and other utility functions.
Implementation of H(grad) basis on the quadrilateral that is templated on H(grad) on the line...
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
virtual BasisPtr< typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
IntegratedLegendreBasis_HGRAD_PYR(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line: e...
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
H(grad) basis on the triangle based on integrated Legendre polynomials.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
ordinal_type getDegree() const
Returns the degree of the basis.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
OrdinalTypeArray2DHost fieldOrdinalH1PolynomialDegree_
H^1 polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
Functor for computing values for the IntegratedLegendreBasis_HGRAD_PYR class.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
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)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Header file for the abstract base class Intrepid2::Basis.
virtual bool requireOrientation() const override
True if orientation is required.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.