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_;
86 int numFields_, numPoints_;
88 size_t fad_size_output_;
90 static const int numVertices = 3;
91 static const int numEdges = 3;
92 static const int numFaceFamilies = 2;
93 const int edge_start_[numEdges] = {0,1,0};
94 const int edge_end_[numEdges] = {1,2,2};
95 const int face_family_start_[numFaceFamilies] = {0,1};
96 const int face_family_middle_[numFaceFamilies] = {1,2};
97 const int face_family_end_ [numFaceFamilies] = {2,0};
100 : opType_(opType), output_(output), inputPoints_(inputPoints),
101 polyOrder_(polyOrder),
104 numFields_ = output.extent_int(0);
105 numPoints_ = output.extent_int(1);
106 const int expectedCardinality = 3 * polyOrder_ + polyOrder_ * (polyOrder-1);
108 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
109 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument,
"output field size does not match basis cardinality");
112 KOKKOS_INLINE_FUNCTION
113 void operator()(
const TeamMember & teamMember )
const
115 auto pointOrdinal = teamMember.league_rank();
116 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
117 if (fad_size_output_ > 0) {
118 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
119 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
120 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
121 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
124 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
125 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
126 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
127 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
130 const auto & x = inputPoints_(pointOrdinal,0);
131 const auto & y = inputPoints_(pointOrdinal,1);
134 const PointScalar lambda[3] = {1. - x - y, x, y};
135 const PointScalar lambda_dx[3] = {-1., 1., 0.};
136 const PointScalar lambda_dy[3] = {-1., 0., 1.};
138 const int num1DEdgeFunctions = polyOrder_;
145 int fieldOrdinalOffset = 0;
146 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
148 const auto & s0 = lambda [edge_start_[edgeOrdinal]];
149 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
150 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
152 const auto & s1 = lambda [ edge_end_[edgeOrdinal]];
153 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
154 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
156 Polynomials::shiftedScaledLegendreValues(edge_field_values_at_point, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
157 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
159 const auto & legendreValue = edge_field_values_at_point(edgeFunctionOrdinal);
160 const PointScalar xWeight = s0 * s1_dx - s1 * s0_dx;
161 const PointScalar yWeight = s0 * s1_dy - s1 * s0_dy;
162 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = legendreValue * xWeight;
163 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = legendreValue * yWeight;
165 fieldOrdinalOffset += num1DEdgeFunctions;
171 const double jacobiScaling = 1.0;
173 const int max_ij_sum = polyOrder_ - 1;
176 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
177 for (
int familyOrdinal=1; familyOrdinal<=2; familyOrdinal++)
179 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
180 const auto &s2 = lambda[ face_family_end_[familyOrdinal-1]];
181 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
183 for (
int i=0; i<ij_sum; i++)
185 const int j = ij_sum - i;
187 const int edgeBasisOrdinal = i + (familyOrdinal-1)*num1DEdgeFunctions;
188 const auto & edgeValue_x = output_(edgeBasisOrdinal,pointOrdinal,0);
189 const auto & edgeValue_y = output_(edgeBasisOrdinal,pointOrdinal,1);
190 const double alpha = i*2.0 + 1;
192 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-1, s2, jacobiScaling);
193 const auto & jacobiValue = jacobi_values_at_point(j);
194 output_(fieldOrdinal,pointOrdinal,0) = edgeValue_x * jacobiValue;
195 output_(fieldOrdinal,pointOrdinal,1) = edgeValue_y * jacobiValue;
207 int fieldOrdinalOffset = 0;
214 auto & P_i = edge_field_values_at_point;
215 auto & L_2ip1_j = jacobi_values_at_point;
216 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
218 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
219 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
221 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
222 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
223 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
224 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
226 const OutputScalar grad_s0_cross_grad_s1 = s0_dx * s1_dy - s1_dx * s0_dy;
228 Polynomials::shiftedScaledLegendreValues(P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
229 for (
int i=0; i<num1DEdgeFunctions; i++)
231 output_(i+fieldOrdinalOffset,pointOrdinal) = (i+2) * P_i(i) * grad_s0_cross_grad_s1;
233 fieldOrdinalOffset += num1DEdgeFunctions;
254 auto & P_2ip1_j = other_values_at_point;
257 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
258 for (
int familyOrdinal=1; familyOrdinal<=2; familyOrdinal++)
260 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
262 const auto &s0_index = face_family_start_ [familyOrdinal-1];
263 const auto &s1_index = face_family_middle_[familyOrdinal-1];
264 const auto &s2_index = face_family_end_ [familyOrdinal-1];
265 const auto &s0 = lambda[s0_index];
266 const auto &s1 = lambda[s1_index];
267 const auto &s2 = lambda[s2_index];
268 const double jacobiScaling = 1.0;
270 const auto & s0_dx = lambda_dx[s0_index];
271 const auto & s0_dy = lambda_dy[s0_index];
272 const auto & s1_dx = lambda_dx[s1_index];
273 const auto & s1_dy = lambda_dy[s1_index];
274 const auto & s2_dx = lambda_dx[s2_index];
275 const auto & s2_dy = lambda_dy[s2_index];
277 const OutputScalar grad_s0_cross_grad_s1 = s0_dx * s1_dy - s1_dx * s0_dy;
279 Polynomials::shiftedScaledLegendreValues (P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
283 const PointScalar xEdgeWeight = s0 * s1_dx - s1 * s0_dx;
284 const PointScalar yEdgeWeight = s0 * s1_dy - s1 * s0_dy;
285 OutputScalar grad_s2_cross_xy_edgeWeight = s2_dx * yEdgeWeight - xEdgeWeight * s2_dy;
287 const int max_ij_sum = polyOrder_ - 1;
288 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
290 for (
int i=0; i<ij_sum; i++)
292 const int j = ij_sum - i;
293 const OutputScalar edgeCurl = (i+2.) * P_i(i) * grad_s0_cross_grad_s1;
295 const double alpha = i*2.0 + 1;
297 Polynomials::shiftedScaledJacobiValues(P_2ip1_j, alpha, polyOrder_-1, PointScalar(s2), jacobiScaling);
298 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ip1_j, alpha, polyOrder_-1, s2, jacobiScaling);
300 const PointScalar & edgeValue = P_i(i);
301 output_(fieldOrdinal,pointOrdinal) = L_2ip1_j(j) * edgeCurl + P_2ip1_j(j-1) * edgeValue * grad_s2_cross_xy_edgeWeight;
320 INTREPID2_TEST_FOR_ABORT(
true,
321 ">>> ERROR: (Intrepid2::Hierarchical_HCURL_TRI_Functor) Unsupported differential operator");
331 size_t team_shmem_size (
int team_size)
const
334 size_t shmem_size = 0;
335 if (fad_size_output_ > 0)
336 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
338 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);