48 #ifndef Intrepid2_HierarchicalBasis_HCURL_TET_h
49 #define Intrepid2_HierarchicalBasis_HCURL_TET_h
51 #include <Intrepid2_config.h>
53 #include <Kokkos_DynRankView.hpp>
68 template<
class DeviceType,
class OutputScalar,
class PointScalar,
69 class OutputFieldType,
class InputPointsType>
72 using ExecutionSpace =
typename DeviceType::execution_space;
73 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
74 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
75 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
77 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
78 using TeamMember =
typename TeamPolicy::member_type;
82 OutputFieldType output_;
83 InputPointsType inputPoints_;
85 ordinal_type polyOrder_;
86 ordinal_type numFields_, numPoints_;
88 size_t fad_size_output_;
90 static constexpr ordinal_type numVertices = 4;
91 static constexpr ordinal_type numEdges = 6;
92 static constexpr ordinal_type numEdgesPerFace = 3;
93 static constexpr ordinal_type numFaceFamilies = 2;
94 static constexpr ordinal_type numFaces = 4;
95 static constexpr ordinal_type numVerticesPerFace = 3;
96 static constexpr ordinal_type numInteriorFamilies = 3;
99 const ordinal_type face_vertices[numFaces*numVerticesPerFace] = {0,1,2,
108 const ordinal_type face_edges[numFaces * numEdgesPerFace] = {0,1,2,
114 const ordinal_type edge_start_[numEdges] = {0,1,0,0,1,2};
115 const ordinal_type edge_end_[numEdges] = {1,2,2,3,3,3};
116 const ordinal_type face_family_start_ [numFaceFamilies] = {0,1};
117 const ordinal_type face_family_middle_[numFaceFamilies] = {1,2};
118 const ordinal_type face_family_end_ [numFaceFamilies] = {2,0};
120 const ordinal_type numEdgeFunctions_;
121 const ordinal_type numFaceFunctionsPerFace_;
122 const ordinal_type numFaceFunctions_;
123 const ordinal_type numInteriorFunctionsPerFamily_;
124 const ordinal_type numInteriorFunctions_;
127 const ordinal_type faceOrdinalForInterior_[numInteriorFamilies] = {0,2,3};
128 const ordinal_type faceFamilyForInterior_[numInteriorFamilies] = {0,0,1};
129 const ordinal_type interiorCoordinateOrdinal_[numInteriorFamilies] = {3,0,1};
131 KOKKOS_INLINE_FUNCTION
132 ordinal_type dofOrdinalForFace(
const ordinal_type &faceOrdinal,
133 const ordinal_type &zeroBasedFaceFamily,
134 const ordinal_type &i,
135 const ordinal_type &j)
const
138 const ordinal_type faceDofOffset = numEdgeFunctions_ + faceOrdinal * numFaceFunctionsPerFace_;
144 const ordinal_type max_ij_sum = polyOrder_ - 1;
146 ordinal_type fieldOrdinal = faceDofOffset + zeroBasedFaceFamily;
148 for (ordinal_type ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
150 for (ordinal_type ii=0; ii<ij_sum; ii++)
153 const ordinal_type jj = ij_sum - ii;
154 if ( (ii == i) && (jj == j))
159 fieldOrdinal += numFaceFamilies;
166 : opType_(opType), output_(output), inputPoints_(inputPoints),
167 polyOrder_(polyOrder),
168 fad_size_output_(getScalarDimensionForView(output)),
169 numEdgeFunctions_(polyOrder * numEdges),
170 numFaceFunctionsPerFace_(polyOrder * (polyOrder-1)),
171 numFaceFunctions_(numFaceFunctionsPerFace_*numFaces),
172 numInteriorFunctionsPerFamily_((polyOrder > 2) ? (polyOrder-2)*(polyOrder-1)*polyOrder/6 : 0),
173 numInteriorFunctions_(numInteriorFunctionsPerFamily_ * numInteriorFamilies)
175 numFields_ = output.extent_int(0);
176 numPoints_ = output.extent_int(1);
178 const ordinal_type expectedCardinality = numEdgeFunctions_ + numFaceFunctions_ + numInteriorFunctions_;
184 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
185 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument,
"output field size does not match basis cardinality");
188 KOKKOS_INLINE_FUNCTION
189 void computeEdgeLegendre(OutputScratchView &P,
190 const ordinal_type &edgeOrdinal,
191 const PointScalar* lambda)
const
193 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
194 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
196 Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
199 KOKKOS_INLINE_FUNCTION
200 void edgeFunctionValue(OutputScalar &edgeValue_x,
201 OutputScalar &edgeValue_y,
202 OutputScalar &edgeValue_z,
203 const ordinal_type &edgeOrdinal,
204 OutputScratchView &P,
205 const ordinal_type &i,
206 const PointScalar* lambda,
207 const PointScalar* lambda_dx,
208 const PointScalar* lambda_dy,
209 const PointScalar* lambda_dz
212 const auto & s0 = lambda [edge_start_[edgeOrdinal]];
213 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
214 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
215 const auto & s0_dz = lambda_dz[edge_start_[edgeOrdinal]];
217 const auto & s1 = lambda [ edge_end_[edgeOrdinal]];
218 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
219 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
220 const auto & s1_dz = lambda_dz[ edge_end_[edgeOrdinal]];
222 const auto & P_i = P(i);
223 const PointScalar xWeight = s0 * s1_dx - s1 * s0_dx;
224 const PointScalar yWeight = s0 * s1_dy - s1 * s0_dy;
225 const PointScalar zWeight = s0 * s1_dz - s1 * s0_dz;
226 edgeValue_x = P_i * xWeight;
227 edgeValue_y = P_i * yWeight;
228 edgeValue_z = P_i * zWeight;
231 KOKKOS_INLINE_FUNCTION
232 void computeFaceIntegratedJacobi(OutputScratchView &L_2ip1,
233 const ordinal_type &zeroBasedFaceOrdinal,
234 const ordinal_type &zeroBasedFamilyOrdinal,
235 const ordinal_type &i,
236 const PointScalar* lambda)
const
238 const auto &s0_vertex_number = face_family_start_ [zeroBasedFamilyOrdinal];
239 const auto &s1_vertex_number = face_family_middle_[zeroBasedFamilyOrdinal];
240 const auto &s2_vertex_number = face_family_end_ [zeroBasedFamilyOrdinal];
243 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + s0_vertex_number];
244 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + s1_vertex_number];
245 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + s2_vertex_number];
247 const auto & s0 = lambda[s0_index];
248 const auto & s1 = lambda[s1_index];
249 const auto & s2 = lambda[s2_index];
250 const PointScalar jacobiScaling = s0 + s1 + s2;
252 const double alpha = i*2.0 + 1;
253 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);
256 KOKKOS_INLINE_FUNCTION
257 void faceFunctionValue(OutputScalar &value_x,
258 OutputScalar &value_y,
259 OutputScalar &value_z,
260 const ordinal_type &j,
261 const OutputScratchView &L_2ip1,
262 const OutputScalar &edgeValue_x,
263 const OutputScalar &edgeValue_y,
264 const OutputScalar &edgeValue_z,
265 const PointScalar* lambda)
const
267 const auto & L_2ip1_j = L_2ip1(j);
268 value_x = edgeValue_x * L_2ip1_j;
269 value_y = edgeValue_y * L_2ip1_j;
270 value_z = edgeValue_z * L_2ip1_j;
273 KOKKOS_INLINE_FUNCTION
274 void operator()(
const TeamMember & teamMember )
const
276 const ordinal_type numFunctionsPerFace = polyOrder_ * (polyOrder_ - 1);
277 auto pointOrdinal = teamMember.league_rank();
278 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
279 if (fad_size_output_ > 0) {
280 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
281 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
282 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
283 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
286 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
287 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
288 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
289 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
292 const auto & x = inputPoints_(pointOrdinal,0);
293 const auto & y = inputPoints_(pointOrdinal,1);
294 const auto & z = inputPoints_(pointOrdinal,2);
297 const PointScalar lambda[4] = {1. - x - y - z, x, y, z};
298 const PointScalar lambda_dx[4] = {-1., 1., 0., 0.};
299 const PointScalar lambda_dy[4] = {-1., 0., 1., 0.};
300 const PointScalar lambda_dz[4] = {-1., 0., 0., 1.};
302 const int num1DEdgeFunctions = polyOrder_;
311 auto & P = edge_field_values_at_point;
313 int fieldOrdinalOffset = 0;
314 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
316 computeEdgeLegendre(P, edgeOrdinal, lambda);
318 for (
int i=0; i<num1DEdgeFunctions; i++)
320 auto &output_x = output_(i+fieldOrdinalOffset,pointOrdinal,0);
321 auto &output_y = output_(i+fieldOrdinalOffset,pointOrdinal,1);
322 auto &output_z = output_(i+fieldOrdinalOffset,pointOrdinal,2);
324 edgeFunctionValue(output_x, output_y, output_z,
326 lambda, lambda_dx, lambda_dy, lambda_dz);
328 fieldOrdinalOffset += num1DEdgeFunctions;
334 auto & L_2ip1 = jacobi_values_at_point;
337 const int max_ij_sum = polyOrder_ - 1;
340 int faceFieldOrdinalOffset = fieldOrdinalOffset;
341 for (
int faceOrdinal = 0; faceOrdinal < numFaces; faceOrdinal++) {
342 for (
int familyOrdinal=1; familyOrdinal<=numFaceFamilies; familyOrdinal++)
344 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
346 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
348 for (
int i=0; i<ij_sum; i++)
350 computeFaceIntegratedJacobi(L_2ip1, faceOrdinal, familyOrdinal-1, i, lambda);
352 const int j = ij_sum - i;
354 const int faceEdgeOrdinal = familyOrdinal-1;
355 const int volumeEdgeOrdinal = face_edges[faceOrdinal * numEdgesPerFace + faceEdgeOrdinal];
356 const int edgeBasisOrdinal = i + volumeEdgeOrdinal*num1DEdgeFunctions;
357 const auto & edgeValue_x = output_(edgeBasisOrdinal,pointOrdinal,0);
358 const auto & edgeValue_y = output_(edgeBasisOrdinal,pointOrdinal,1);
359 const auto & edgeValue_z = output_(edgeBasisOrdinal,pointOrdinal,2);
361 auto & output_x = output_(fieldOrdinal,pointOrdinal,0);
362 auto & output_y = output_(fieldOrdinal,pointOrdinal,1);
363 auto & output_z = output_(fieldOrdinal,pointOrdinal,2);
365 faceFunctionValue(output_x, output_y, output_z, j, L_2ip1, edgeValue_x, edgeValue_y, edgeValue_z, lambda);
367 fieldOrdinal += numFaceFamilies;
370 fieldOrdinalOffset = fieldOrdinal - numFaceFamilies + 1;
372 faceFieldOrdinalOffset += numFunctionsPerFace;
378 const int interiorFieldOrdinalOffset = fieldOrdinalOffset;
379 const int min_ijk_sum = 2;
380 const int max_ijk_sum = polyOrder_-1;
383 const auto & L_2ipj = jacobi_values_at_point;
384 for (
int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
389 const auto & lambda_m = lambda[interiorCoordinateOrdinal_[interiorFamilyOrdinal-1]];
390 const PointScalar jacobiScaling = 1.0;
392 ordinal_type fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
394 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
395 const ordinal_type relatedFaceFamily = faceFamilyForInterior_ [interiorFamilyOrdinal-1];
397 for (
int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
399 for (
int i=0; i<ijk_sum-1; i++)
401 for (
int j=1; j<ijk_sum-i; j++)
404 const ordinal_type faceDofOrdinal = dofOrdinalForFace(relatedFaceOrdinal, relatedFaceFamily, i, j);
406 const double alpha = 2 * (i + j);
408 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipj, alpha, polyOrder_-1, lambda_m, jacobiScaling);
410 const int k = ijk_sum - i - j;
411 const auto & L_k = L_2ipj(k);
412 for (
int d=0; d<3; d++)
414 const auto & E_face_d = output_(faceDofOrdinal,pointOrdinal,d);
415 output_(fieldOrdinal,pointOrdinal,d) = L_k * E_face_d;
417 fieldOrdinal += numInteriorFamilies;
421 fieldOrdinalOffset = fieldOrdinal - numInteriorFamilies + 1;
430 int fieldOrdinalOffset = 0;
437 auto & P_i = edge_field_values_at_point;
438 auto & L_2ip1_j = jacobi_values_at_point;
439 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
441 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
442 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
444 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
445 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
446 const auto & s0_dz = lambda_dz[edge_start_[edgeOrdinal]];
447 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
448 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
449 const auto & s1_dz = lambda_dz[ edge_end_[edgeOrdinal]];
451 const OutputScalar grad_s0_cross_grad_s1[3] = {s0_dy * s1_dz - s0_dz * s1_dy,
452 s0_dz * s1_dx - s0_dx * s1_dz,
453 s0_dx * s1_dy - s0_dy * s1_dx};
455 Polynomials::shiftedScaledLegendreValues(P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
456 for (
int i=0; i<num1DEdgeFunctions; i++)
458 for (
int d=0; d<3; d++)
460 output_(i+fieldOrdinalOffset,pointOrdinal,d) = (i+2) * P_i(i) * grad_s0_cross_grad_s1[d];
463 fieldOrdinalOffset += num1DEdgeFunctions;
484 auto & P_2ip1_j = other_values_at_point;
485 auto & L_2ip1_j_dt = other_values2_at_point;
488 int faceFieldOrdinalOffset = fieldOrdinalOffset;
489 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
491 for (
int familyOrdinal=1; familyOrdinal<=numFaceFamilies; familyOrdinal++)
493 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
495 const auto &s0_vertex_number = face_family_start_ [familyOrdinal-1];
496 const auto &s1_vertex_number = face_family_middle_[familyOrdinal-1];
497 const auto &s2_vertex_number = face_family_end_ [familyOrdinal-1];
500 const auto &s0_index = face_vertices[faceOrdinal * numVerticesPerFace + s0_vertex_number];
501 const auto &s1_index = face_vertices[faceOrdinal * numVerticesPerFace + s1_vertex_number];
502 const auto &s2_index = face_vertices[faceOrdinal * numVerticesPerFace + s2_vertex_number];
504 const auto & s0 = lambda[s0_index];
505 const auto & s1 = lambda[s1_index];
506 const auto & s2 = lambda[s2_index];
507 const PointScalar jacobiScaling = s0 + s1 + s2;
509 const auto & s0_dx = lambda_dx[s0_index];
510 const auto & s0_dy = lambda_dy[s0_index];
511 const auto & s0_dz = lambda_dz[s0_index];
512 const auto & s1_dx = lambda_dx[s1_index];
513 const auto & s1_dy = lambda_dy[s1_index];
514 const auto & s1_dz = lambda_dz[s1_index];
515 const auto & s2_dx = lambda_dx[s2_index];
516 const auto & s2_dy = lambda_dy[s2_index];
517 const auto & s2_dz = lambda_dz[s2_index];
519 const PointScalar grad_s2[3] = {s2_dx, s2_dy, s2_dz};
520 const PointScalar gradJacobiScaling[3] = {s0_dx + s1_dx + s2_dx,
521 s0_dy + s1_dy + s2_dy,
522 s0_dz + s1_dz + s2_dz};
524 const PointScalar grad_s0_cross_grad_s1[3] = {s0_dy * s1_dz - s0_dz * s1_dy,
525 s0_dz * s1_dx - s0_dx * s1_dz,
526 s0_dx * s1_dy - s0_dy * s1_dx};
528 const PointScalar s0_grad_s1_minus_s1_grad_s0[3] = {s0 * s1_dx - s1 * s0_dx,
529 s0 * s1_dy - s1 * s0_dy,
530 s0 * s1_dz - s1 * s0_dz};
532 Polynomials::shiftedScaledLegendreValues (P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
539 const int max_ij_sum = polyOrder_ - 1;
540 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
542 for (
int i=0; i<ij_sum; i++)
544 const int j = ij_sum - i;
546 const double alpha = i*2.0 + 1;
548 Polynomials::shiftedScaledJacobiValues (P_2ip1_j, alpha, polyOrder_-1, PointScalar(s2), jacobiScaling);
549 Polynomials::shiftedScaledIntegratedJacobiValues (L_2ip1_j, alpha, polyOrder_-1, PointScalar(s2), jacobiScaling);
550 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2ip1_j_dt, alpha, polyOrder_-1, PointScalar(s2), jacobiScaling);
552 const PointScalar & edgeValue = P_i(i);
554 PointScalar grad_L_2ip1_j[3];
555 for (
int d=0; d<3; d++)
557 grad_L_2ip1_j[d] = P_2ip1_j(j-1) * grad_s2[d]
558 + L_2ip1_j_dt(j) * gradJacobiScaling[d];
561 const PointScalar grad_L_2ip1_j_cross_E_i[3] = { grad_L_2ip1_j[1] * edgeValue * s0_grad_s1_minus_s1_grad_s0[2] - grad_L_2ip1_j[2] * edgeValue * s0_grad_s1_minus_s1_grad_s0[1],
562 grad_L_2ip1_j[2] * edgeValue * s0_grad_s1_minus_s1_grad_s0[0] - grad_L_2ip1_j[0] * edgeValue * s0_grad_s1_minus_s1_grad_s0[2],
563 grad_L_2ip1_j[0] * edgeValue * s0_grad_s1_minus_s1_grad_s0[1] - grad_L_2ip1_j[1] * edgeValue * s0_grad_s1_minus_s1_grad_s0[0] };
565 for (
int d=0; d<3; d++)
567 const OutputScalar edgeCurl_d = (i+2.) * P_i(i) * grad_s0_cross_grad_s1[d];
568 output_(fieldOrdinal,pointOrdinal,d) = L_2ip1_j(j) * edgeCurl_d
569 + grad_L_2ip1_j_cross_E_i[d];
572 fieldOrdinal += numFaceFamilies;
575 fieldOrdinalOffset = fieldOrdinal - numFaceFamilies + 1;
577 faceFieldOrdinalOffset += numFunctionsPerFace;
583 auto & L_2ipj = jacobi_values_at_point;
584 auto & P_2ipj = other_values_at_point;
585 auto & L_2ip1 = edge_field_values_at_point;
586 auto & P = other_values2_at_point;
588 const int interiorFieldOrdinalOffset = fieldOrdinalOffset;
589 const int min_ijk_sum = 2;
590 const int max_ijk_sum = polyOrder_-1;
591 for (
int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
595 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
596 const ordinal_type relatedFaceFamily = faceFamilyForInterior_ [interiorFamilyOrdinal-1];
599 const auto & m = interiorCoordinateOrdinal_[interiorFamilyOrdinal-1];
600 const auto & lambda_m = lambda[m];
601 const PointScalar jacobiScaling = 1.0;
603 ordinal_type fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
605 for (
int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
607 for (
int i=0; i<ijk_sum-1; i++)
609 computeFaceIntegratedJacobi(L_2ip1, relatedFaceOrdinal, relatedFaceFamily, i, lambda);
611 const ordinal_type faceEdgeOrdinal = relatedFaceFamily;
612 const int volumeEdgeOrdinal = face_edges[relatedFaceOrdinal * numEdgesPerFace + faceEdgeOrdinal];
613 computeEdgeLegendre(P, volumeEdgeOrdinal, lambda);
615 OutputScalar edgeValue[3];
616 edgeFunctionValue(edgeValue[0], edgeValue[1], edgeValue[2], volumeEdgeOrdinal, P, i, lambda, lambda_dx, lambda_dy, lambda_dz);
618 for (
int j=1; j<ijk_sum-i; j++)
621 const ordinal_type faceDofOrdinal = dofOrdinalForFace(relatedFaceOrdinal, relatedFaceFamily, i, j);
623 const double alpha = 2 * (i + j);
625 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipj, alpha, polyOrder_-1, lambda_m, jacobiScaling);
626 Polynomials::shiftedScaledJacobiValues (P_2ipj, alpha, polyOrder_-1, lambda_m, jacobiScaling);
631 const int k = ijk_sum - i - j;
632 const auto & L_k = L_2ipj(k);
633 const auto & P_km1 = P_2ipj(k-1);
635 const PointScalar grad_L_k[3] = {P_km1 * lambda_dx[m],
636 P_km1 * lambda_dy[m],
637 P_km1 * lambda_dz[m]};
640 OutputScalar E_face[3];
641 faceFunctionValue(E_face[0], E_face[1], E_face[2], j, L_2ip1, edgeValue[0], edgeValue[1], edgeValue[2], lambda);
643 PointScalar grad_L_k_cross_E_face[3] = {grad_L_k[1] * E_face[2] - grad_L_k[2] * E_face[1],
644 grad_L_k[2] * E_face[0] - grad_L_k[0] * E_face[2],
645 grad_L_k[0] * E_face[1] - grad_L_k[1] * E_face[0]};
646 for (
int d=0; d<3; d++)
648 const auto & curl_E_face_d = output_(faceDofOrdinal,pointOrdinal,d);
649 output_(fieldOrdinal,pointOrdinal,d) = L_k * curl_E_face_d + grad_L_k_cross_E_face[d];
652 fieldOrdinal += numInteriorFamilies;
656 fieldOrdinalOffset = fieldOrdinal - numInteriorFamilies + 1;
672 INTREPID2_TEST_FOR_ABORT(
true,
673 ">>> ERROR: (Intrepid2::Hierarchical_HCURL_TET_Functor) Unsupported differential operator");
676 device_assert(
false);
683 size_t team_shmem_size (
int team_size)
const
686 size_t shmem_size = 0;
687 if (fad_size_output_ > 0)
688 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
690 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
710 template<
typename DeviceType,
711 typename OutputScalar = double,
712 typename PointScalar = double,
713 bool useCGBasis =
true>
715 :
public Basis<DeviceType,OutputScalar,PointScalar>
738 polyOrder_(polyOrder)
740 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
744 const int numEdgeFunctions = polyOrder * numEdges;
745 const int numFaceFunctions = polyOrder * (polyOrder-1) * numFaces;
746 const int numInteriorFunctionsPerFamily = (polyOrder > 2) ? (polyOrder-2)*(polyOrder-1)*polyOrder/6 : 0;
747 const int numInteriorFunctions = numInteriorFunctionsPerFamily * 3;
748 this->
basisCardinality_ = numEdgeFunctions + numFaceFunctions + numInteriorFunctions;
755 const int degreeLength = 1;
758 int fieldOrdinalOffset = 0;
763 const int numFunctionsPerEdge = polyOrder;
764 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
766 for (
int i=0; i<numFunctionsPerEdge; i++)
770 fieldOrdinalOffset += numFunctionsPerEdge;
772 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != numEdgeFunctions, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
775 const int max_ij_sum = polyOrder-1;
776 int faceFieldOrdinalOffset = fieldOrdinalOffset;
777 const int numFaceFamilies = 2;
778 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
780 for (
int faceFamilyOrdinal=1; faceFamilyOrdinal<=numFaceFamilies; faceFamilyOrdinal++)
783 int fieldOrdinal = faceFieldOrdinalOffset + faceFamilyOrdinal - 1;
784 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
786 for (
int i=0; i<ij_sum; i++)
789 fieldOrdinal += numFaceFamilies;
792 fieldOrdinalOffset = fieldOrdinal - numFaceFamilies + 1;
794 faceFieldOrdinalOffset += numFaceFunctions / numFaces;
796 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != numEdgeFunctions + numFaceFunctions, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
798 const int numInteriorFamilies = 3;
799 const int interiorFieldOrdinalOffset = fieldOrdinalOffset;
800 const int min_ijk_sum = 2;
801 const int max_ijk_sum = polyOrder-1;
802 for (
int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
805 int fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
806 for (
int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
808 for (
int i=0; i<ijk_sum-1; i++)
810 for (
int j=1; j<ijk_sum-i; j++)
813 fieldOrdinal += numInteriorFamilies;
817 fieldOrdinalOffset = fieldOrdinal - numInteriorFamilies + 1;
820 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
827 const int intrepid2FaceOrdinals[4] {3,0,1,2};
832 const ordinal_type tagSize = 4;
833 const ordinal_type posScDim = 0;
834 const ordinal_type posScOrd = 1;
835 const ordinal_type posDfOrd = 2;
838 const ordinal_type edgeDim = 1, faceDim = 2, volumeDim = 3;
843 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
845 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
847 tagView(tagNumber*tagSize+0) = edgeDim;
848 tagView(tagNumber*tagSize+1) = edgeOrdinal;
849 tagView(tagNumber*tagSize+2) = functionOrdinal;
850 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge;
854 const int numFunctionsPerFace = numFaceFunctions / numFaces;
855 for (
int faceOrdinalESEAS=0; faceOrdinalESEAS<numFaces; faceOrdinalESEAS++)
857 int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
858 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
860 tagView(tagNumber*tagSize+0) = faceDim;
861 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2;
862 tagView(tagNumber*tagSize+2) = functionOrdinal;
863 tagView(tagNumber*tagSize+3) = numFunctionsPerFace;
869 for (
int functionOrdinal=0; functionOrdinal<numInteriorFunctions; functionOrdinal++)
871 tagView(tagNumber*tagSize+0) = volumeDim;
872 tagView(tagNumber*tagSize+1) = 0;
873 tagView(tagNumber*tagSize+2) = functionOrdinal;
874 tagView(tagNumber*tagSize+3) = numInteriorFunctions;
882 for (ordinal_type i=0;i<cardinality;++i) {
883 tagView(i*tagSize+0) = volumeDim;
884 tagView(i*tagSize+1) = 0;
885 tagView(i*tagSize+2) = i;
886 tagView(i*tagSize+3) = cardinality;
908 return "Intrepid2_HierarchicalBasis_HCURL_TET";
941 const EOperator operatorType = OPERATOR_VALUE )
const override
943 auto numPoints = inputPoints.extent_int(0);
947 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
949 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
950 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
951 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
952 const int teamSize = 1;
954 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
955 Kokkos::parallel_for(
"Hierarchical_HCURL_TET_Functor", policy , functor);
966 BasisPtr<DeviceType,OutputScalar,PointScalar>
973 return Teuchos::rcp(
new HVOL_Line(this->
basisDegree_-1));
975 else if (subCellDim == 2)
979 INTREPID2_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Input parameters out of bounds");
986 virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>
988 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
990 return Teuchos::rcp(
new HostBasisType(polyOrder_) );
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
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.
const char * getName() const override
Returns basis name.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
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...
virtual bool requireOrientation() const override
True if orientation is required.
H(curl) basis on the triangle using a construction involving Legendre and integrated Jacobi polynomia...
For mathematical details of the construction, see:
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.
For mathematical details of the construction, see:
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
H(vol) basis on the line based on Legendre polynomials.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
ordinal_type getDegree() const
Returns the degree of the basis.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
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...
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
HierarchicalBasis_HCURL_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
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.
Functor for computing values for the HierarchicalBasis_HCURL_TET class.
Header file for the abstract base class Intrepid2::Basis.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.