49 #ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h
50 #define Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h
52 #include <Kokkos_DynRankView.hpp>
54 #include <Intrepid2_config.h>
69 template<
class DeviceType,
class OutputScalar,
class PointScalar,
70 class OutputFieldType,
class InputPointsType>
73 using ExecutionSpace =
typename DeviceType::execution_space;
74 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
75 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
76 using OutputScratchView2D = Kokkos::View<OutputScalar**,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
77 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
79 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
80 using TeamMember =
typename TeamPolicy::member_type;
84 OutputFieldType output_;
85 InputPointsType inputPoints_;
88 bool defineVertexFunctions_;
89 int numFields_, numPoints_;
91 size_t fad_size_output_;
93 static const int numVertices = 4;
94 static const int numEdges = 6;
96 const int edge_start_[numEdges] = {0,1,0,0,1,2};
97 const int edge_end_[numEdges] = {1,2,2,3,3,3};
99 static const int numFaces = 4;
100 const int face_vertex_0[numFaces] = {0,0,1,0};
101 const int face_vertex_1[numFaces] = {1,1,2,2};
102 const int face_vertex_2[numFaces] = {2,3,3,3};
106 const int face_ordinal_of_first_edge[numFaces] = {0,0,1,2};
109 int polyOrder,
bool defineVertexFunctions)
110 : opType_(opType), output_(output), inputPoints_(inputPoints),
111 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
112 fad_size_output_(getScalarDimensionForView(output))
114 numFields_ = output.extent_int(0);
115 numPoints_ = output.extent_int(1);
116 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
117 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)*(polyOrder_+3)/6, std::invalid_argument,
"output field size does not match basis cardinality");
120 KOKKOS_INLINE_FUNCTION
121 void operator()(
const TeamMember & teamMember )
const
123 const int numFaceBasisFunctionsPerFace = (polyOrder_-2) * (polyOrder_-1) / 2;
124 const int numInteriorBasisFunctions = (polyOrder_-3) * (polyOrder_-2) * (polyOrder_-1) / 6;
126 auto pointOrdinal = teamMember.league_rank();
127 OutputScratchView legendre_values1_at_point, legendre_values2_at_point;
128 OutputScratchView2D jacobi_values1_at_point, jacobi_values2_at_point, jacobi_values3_at_point;
129 const int numAlphaValues = (polyOrder_-1 > 1) ? (polyOrder_-1) : 1;
130 if (fad_size_output_ > 0) {
131 legendre_values1_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
132 legendre_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
133 jacobi_values1_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
134 jacobi_values2_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
135 jacobi_values3_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
138 legendre_values1_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
139 legendre_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
140 jacobi_values1_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
141 jacobi_values2_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
142 jacobi_values3_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
145 const auto & x = inputPoints_(pointOrdinal,0);
146 const auto & y = inputPoints_(pointOrdinal,1);
147 const auto & z = inputPoints_(pointOrdinal,2);
150 const PointScalar lambda[numVertices] = {1. - x - y - z, x, y, z};
151 const PointScalar lambda_dx[numVertices] = {-1., 1., 0., 0.};
152 const PointScalar lambda_dy[numVertices] = {-1., 0., 1., 0.};
153 const PointScalar lambda_dz[numVertices] = {-1., 0., 0., 1.};
155 const int num1DEdgeFunctions = polyOrder_ - 1;
162 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
164 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
166 if (!defineVertexFunctions_)
170 output_(0,pointOrdinal) = 1.0;
174 int fieldOrdinalOffset = numVertices;
175 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
177 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
178 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
180 Polynomials::shiftedScaledIntegratedLegendreValues(legendre_values1_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
181 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
184 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = legendre_values1_at_point(edgeFunctionOrdinal+2);
186 fieldOrdinalOffset += num1DEdgeFunctions;
192 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
194 const auto & s0 = lambda[face_vertex_0[faceOrdinal]];
195 const auto & s1 = lambda[face_vertex_1[faceOrdinal]];
196 const auto & s2 = lambda[face_vertex_2[faceOrdinal]];
197 const PointScalar jacobiScaling = s0 + s1 + s2;
200 for (
int n=2; n<=polyOrder_; n++)
202 const double alpha = n*2;
203 const int alphaOrdinal = n-2;
204 using Kokkos::subview;
206 auto jacobi_alpha = subview(jacobi_values1_at_point, alphaOrdinal, ALL);
207 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
210 const int edgeOrdinal = face_ordinal_of_first_edge[faceOrdinal];
211 int localFaceBasisOrdinal = 0;
212 for (
int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
214 for (
int i=2; i<totalPolyOrder; i++)
216 const int edgeBasisOrdinal = edgeOrdinal*num1DEdgeFunctions + i-2 + numVertices;
217 const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
218 const int alphaOrdinal = i-2;
220 const int j = totalPolyOrder - i;
221 const auto & jacobiValue = jacobi_values1_at_point(alphaOrdinal,j);
222 const int fieldOrdinal = fieldOrdinalOffset + localFaceBasisOrdinal;
223 output_(fieldOrdinal,pointOrdinal) = edgeValue * jacobiValue;
225 localFaceBasisOrdinal++;
228 fieldOrdinalOffset += numFaceBasisFunctionsPerFace;
232 for (
int n=3; n<=polyOrder_; n++)
234 const double alpha = n*2;
235 const double jacobiScaling = 1.0;
236 const int alphaOrdinal = n-3;
237 using Kokkos::subview;
239 auto jacobi_alpha = subview(jacobi_values1_at_point, alphaOrdinal, ALL);
240 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-3, lambda[3], jacobiScaling);
245 const int min_ij = min_i + min_j;
246 const int min_ijk = min_ij + min_k;
247 int localInteriorBasisOrdinal = 0;
248 for (
int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= polyOrder_; totalPolyOrder_ijk++)
250 int localFaceBasisOrdinal = 0;
251 for (
int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
253 for (
int i=2; i <= totalPolyOrder_ij-min_j; i++)
255 const int j = totalPolyOrder_ij - i;
256 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
257 const int faceBasisOrdinal = numEdges*num1DEdgeFunctions + numVertices + localFaceBasisOrdinal;
258 const auto & faceValue = output_(faceBasisOrdinal,pointOrdinal);
259 const int alphaOrdinal = (i+j)-3;
260 localFaceBasisOrdinal++;
262 const int fieldOrdinal = fieldOrdinalOffset + localInteriorBasisOrdinal;
263 const auto & jacobiValue = jacobi_values1_at_point(alphaOrdinal,k);
264 output_(fieldOrdinal,pointOrdinal) = faceValue * jacobiValue;
265 localInteriorBasisOrdinal++;
269 fieldOrdinalOffset += numInteriorBasisFunctions;
276 if (defineVertexFunctions_)
280 output_(0,pointOrdinal,0) = -1.0;
281 output_(0,pointOrdinal,1) = -1.0;
282 output_(0,pointOrdinal,2) = -1.0;
288 output_(0,pointOrdinal,0) = 0.0;
289 output_(0,pointOrdinal,1) = 0.0;
290 output_(0,pointOrdinal,2) = 0.0;
293 output_(1,pointOrdinal,0) = 1.0;
294 output_(1,pointOrdinal,1) = 0.0;
295 output_(1,pointOrdinal,2) = 0.0;
297 output_(2,pointOrdinal,0) = 0.0;
298 output_(2,pointOrdinal,1) = 1.0;
299 output_(2,pointOrdinal,2) = 0.0;
301 output_(3,pointOrdinal,0) = 0.0;
302 output_(3,pointOrdinal,1) = 0.0;
303 output_(3,pointOrdinal,2) = 1.0;
306 int fieldOrdinalOffset = numVertices;
318 auto & P_i_minus_1 = legendre_values1_at_point;
319 auto & L_i_dt = legendre_values2_at_point;
320 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
322 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
323 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
325 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
326 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
327 const auto & s0_dz = lambda_dz[edge_start_[edgeOrdinal]];
328 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
329 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
330 const auto & s1_dz = lambda_dz[ edge_end_[edgeOrdinal]];
332 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
333 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
334 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
337 const int i = edgeFunctionOrdinal+2;
338 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
339 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
340 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,2) = P_i_minus_1(i-1) * s1_dz + L_i_dt(i) * (s1_dz + s0_dz);
342 fieldOrdinalOffset += num1DEdgeFunctions;
363 auto & L_i = legendre_values2_at_point;
364 auto & L_2i_j_dt = jacobi_values1_at_point;
365 auto & L_2i_j = jacobi_values2_at_point;
366 auto & P_2i_j_minus_1 = jacobi_values3_at_point;
368 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
370 const auto & s0 = lambda[face_vertex_0[faceOrdinal]];
371 const auto & s1 = lambda[face_vertex_1[faceOrdinal]];
372 const auto & s2 = lambda[face_vertex_2[faceOrdinal]];
373 Polynomials::shiftedScaledIntegratedLegendreValues(L_i, polyOrder_, s1, s0+s1);
375 const PointScalar jacobiScaling = s0 + s1 + s2;
378 for (
int n=2; n<=polyOrder_; n++)
380 const double alpha = n*2;
381 const int alphaOrdinal = n-2;
382 using Kokkos::subview;
384 auto L_2i_j_dt_alpha = subview(L_2i_j_dt, alphaOrdinal, ALL);
385 auto L_2i_j_alpha = subview(L_2i_j, alphaOrdinal, ALL);
386 auto P_2i_j_minus_1_alpha = subview(P_2i_j_minus_1, alphaOrdinal, ALL);
387 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2i_j_dt_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
388 Polynomials::shiftedScaledIntegratedJacobiValues (L_2i_j_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
389 Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1_alpha, alpha, polyOrder_-1, s2, jacobiScaling);
392 const int edgeOrdinal = face_ordinal_of_first_edge[faceOrdinal];
393 int localFaceOrdinal = 0;
394 for (
int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
396 for (
int i=2; i<totalPolyOrder; i++)
398 const int edgeBasisOrdinal = edgeOrdinal*num1DEdgeFunctions + i-2 + numVertices;
399 const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
400 const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
401 const auto & grad_L_i_dz = output_(edgeBasisOrdinal,pointOrdinal,2);
403 const int alphaOrdinal = i-2;
405 const auto & s0_dx = lambda_dx[face_vertex_0[faceOrdinal]];
406 const auto & s0_dy = lambda_dy[face_vertex_0[faceOrdinal]];
407 const auto & s0_dz = lambda_dz[face_vertex_0[faceOrdinal]];
408 const auto & s1_dx = lambda_dx[face_vertex_1[faceOrdinal]];
409 const auto & s1_dy = lambda_dy[face_vertex_1[faceOrdinal]];
410 const auto & s1_dz = lambda_dz[face_vertex_1[faceOrdinal]];
411 const auto & s2_dx = lambda_dx[face_vertex_2[faceOrdinal]];
412 const auto & s2_dy = lambda_dy[face_vertex_2[faceOrdinal]];
413 const auto & s2_dz = lambda_dz[face_vertex_2[faceOrdinal]];
415 int j = totalPolyOrder - i;
418 auto & l_2i_j = L_2i_j(alphaOrdinal,j);
420 auto & l_2i_j_dt = L_2i_j_dt(alphaOrdinal,j);
421 auto & p_2i_j_minus_1 = P_2i_j_minus_1(alphaOrdinal,j-1);
423 const OutputScalar basisValue_dx = l_2i_j * grad_L_i_dx + l_i * (p_2i_j_minus_1 * s2_dx + l_2i_j_dt * (s0_dx + s1_dx + s2_dx));
424 const OutputScalar basisValue_dy = l_2i_j * grad_L_i_dy + l_i * (p_2i_j_minus_1 * s2_dy + l_2i_j_dt * (s0_dy + s1_dy + s2_dy));
425 const OutputScalar basisValue_dz = l_2i_j * grad_L_i_dz + l_i * (p_2i_j_minus_1 * s2_dz + l_2i_j_dt * (s0_dz + s1_dz + s2_dz));
427 const int fieldOrdinal = fieldOrdinalOffset + localFaceOrdinal;
429 output_(fieldOrdinal,pointOrdinal,0) = basisValue_dx;
430 output_(fieldOrdinal,pointOrdinal,1) = basisValue_dy;
431 output_(fieldOrdinal,pointOrdinal,2) = basisValue_dz;
436 fieldOrdinalOffset += numFaceBasisFunctionsPerFace;
455 auto & L_alpha = jacobi_values1_at_point;
456 auto & P_alpha = jacobi_values2_at_point;
460 const auto & s0 = lambda[0];
461 const auto & s1 = lambda[1];
462 const auto & s2 = lambda[2];
464 Polynomials::shiftedScaledIntegratedLegendreValues(legendre_values1_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
467 const PointScalar jacobiScaling = s0 + s1 + s2;
468 for (
int n=2; n<=polyOrder_; n++)
470 const double alpha = n*2;
471 const int alphaOrdinal = n-2;
472 using Kokkos::subview;
474 auto jacobi_alpha = subview(jacobi_values3_at_point, alphaOrdinal, ALL);
475 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
480 for (
int n=3; n<=polyOrder_; n++)
482 const double alpha = n*2;
483 const double jacobiScaling = 1.0;
484 const int alphaOrdinal = n-3;
485 using Kokkos::subview;
489 auto L = subview(L_alpha, alphaOrdinal, ALL);
490 auto P = subview(P_alpha, alphaOrdinal, ALL);
491 Polynomials::shiftedScaledIntegratedJacobiValues(L, alpha, polyOrder_-3, lambda[3], jacobiScaling);
492 Polynomials::shiftedScaledJacobiValues (P, alpha, polyOrder_-3, lambda[3], jacobiScaling);
498 const int min_ij = min_i + min_j;
499 const int min_ijk = min_ij + min_k;
500 int localInteriorBasisOrdinal = 0;
501 for (
int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= polyOrder_; totalPolyOrder_ijk++)
503 int localFaceBasisOrdinal = 0;
504 for (
int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
506 for (
int i=2; i <= totalPolyOrder_ij-min_j; i++)
508 const int j = totalPolyOrder_ij - i;
509 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
511 const int faceBasisOrdinal = numEdges*num1DEdgeFunctions + numVertices + localFaceBasisOrdinal;
513 const auto & faceValue_dx = output_(faceBasisOrdinal,pointOrdinal,0);
514 const auto & faceValue_dy = output_(faceBasisOrdinal,pointOrdinal,1);
515 const auto & faceValue_dz = output_(faceBasisOrdinal,pointOrdinal,2);
518 OutputScalar faceValue;
520 const auto & edgeValue = legendre_values1_at_point(i);
521 const int alphaOrdinal = i-2;
522 const auto & jacobiValue = jacobi_values3_at_point(alphaOrdinal,j);
523 faceValue = edgeValue * jacobiValue;
525 localFaceBasisOrdinal++;
527 const int alphaOrdinal = (i+j)-3;
529 const int fieldOrdinal = fieldOrdinalOffset + localInteriorBasisOrdinal;
530 const auto & integratedJacobiValue = L_alpha(alphaOrdinal,k);
531 const auto & jacobiValue = P_alpha(alphaOrdinal,k-1);
532 output_(fieldOrdinal,pointOrdinal,0) = integratedJacobiValue * faceValue_dx + faceValue * jacobiValue * lambda_dx[3];
533 output_(fieldOrdinal,pointOrdinal,1) = integratedJacobiValue * faceValue_dy + faceValue * jacobiValue * lambda_dy[3];
534 output_(fieldOrdinal,pointOrdinal,2) = integratedJacobiValue * faceValue_dz + faceValue * jacobiValue * lambda_dz[3];
536 localInteriorBasisOrdinal++;
540 fieldOrdinalOffset += numInteriorBasisFunctions;
552 INTREPID2_TEST_FOR_ABORT(
true,
553 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
556 device_assert(
false);
563 size_t team_shmem_size (
int team_size)
const
570 const int numAlphaValues = std::max(polyOrder_-1, 1);
571 size_t shmem_size = 0;
572 if (fad_size_output_ > 0)
575 shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
577 shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1, fad_size_output_);
582 shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1);
584 shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1);
608 template<
typename DeviceType,
609 typename OutputScalar = double,
610 typename PointScalar = double,
611 bool defineVertexFunctions =
true>
613 :
public Basis<DeviceType,OutputScalar,PointScalar>
629 EPointType pointType_;
643 polyOrder_(polyOrder),
644 pointType_(pointType)
646 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,
"PointType not supported");
649 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
654 const int degreeLength = 1;
658 int fieldOrdinalOffset = 0;
661 const int numFunctionsPerVertex = 1;
662 const int numVertexFunctions = numVertices * numFunctionsPerVertex;
663 for (
int i=0; i<numVertexFunctions; i++)
670 if (!defineVertexFunctions)
675 fieldOrdinalOffset += numVertexFunctions;
678 const int numFunctionsPerEdge = polyOrder - 1;
680 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
682 for (
int i=0; i<numFunctionsPerEdge; i++)
687 fieldOrdinalOffset += numFunctionsPerEdge;
691 const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
692 const int numFaces = 4;
693 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
695 for (
int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
697 const int totalFaceDofs = (totalPolyOrder-2)*(totalPolyOrder-1)/2;
698 const int totalFaceDofsPrevious = (totalPolyOrder-3)*(totalPolyOrder-2)/2;
699 const int faceDofsForPolyOrder = totalFaceDofs - totalFaceDofsPrevious;
700 for (
int i=0; i<faceDofsForPolyOrder; i++)
704 fieldOrdinalOffset++;
710 const int numFunctionsPerVolume = ((polyOrder-1)*(polyOrder-2)*(polyOrder-3))/6;
711 const int numVolumes = 1;
712 for (
int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
714 for (
int totalPolyOrder=4; totalPolyOrder<=polyOrder_; totalPolyOrder++)
716 const int totalInteriorDofs = (totalPolyOrder-3)*(totalPolyOrder-2)*(totalPolyOrder-1)/6;
717 const int totalInteriorDofsPrevious = (totalPolyOrder-4)*(totalPolyOrder-3)*(totalPolyOrder-2)/6;
718 const int interiorDofsForPolyOrder = totalInteriorDofs - totalInteriorDofsPrevious;
720 for (
int i=0; i<interiorDofsForPolyOrder; i++)
724 fieldOrdinalOffset++;
729 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
736 const int intrepid2FaceOrdinals[4] {3,0,1,2};
741 const ordinal_type tagSize = 4;
742 const ordinal_type posScDim = 0;
743 const ordinal_type posScOrd = 1;
744 const ordinal_type posDfOrd = 2;
747 const int vertexDim = 0, edgeDim = 1, faceDim = 2, volumeDim = 3;
749 if (defineVertexFunctions) {
752 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
754 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerVertex; functionOrdinal++)
756 tagView(tagNumber*tagSize+0) = vertexDim;
757 tagView(tagNumber*tagSize+1) = vertexOrdinal;
758 tagView(tagNumber*tagSize+2) = functionOrdinal;
759 tagView(tagNumber*tagSize+3) = numFunctionsPerVertex;
763 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
765 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
767 tagView(tagNumber*tagSize+0) = edgeDim;
768 tagView(tagNumber*tagSize+1) = edgeOrdinal;
769 tagView(tagNumber*tagSize+2) = functionOrdinal;
770 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge;
774 for (
int faceOrdinalESEAS=0; faceOrdinalESEAS<numFaces; faceOrdinalESEAS++)
776 int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
777 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
779 tagView(tagNumber*tagSize+0) = faceDim;
780 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2;
781 tagView(tagNumber*tagSize+2) = functionOrdinal;
782 tagView(tagNumber*tagSize+3) = numFunctionsPerFace;
786 for (
int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
788 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerVolume; functionOrdinal++)
790 tagView(tagNumber*tagSize+0) = volumeDim;
791 tagView(tagNumber*tagSize+1) = volumeOrdinal;
792 tagView(tagNumber*tagSize+2) = functionOrdinal;
793 tagView(tagNumber*tagSize+3) = numFunctionsPerVolume;
797 INTREPID2_TEST_FOR_EXCEPTION(tagNumber != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis tag enumeration is incorrect");
800 for (ordinal_type i=0;i<cardinality;++i) {
801 tagView(i*tagSize+0) = volumeDim;
802 tagView(i*tagSize+1) = 0;
803 tagView(i*tagSize+2) = i;
804 tagView(i*tagSize+3) = cardinality;
826 return "Intrepid2_IntegratedLegendreBasis_HGRAD_TET";
859 const EOperator operatorType = OPERATOR_VALUE )
const override
861 auto numPoints = inputPoints.extent_int(0);
865 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
867 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
868 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
869 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
870 const int teamSize = 1;
872 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
873 Kokkos::parallel_for(
"Hierarchical_HGRAD_TET_Functor", policy, functor);
884 BasisPtr<DeviceType,OutputScalar,PointScalar>
886 if(subCellDim == 1) {
887 return Teuchos::rcp(
new
890 }
else if(subCellDim == 2) {
891 return Teuchos::rcp(
new
895 INTREPID2_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Input parameters out of bounds");
902 virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>
904 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
906 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.
virtual void getValues(const ExecutionSpace &, OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
EBasis basisType_
Type of the basis.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
IntegratedLegendreBasis_HGRAD_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
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...
Free functions, callable from device code, that implement various polynomials useful in basis definit...
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
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.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
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...
Functor for computing values for the IntegratedLegendreBasis_HGRAD_TET class.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
virtual bool requireOrientation() const override
True if orientation is required.
H(grad) basis on the triangle based on integrated Legendre polynomials.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
ordinal_type getDegree() const
Returns the degree of the basis.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
const char * getName() const override
Returns basis name.
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...
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
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.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Header file for the abstract base class Intrepid2::Basis.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.