Intrepid2
Intrepid2_HierarchicalBasis_HCURL_TET.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Mauro Perego (mperego@sandia.gov) or
38// Nate Roberts (nvrober@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
48#ifndef Intrepid2_HierarchicalBasis_HCURL_TET_h
49#define Intrepid2_HierarchicalBasis_HCURL_TET_h
50
51#include <Intrepid2_config.h>
52
53#include <Kokkos_DynRankView.hpp>
54
55#include "Intrepid2_Basis.hpp"
59#include "Intrepid2_Utils.hpp"
60
61namespace Intrepid2
62{
68 template<class DeviceType, class OutputScalar, class PointScalar,
69 class OutputFieldType, class InputPointsType>
71 {
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>>;
76
77 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
78 using TeamMember = typename TeamPolicy::member_type;
79
80 EOperator opType_;
81
82 OutputFieldType output_; // F,P
83 InputPointsType inputPoints_; // P,D
84
85 ordinal_type polyOrder_;
86 ordinal_type numFields_, numPoints_;
87
88 size_t fad_size_output_;
89
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;
97
98 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
99 const ordinal_type face_vertices[numFaces*numVerticesPerFace] = {0,1,2, // face 0
100 0,1,3, // face 1
101 1,2,3, // face 2
102 0,2,3 // face 3
103 };
104
105 // index into face_edges with faceOrdinal * numEdgesPerFace + faceEdgeNumber
106 // (entries are the edge numbers in the tetrahedron)
107 // note that the orientation of each face's third edge is reversed relative to the orientation in the volume
108 const ordinal_type face_edges[numFaces * numEdgesPerFace] = {0,1,2, // face 0
109 0,4,3, // face 1
110 1,5,4, // face 2
111 2,5,3}; // face 3
112
113 // the following ordering of the edges matches that used by ESEAS
114 const ordinal_type edge_start_[numEdges] = {0,1,0,0,1,2}; // edge i is from edge_start_[i] to edge_end_[i]
115 const ordinal_type edge_end_[numEdges] = {1,2,2,3,3,3}; // edge i is from edge_start_[i] to edge_end_[i]
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};
119
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_;
125
126 // interior basis functions are computed in terms of certain face basis functions.
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}; // m, where E^b_{ijk} is computed in terms of [L^{2(i+j)}_k](1-lambda_m, lambda_m)
130
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
136 {
137 // determine where the functions for this face start
138 const ordinal_type faceDofOffset = numEdgeFunctions_ + faceOrdinal * numFaceFunctionsPerFace_;
139
140 // rather than deriving a closed formula in terms of i and j (which is potentially error-prone),
141 // we simply step through a for loop much as we do in the basis computations themselves. (This method
142 // is not expected to be called so much as to be worth optimizing.)
143
144 const ordinal_type max_ij_sum = polyOrder_ - 1;
145
146 ordinal_type fieldOrdinal = faceDofOffset + zeroBasedFaceFamily; // families are interleaved on the face.
147
148 for (ordinal_type ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
149 {
150 for (ordinal_type ii=0; ii<ij_sum; ii++)
151 {
152 // j will be ij_sum - i; j >= 1.
153 const ordinal_type jj = ij_sum - ii; // jj >= 1
154 if ( (ii == i) && (jj == j))
155 {
156 // have reached the (i,j) we're looking for
157 return fieldOrdinal;
158 }
159 fieldOrdinal += numFaceFamilies; // increment for the interleaving of face families.
160 }
161 }
162 return -1; // error: not found.
163 }
164
165 Hierarchical_HCURL_TET_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, int polyOrder)
166 : opType_(opType), output_(output), inputPoints_(inputPoints),
167 polyOrder_(polyOrder),
168 fad_size_output_(getScalarDimensionForView(output)),
169 numEdgeFunctions_(polyOrder * numEdges), // 6 edges
170 numFaceFunctionsPerFace_(polyOrder * (polyOrder-1)), // 2 families, each with p*(p-1)/2 functions per face
171 numFaceFunctions_(numFaceFunctionsPerFace_*numFaces), // 4 faces
172 numInteriorFunctionsPerFamily_((polyOrder > 2) ? (polyOrder-2)*(polyOrder-1)*polyOrder/6 : 0), // p choose 3
173 numInteriorFunctions_(numInteriorFunctionsPerFamily_ * numInteriorFamilies) // 3 families of interior functions
174 {
175 numFields_ = output.extent_int(0);
176 numPoints_ = output.extent_int(1);
177
178 const ordinal_type expectedCardinality = numEdgeFunctions_ + numFaceFunctions_ + numInteriorFunctions_;
179
180 // interior family I: computed in terms of face 012 (face ordinal 0), ordinal 0 in face family I. First interior family is computed in terms of the first set of face functions (note that both sets of families are interleaved, so basis ordinal increments are by numInteriorFamilies and numFaceFamilies, respectively).
181 // interior family II: computed in terms of face 123 (face ordinal 2), ordinal 2 in face family I.
182 // interior family III: computed in terms of face 230 (face ordinal 3), ordinal 3 in face family II.
183
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");
186 }
187
188 KOKKOS_INLINE_FUNCTION
189 void computeEdgeLegendre(OutputScratchView &P,
190 const ordinal_type &edgeOrdinal,
191 const PointScalar* lambda) const
192 {
193 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
194 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
195
196 Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
197 }
198
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
210 ) const
211 {
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]];
216
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]];
221
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;
229 }
230
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
237 {
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];
241
242 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
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];
246
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;
251
252 const double alpha = i*2.0 + 1;
253 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);
254 }
255
256 KOKKOS_INLINE_FUNCTION
257 void faceFunctionValue(OutputScalar &value_x,
258 OutputScalar &value_y,
259 OutputScalar &value_z,
260 const ordinal_type &j, // j >= 1
261 const OutputScratchView &L_2ip1, // container in which shiftedScaledIntegratedJacobiValues have been computed for (2i+1) for appropriate face and family
262 const OutputScalar &edgeValue_x,
263 const OutputScalar &edgeValue_y,
264 const OutputScalar &edgeValue_z,
265 const PointScalar* lambda) const
266 {
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;
271 }
272
273 KOKKOS_INLINE_FUNCTION
274 void operator()( const TeamMember & teamMember ) const
275 {
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_);
284 }
285 else {
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_);
290 }
291
292 const auto & x = inputPoints_(pointOrdinal,0);
293 const auto & y = inputPoints_(pointOrdinal,1);
294 const auto & z = inputPoints_(pointOrdinal,2);
295
296 // write as barycentric coordinates:
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.};
301
302 const int num1DEdgeFunctions = polyOrder_; // per edge
303
304 switch (opType_)
305 {
306 case OPERATOR_VALUE:
307 {
308 // edge functions
309
310 // relabel scratch view
311 auto & P = edge_field_values_at_point;
312
313 int fieldOrdinalOffset = 0;
314 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
315 {
316 computeEdgeLegendre(P, edgeOrdinal, lambda);
317
318 for (int i=0; i<num1DEdgeFunctions; i++)
319 {
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);
323
324 edgeFunctionValue(output_x, output_y, output_z,
325 edgeOrdinal, P, i,
326 lambda, lambda_dx, lambda_dy, lambda_dz);
327 }
328 fieldOrdinalOffset += num1DEdgeFunctions;
329 }
330
331 // face functions
332 {
333 // relabel scratch view
334 auto & L_2ip1 = jacobi_values_at_point;
335
336 // these functions multiply the edge functions from the 01 edge by integrated Jacobi functions, appropriately scaled
337 const int max_ij_sum = polyOrder_ - 1;
338
339 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
340 int faceFieldOrdinalOffset = fieldOrdinalOffset;
341 for (int faceOrdinal = 0; faceOrdinal < numFaces; faceOrdinal++) {
342 for (int familyOrdinal=1; familyOrdinal<=numFaceFamilies; familyOrdinal++)
343 {
344 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
345
346 for (int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
347 {
348 for (int i=0; i<ij_sum; i++)
349 {
350 computeFaceIntegratedJacobi(L_2ip1, faceOrdinal, familyOrdinal-1, i, lambda);
351
352 const int j = ij_sum - i; // j >= 1
353 // family 1 involves edge functions from edge (s0,s1) (edgeOrdinal 0 in the face); family 2 involves functions from edge (s1,s2) (edgeOrdinal 1 in the face)
354 const int faceEdgeOrdinal = familyOrdinal-1; // family I: use first edge from face; family II: use second edge
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);
360
361 auto & output_x = output_(fieldOrdinal,pointOrdinal,0);
362 auto & output_y = output_(fieldOrdinal,pointOrdinal,1);
363 auto & output_z = output_(fieldOrdinal,pointOrdinal,2);
364
365 faceFunctionValue(output_x, output_y, output_z, j, L_2ip1, edgeValue_x, edgeValue_y, edgeValue_z, lambda);
366
367 fieldOrdinal += numFaceFamilies; // increment due to the interleaving
368 } // i
369 } // ij_sum
370 fieldOrdinalOffset = fieldOrdinal - numFaceFamilies + 1; // due to the interleaving increment, we've gone numFaceFamilies past the last face ordinal. Set offset to be one past.
371 } // familyOrdinal
372 faceFieldOrdinalOffset += numFunctionsPerFace;
373 } // faceOrdinal
374 } // face functions block
375
376 // interior functions
377 {
378 const int interiorFieldOrdinalOffset = fieldOrdinalOffset;
379 const int min_ijk_sum = 2;
380 const int max_ijk_sum = polyOrder_-1;
381
382 // relabel Jacobi values container:
383 const auto & L_2ipj = jacobi_values_at_point;
384 for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
385 {
386 // following ESEAS, we interleave the interior families. This groups all the interior dofs of a given degree together.
387
388 // lambda_m is used to compute the appropriate weight in terms of Jacobi functions of order k below.
389 const auto & lambda_m = lambda[interiorCoordinateOrdinal_[interiorFamilyOrdinal-1]];
390 const PointScalar jacobiScaling = 1.0;
391
392 ordinal_type fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
393
394 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
395 const ordinal_type relatedFaceFamily = faceFamilyForInterior_ [interiorFamilyOrdinal-1]; // zero-based
396
397 for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
398 {
399 for (int i=0; i<ijk_sum-1; i++)
400 {
401 for (int j=1; j<ijk_sum-i; j++)
402 {
403 // the interior functions are blended face functions. This dof ordinal corresponds to the face function which we blend.
404 const ordinal_type faceDofOrdinal = dofOrdinalForFace(relatedFaceOrdinal, relatedFaceFamily, i, j);
405
406 const double alpha = 2 * (i + j);
407
408 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipj, alpha, polyOrder_-1, lambda_m, jacobiScaling);
409
410 const int k = ijk_sum - i - j;
411 const auto & L_k = L_2ipj(k);
412 for (int d=0; d<3; d++)
413 {
414 const auto & E_face_d = output_(faceDofOrdinal,pointOrdinal,d);
415 output_(fieldOrdinal,pointOrdinal,d) = L_k * E_face_d;
416 }
417 fieldOrdinal += numInteriorFamilies; // increment due to the interleaving.
418 }
419 }
420 }
421 fieldOrdinalOffset = fieldOrdinal - numInteriorFamilies + 1; // due to the interleaving increment, we've gone numInteriorFamilies past the last interior ordinal. Set offset to be one past.
422 }
423 } // interior functions block
424
425 } // end OPERATOR_VALUE
426 break;
427 case OPERATOR_CURL:
428 {
429 // edge functions
430 int fieldOrdinalOffset = 0;
431 /*
432 Per Fuentes et al. (see Appendix E.1, E.2), the curls of the edge functions, are
433 (i+2) * [P_i](s0,s1) * (grad s0 \times grad s1)
434 The P_i we have implemented in shiftedScaledLegendreValues.
435 */
436 // rename the scratch memory to match our usage here:
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++)
440 {
441 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
442 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
443
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]];
450
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};
454
455 Polynomials::shiftedScaledLegendreValues(P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
456 for (int i=0; i<num1DEdgeFunctions; i++)
457 {
458 for (int d=0; d<3; d++)
459 {
460 output_(i+fieldOrdinalOffset,pointOrdinal,d) = (i+2) * P_i(i) * grad_s0_cross_grad_s1[d];
461 }
462 }
463 fieldOrdinalOffset += num1DEdgeFunctions;
464 }
465
466 /*
467 Fuentes et al give the face functions as E^f_{ij}, with curl:
468 [L^{2i+1}_j](s0+s1,s2) curl(E^E_i(s0,s1)) + grad[L^(2i+1)_j](s0+s1,s2) \times E^E_i(s0,s1)
469 where:
470 - E^E_i is the ith edge function on the edge s0 to s1
471 - L^{2i+1}_j is an shifted, scaled integrated Jacobi polynomial.
472 - For family 1, s0s1s2 = 012
473 - For family 2, s0s1s2 = 120
474 - Note that grad[L^(2i+1)_j](s0+s1,s2) is computed as [P^{2i+1}_{j-1}](s0+s1,s2) (grad s2) + [R^{2i+1}_{j-1}] grad (s0+s1+s2),
475 but for triangles (s0+s1+s2) is always 1, so that the grad (s0+s1+s2) is 0.
476 - Here,
477 [P^{2i+1}_{j-1}](s0,s1) = P^{2i+1}_{j-1}(s1,s0+s1)
478 and
479 [R^{2i+1}_{j-1}(s0,s1)] = d/dt L^{2i+1}_j(s1,s0+s1)
480 We have implemented P^{alpha}_{j} as shiftedScaledJacobiValues,
481 and d/dt L^{alpha}_{j} as shiftedScaledIntegratedJacobiValues_dt.
482 */
483 // rename the scratch memory to match our usage here:
484 auto & P_2ip1_j = other_values_at_point;
485 auto & L_2ip1_j_dt = other_values2_at_point;
486
487 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
488 int faceFieldOrdinalOffset = fieldOrdinalOffset;
489 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
490 {
491 for (int familyOrdinal=1; familyOrdinal<=numFaceFamilies; familyOrdinal++)
492 {
493 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
494
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];
498
499 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
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];
503
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;
508
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];
518
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};
523
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};
527
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};
531
532 Polynomials::shiftedScaledLegendreValues (P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
533 // [L^{2i+1}_j](s0+s1,s2) curl(E^E_i(s0,s1)) + grad[L^(2i+1)_j](s0+s1,s2) \times E^E_i(s0,s1)
534 // - Note that grad[L^(2i+1)_j](s0+s1,s2) is computed as [P^{2i+1}_{j-1}](s0+s1,s2) (grad s2) + [R^{2i+1}_{j-1}](s0+s1,s2) grad (s0+s1+s2),
535 // - R^{2i+1}_{j-1}(s0+s1;s0+s1+s2) = d/dt L^{2i+1}_j(s0+s1;s0+s1+s2)
536 // - We have implemented d/dt L^{alpha}_{j} as shiftedScaledIntegratedJacobiValues_dt.
537 // - E^E_i(s0,s1) = [P_i](s0,s1) (s0 grad s1 - s1 grad s0)
538
539 const int max_ij_sum = polyOrder_ - 1;
540 for (int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
541 {
542 for (int i=0; i<ij_sum; i++)
543 {
544 const int j = ij_sum - i; // j >= 1
545
546 const double alpha = i*2.0 + 1;
547
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);
551
552 const PointScalar & edgeValue = P_i(i);
553
554 PointScalar grad_L_2ip1_j[3];
555 for (int d=0; d<3; d++)
556 {
557 grad_L_2ip1_j[d] = P_2ip1_j(j-1) * grad_s2[d] // [P^{2i+1}_{j-1}](s0+s1,s2) (grad s2)
558 + L_2ip1_j_dt(j) * gradJacobiScaling[d]; // [R^{2i+1}_{j-1}](s0+s1,s2) grad (s0+s1+s2)
559 }
560
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] };
564
565 for (int d=0; d<3; d++)
566 {
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 // [L^{2i+1}_j](s0+s1,s2) curl(E^E_i(s0,s1))
569 + grad_L_2ip1_j_cross_E_i[d];
570 }
571
572 fieldOrdinal += numFaceFamilies; // increment due to the interleaving
573 } // i
574 } // ij_sum
575 fieldOrdinalOffset = fieldOrdinal - numFaceFamilies + 1; // due to the interleaving increment, we've gone numFaceFamilies past the last face ordinal. Set offset to be one past.
576 } // familyOrdinal
577 faceFieldOrdinalOffset += numFunctionsPerFace;
578 } // faceOrdinal
579
580 // interior functions
581 {
582 // relabel values containers:
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;
587
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++)
592 {
593 // following ESEAS, we interleave the interior families. This groups all the interior dofs of a given degree together.
594
595 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
596 const ordinal_type relatedFaceFamily = faceFamilyForInterior_ [interiorFamilyOrdinal-1]; // zero-based
597
598 // lambda_m is used to compute the appropriate weight in terms of Jacobi functions of order k below.
599 const auto & m = interiorCoordinateOrdinal_[interiorFamilyOrdinal-1];
600 const auto & lambda_m = lambda[m];
601 const PointScalar jacobiScaling = 1.0;
602
603 ordinal_type fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
604
605 for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
606 {
607 for (int i=0; i<ijk_sum-1; i++)
608 {
609 computeFaceIntegratedJacobi(L_2ip1, relatedFaceOrdinal, relatedFaceFamily, i, lambda);
610 // face family 1 involves edge functions from edge (0,1) (edgeOrdinal 0); family 2 involves functions from edge (1,2) (edgeOrdinal 1)
611 const ordinal_type faceEdgeOrdinal = relatedFaceFamily;
612 const int volumeEdgeOrdinal = face_edges[relatedFaceOrdinal * numEdgesPerFace + faceEdgeOrdinal];
613 computeEdgeLegendre(P, volumeEdgeOrdinal, lambda);
614
615 OutputScalar edgeValue[3];
616 edgeFunctionValue(edgeValue[0], edgeValue[1], edgeValue[2], volumeEdgeOrdinal, P, i, lambda, lambda_dx, lambda_dy, lambda_dz);
617
618 for (int j=1; j<ijk_sum-i; j++)
619 {
620 // the interior functions are blended face functions. This dof ordinal corresponds to the face function which we blend.
621 const ordinal_type faceDofOrdinal = dofOrdinalForFace(relatedFaceOrdinal, relatedFaceFamily, i, j);
622
623 const double alpha = 2 * (i + j);
624
625 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipj, alpha, polyOrder_-1, lambda_m, jacobiScaling);
626 Polynomials::shiftedScaledJacobiValues (P_2ipj, alpha, polyOrder_-1, lambda_m, jacobiScaling);
627
628 // gradient of [L^{2(i+j)}_k](t0,t1) = [P^{2(i+j)}_{k-1}](t0,t1) grad t1 + [R^{2(i+j}_k](t0,t1) grad (t0+t1).
629 // we have t0 = lambda_m, t1 = 1 - lambda_m, so grad (t0 + t1) = 0.
630
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);
634
635 const PointScalar grad_L_k[3] = {P_km1 * lambda_dx[m],
636 P_km1 * lambda_dy[m],
637 P_km1 * lambda_dz[m]};
638
639 // compute E_face -- OPERATOR_VALUE for the face function corresponding to this interior function
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);
642
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++)
647 {
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];
650 }
651
652 fieldOrdinal += numInteriorFamilies; // increment due to the interleaving.
653 }
654 }
655 }
656 fieldOrdinalOffset = fieldOrdinal - numInteriorFamilies + 1; // due to the interleaving increment, we've gone numInteriorFamilies past the last interior ordinal. Set offset to be one past.
657 }
658 } // interior functions block
659 } // OPERATOR_CURL
660 break;
661 case OPERATOR_GRAD:
662 case OPERATOR_D1:
663 case OPERATOR_D2:
664 case OPERATOR_D3:
665 case OPERATOR_D4:
666 case OPERATOR_D5:
667 case OPERATOR_D6:
668 case OPERATOR_D7:
669 case OPERATOR_D8:
670 case OPERATOR_D9:
671 case OPERATOR_D10:
672 INTREPID2_TEST_FOR_ABORT( true,
673 ">>> ERROR: (Intrepid2::Hierarchical_HCURL_TET_Functor) Unsupported differential operator");
674 default:
675 // unsupported operator type
676 device_assert(false);
677 }
678 }
679
680 // Provide the shared memory capacity.
681 // This function takes the team_size as an argument,
682 // which allows team_size-dependent allocations.
683 size_t team_shmem_size (int team_size) const
684 {
685 // we will use shared memory to create a fast buffer for basis computations
686 size_t shmem_size = 0;
687 if (fad_size_output_ > 0)
688 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
689 else
690 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
691
692 return shmem_size;
693 }
694 };
695
710 template<typename DeviceType,
711 typename OutputScalar = double,
712 typename PointScalar = double,
713 bool useCGBasis = true> // if useCGBasis is false, all basis functions will be associated with the interior
715 : public Basis<DeviceType,OutputScalar,PointScalar>
716 {
717 public:
719
722
723 using typename BasisBase::OutputViewType;
724 using typename BasisBase::PointViewType;
725 using typename BasisBase::ScalarViewType;
726
727 using typename BasisBase::ExecutionSpace;
728
729 protected:
730 int polyOrder_; // the maximum order of the polynomial
731 public:
736 HierarchicalBasis_HCURL_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
737 :
738 polyOrder_(polyOrder)
739 {
740 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
741 const int numEdges = this->basisCellTopology_.getEdgeCount();
742 const int numFaces = this->basisCellTopology_.getFaceCount();
743
744 const int numEdgeFunctions = polyOrder * numEdges;
745 const int numFaceFunctions = polyOrder * (polyOrder-1) * numFaces; // 4 faces; 2 families, each with p*(p-1)/2 functions per face
746 const int numInteriorFunctionsPerFamily = (polyOrder > 2) ? (polyOrder-2)*(polyOrder-1)*polyOrder/6 : 0; // (p choose 3)
747 const int numInteriorFunctions = numInteriorFunctionsPerFamily * 3; // 3 families of interior functions
748 this->basisCardinality_ = numEdgeFunctions + numFaceFunctions + numInteriorFunctions;
749 this->basisDegree_ = polyOrder;
750
751 this->basisType_ = BASIS_FEM_HIERARCHICAL;
752 this->basisCoordinates_ = COORDINATES_CARTESIAN;
753 this->functionSpace_ = FUNCTION_SPACE_HCURL;
754
755 const int degreeLength = 1;
756 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Hierarchical H(curl) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
757
758 int fieldOrdinalOffset = 0;
759 // **** vertex functions **** //
760 // no vertex functions in H(curl)
761
762 // **** edge functions **** //
763 const int numFunctionsPerEdge = polyOrder; // p functions associated with each edge
764 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
765 {
766 for (int i=0; i<numFunctionsPerEdge; i++)
767 {
768 this->fieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+1; // the multiplicands involving the gradients of the vertex functions are first degree polynomials; hence the +1 (the remaining multiplicands are order i = 0,…,p-1).
769 }
770 fieldOrdinalOffset += numFunctionsPerEdge;
771 }
772 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != numEdgeFunctions, std::invalid_argument, "Internal error: basis enumeration is incorrect");
773
774 // **** face functions **** //
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++)
779 {
780 for (int faceFamilyOrdinal=1; faceFamilyOrdinal<=numFaceFamilies; faceFamilyOrdinal++)
781 {
782 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
783 int fieldOrdinal = faceFieldOrdinalOffset + faceFamilyOrdinal - 1;
784 for (int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
785 {
786 for (int i=0; i<ij_sum; i++)
787 {
788 this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ij_sum+1;
789 fieldOrdinal += numFaceFamilies; // increment due to the interleaving.
790 }
791 }
792 fieldOrdinalOffset = fieldOrdinal - numFaceFamilies + 1; // due to the interleaving increment, we've gone numFaceFamilies past the last face ordinal. Set offset to be one past.
793 }
794 faceFieldOrdinalOffset += numFaceFunctions / numFaces;
795 }
796 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != numEdgeFunctions + numFaceFunctions, std::invalid_argument, "Internal error: basis enumeration is incorrect");
797
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++)
803 {
804 // following ESEAS, we interleave the interior families. This groups all the interior dofs of a given degree together.
805 int fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
806 for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
807 {
808 for (int i=0; i<ijk_sum-1; i++)
809 {
810 for (int j=1; j<ijk_sum-i; j++)
811 {
812 this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ijk_sum+1;
813 fieldOrdinal += numInteriorFamilies; // increment due to the interleaving.
814 }
815 }
816 }
817 fieldOrdinalOffset = fieldOrdinal - numInteriorFamilies + 1; // due to the interleaving increment, we've gone numFaceFamilies past the last face ordinal. Set offset to be one past.
818 }
819
820 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
821
822 // initialize tags
823 {
824 // ESEAS numbers tetrahedron faces differently from Intrepid2
825 // ESEAS: 012, 013, 123, 023
826 // Intrepid2: 013, 123, 032, 021
827 const int intrepid2FaceOrdinals[4] {3,0,1,2}; // index is the ESEAS face ordinal; value is the intrepid2 ordinal
828
829 const auto & cardinality = this->basisCardinality_;
830
831 // Basis-dependent initializations
832 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
833 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
834 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
835 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
836
837 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
838 const ordinal_type edgeDim = 1, faceDim = 2, volumeDim = 3;
839
840 if (useCGBasis) {
841 {
842 int tagNumber = 0;
843 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
844 {
845 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
846 {
847 tagView(tagNumber*tagSize+0) = edgeDim; // edge dimension
848 tagView(tagNumber*tagSize+1) = edgeOrdinal; // edge id
849 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
850 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge; // total number of dofs on this edge
851 tagNumber++;
852 }
853 }
854 const int numFunctionsPerFace = numFaceFunctions / numFaces;
855 for (int faceOrdinalESEAS=0; faceOrdinalESEAS<numFaces; faceOrdinalESEAS++)
856 {
857 int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
858 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
859 {
860 tagView(tagNumber*tagSize+0) = faceDim; // face dimension
861 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2; // face id
862 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
863 tagView(tagNumber*tagSize+3) = numFunctionsPerFace; // total number of dofs on this face
864 tagNumber++;
865 }
866 }
867
868 // interior
869 for (int functionOrdinal=0; functionOrdinal<numInteriorFunctions; functionOrdinal++)
870 {
871 tagView(tagNumber*tagSize+0) = volumeDim; // interior dimension
872 tagView(tagNumber*tagSize+1) = 0; // volume id
873 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
874 tagView(tagNumber*tagSize+3) = numInteriorFunctions; // total number of interior dofs
875 tagNumber++;
876 }
877 }
878 }
879 else
880 {
881 // DG basis: all functions are associated with interior
882 for (ordinal_type i=0;i<cardinality;++i) {
883 tagView(i*tagSize+0) = volumeDim; // face dimension
884 tagView(i*tagSize+1) = 0; // interior/volume id
885 tagView(i*tagSize+2) = i; // local dof id
886 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
887 }
888 }
889
890 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
891 // tags are constructed on host
893 this->ordinalToTag_,
894 tagView,
895 this->basisCardinality_,
896 tagSize,
897 posScDim,
898 posScOrd,
899 posDfOrd);
900 }
901 }
902
907 const char* getName() const override {
908 return "Intrepid2_HierarchicalBasis_HCURL_TET";
909 }
910
913 virtual bool requireOrientation() const override {
914 return (this->getDegree() > 2);
915 }
916
917 // since the getValues() below only overrides the FEM variant, we specify that
918 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
919 // (It's an error to use the FVD variant on this basis.)
921
940 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
941 const EOperator operatorType = OPERATOR_VALUE ) const override
942 {
943 auto numPoints = inputPoints.extent_int(0);
944
946
947 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
948
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; // because of the way the basis functions are computed, we don't have a second level of parallelism...
953
954 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
955 Kokkos::parallel_for("Hierarchical_HCURL_TET_Functor", policy , functor);
956 }
957
967 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
968 {
971 if (subCellDim == 1)
972 {
973 return Teuchos::rcp(new HVOL_Line(this->basisDegree_-1));
974 }
975 else if (subCellDim == 2)
976 {
977 return Teuchos::rcp(new HCURL_Tri(this->basisDegree_));
978 }
979 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
980 }
981
987 getHostBasis() const override {
988 using HostDeviceType = typename Kokkos::HostSpace::device_type;
990 return Teuchos::rcp( new HostBasisType(polyOrder_) );
991 }
992 };
993} // end namespace Intrepid2
994
995#endif /* Intrepid2_HierarchicalBasis_HCURL_TET_h */
Header file for the abstract base class Intrepid2::Basis.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
H(curl) basis on the triangle using a construction involving Legendre and integrated Jacobi polynomia...
H(vol) basis on the line based on Legendre polynomials.
Free functions, callable from device code, that implement various polynomials useful in basis definit...
Header function for Intrepid2::Util class and other utility functions.
KOKKOS_INLINE_FUNCTION constexpr unsigned getScalarDimensionForView(const ViewType &view)
Returns the size of the Scalar dimension for the View. This is 0 for non-AD types....
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
EBasis basisType_
Type of the basis.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
ordinal_type getDegree() const
Returns the degree of the basis.
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< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
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.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
EFunctionSpace functionSpace_
The function space in which the basis is defined.
For mathematical details of the construction, see:
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
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...
virtual bool requireOrientation() const override
True if orientation is required.
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 > OrdinalTypeArray2DHost
View type for 2d host array.
const char * getName() const override
Returns basis name.
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.
For mathematical details of the construction, see:
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line.
Functor for computing values for the HierarchicalBasis_HCURL_TET class.