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 bool defineVertexFunctions_;
87 int numFields_, numPoints_;
89 size_t fad_size_output_;
91 static const int numVertices = 3;
92 static const int numEdges = 3;
93 const int edge_start_[numEdges] = {0,1,0};
94 const int edge_end_[numEdges] = {1,2,2};
97 int polyOrder,
bool defineVertexFunctions)
98 : opType_(opType), output_(output), inputPoints_(inputPoints),
99 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
102 numFields_ = output.extent_int(0);
103 numPoints_ = output.extent_int(1);
104 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
105 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument,
"output field size does not match basis cardinality");
108 KOKKOS_INLINE_FUNCTION
109 void operator()(
const TeamMember & teamMember )
const
111 auto pointOrdinal = teamMember.league_rank();
112 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
113 if (fad_size_output_ > 0) {
114 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
115 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
116 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
117 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
120 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
121 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
122 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
123 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
126 const auto & x = inputPoints_(pointOrdinal,0);
127 const auto & y = inputPoints_(pointOrdinal,1);
130 const PointScalar lambda[3] = {1. - x - y, x, y};
131 const PointScalar lambda_dx[3] = {-1., 1., 0.};
132 const PointScalar lambda_dy[3] = {-1., 0., 1.};
134 const int num1DEdgeFunctions = polyOrder_ - 1;
141 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
143 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
145 if (!defineVertexFunctions_)
149 output_(0,pointOrdinal) = 1.0;
153 int fieldOrdinalOffset = 3;
154 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
156 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
157 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
159 Polynomials::shiftedScaledIntegratedLegendreValues(edge_field_values_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
160 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
163 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = edge_field_values_at_point(edgeFunctionOrdinal+2);
165 fieldOrdinalOffset += num1DEdgeFunctions;
171 const double jacobiScaling = 1.0;
173 const int max_ij_sum = polyOrder_;
176 const int min_ij_sum = min_i + min_j;
177 for (
int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
179 for (
int i=min_i; i<=ij_sum-min_j; i++)
181 const int j = ij_sum - i;
182 const int edgeBasisOrdinal = i+numVertices-2;
183 const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
184 const double alpha = i*2.0;
186 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-2, lambda[2], jacobiScaling);
187 const auto & jacobiValue = jacobi_values_at_point(j);
188 output_(fieldOrdinalOffset,pointOrdinal) = edgeValue * jacobiValue;
189 fieldOrdinalOffset++;
199 if (defineVertexFunctions_)
203 output_(0,pointOrdinal,0) = -1.0;
204 output_(0,pointOrdinal,1) = -1.0;
210 output_(0,pointOrdinal,0) = 0.0;
211 output_(0,pointOrdinal,1) = 0.0;
214 output_(1,pointOrdinal,0) = 1.0;
215 output_(1,pointOrdinal,1) = 0.0;
217 output_(2,pointOrdinal,0) = 0.0;
218 output_(2,pointOrdinal,1) = 1.0;
221 int fieldOrdinalOffset = 3;
233 auto & P_i_minus_1 = edge_field_values_at_point;
234 auto & L_i_dt = jacobi_values_at_point;
235 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
237 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
238 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
240 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
241 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
242 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
243 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
245 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
246 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
247 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
250 const int i = edgeFunctionOrdinal+2;
251 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
252 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
254 fieldOrdinalOffset += num1DEdgeFunctions;
275 auto & P_2i_j_minus_1 = edge_field_values_at_point;
276 auto & L_2i_j_dt = jacobi_values_at_point;
277 auto & L_i = other_values_at_point;
278 auto & L_2i_j = other_values2_at_point;
281 const double jacobiScaling = 1.0;
283 const int max_ij_sum = polyOrder_;
286 const int min_ij_sum = min_i + min_j;
287 for (
int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
289 for (
int i=min_i; i<=ij_sum-min_j; i++)
291 const int j = ij_sum - i;
293 const int edgeBasisOrdinal = i+numVertices-2;
294 const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
295 const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
297 const double alpha = i*2.0;
299 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, lambda[1], lambda[0]+lambda[1]);
300 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2i_j_dt, alpha, polyOrder_, lambda[2], jacobiScaling);
301 Polynomials::shiftedScaledIntegratedJacobiValues ( L_2i_j, alpha, polyOrder_, lambda[2], jacobiScaling);
302 Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1, alpha, polyOrder_-1, lambda[2], jacobiScaling);
304 const auto & s0_dx = lambda_dx[0];
305 const auto & s0_dy = lambda_dy[0];
306 const auto & s1_dx = lambda_dx[1];
307 const auto & s1_dy = lambda_dy[1];
308 const auto & s2_dx = lambda_dx[2];
309 const auto & s2_dy = lambda_dy[2];
311 const OutputScalar basisValue_dx = L_2i_j(j) * grad_L_i_dx + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dx + L_2i_j_dt(j) * (s0_dx + s1_dx + s2_dx));
312 const OutputScalar basisValue_dy = L_2i_j(j) * grad_L_i_dy + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dy + L_2i_j_dt(j) * (s0_dy + s1_dy + s2_dy));
314 output_(fieldOrdinalOffset,pointOrdinal,0) = basisValue_dx;
315 output_(fieldOrdinalOffset,pointOrdinal,1) = basisValue_dy;
316 fieldOrdinalOffset++;
331 INTREPID2_TEST_FOR_ABORT(
true,
332 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
342 size_t team_shmem_size (
int team_size)
const
345 size_t shmem_size = 0;
346 if (fad_size_output_ > 0)
347 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
349 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);