Intrepid2
Intrepid2_HierarchicalBasis_HDIV_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_HDIV_TET_h
49#define Intrepid2_HierarchicalBasis_HDIV_TET_h
50
51#include <Kokkos_DynRankView.hpp>
52
53#include <Intrepid2_config.h>
54
55#include "Intrepid2_Basis.hpp"
58#include "Intrepid2_Utils.hpp"
59
60namespace Intrepid2
61{
67 template<class DeviceType, class OutputScalar, class PointScalar,
68 class OutputFieldType, class InputPointsType>
70 {
71 using ExecutionSpace = typename DeviceType::execution_space;
72 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
73 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
74 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
75
76 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
77 using TeamMember = typename TeamPolicy::member_type;
78
79 EOperator opType_;
80
81 OutputFieldType output_; // F,P
82 InputPointsType inputPoints_; // P,D
83
84 ordinal_type polyOrder_;
85 ordinal_type numFields_, numPoints_;
86
87 size_t fad_size_output_;
88
89 static constexpr ordinal_type numVertices = 4;
90 static constexpr ordinal_type numFaces = 4;
91 static constexpr ordinal_type numVerticesPerFace = 3;
92 static constexpr ordinal_type numInteriorFamilies = 3;
93
94 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
95 const ordinal_type face_vertices[numFaces*numVerticesPerFace] = {0,1,2, // face 0
96 0,1,3, // face 1
97 1,2,3, // face 2
98 0,2,3 // face 3
99 };
100
101 const ordinal_type numFaceFunctionsPerFace_;
102 const ordinal_type numFaceFunctions_;
103 const ordinal_type numInteriorFunctionsPerFamily_;
104 const ordinal_type numInteriorFunctions_;
105
106 // interior basis functions are computed in terms of certain face basis functions.
107 const ordinal_type faceOrdinalForInterior_[numInteriorFamilies] = {0,2,3};
108 const ordinal_type faceFamilyForInterior_[numInteriorFamilies] = {0,0,1};
109 const ordinal_type interiorCoordinateOrdinal_[numInteriorFamilies] = {3,0,1}; // m, where V^b_{ijk} is computed in terms of [L^{2(i+j)}_k](1-lambda_m, lambda_m)
110
111 //
112 const ordinal_type interior_face_family_start_ [numInteriorFamilies] = {0,0,1};
113 const ordinal_type interior_face_family_middle_[numInteriorFamilies] = {1,1,2};
114 const ordinal_type interior_face_family_end_ [numInteriorFamilies] = {2,2,0};
115
116 KOKKOS_INLINE_FUNCTION
117 ordinal_type dofOrdinalForFace(const ordinal_type &faceOrdinal,
118 const ordinal_type &zeroBasedFaceFamily,
119 const ordinal_type &i,
120 const ordinal_type &j) const
121 {
122 // determine where the functions for this face start
123 const ordinal_type faceDofOffset = faceOrdinal * numFaceFunctionsPerFace_;
124
125 // rather than deriving a closed formula in terms of i and j (which is potentially error-prone),
126 // we simply step through a for loop much as we do in the basis computations themselves. (This method
127 // is not expected to be called so much as to be worth optimizing.)
128
129 const ordinal_type max_ij_sum = polyOrder_ - 1;
130
131 ordinal_type fieldOrdinal = faceDofOffset + zeroBasedFaceFamily; // families are interleaved on the face.
132
133 for (ordinal_type ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
134 {
135 for (ordinal_type ii=0; ii<ij_sum; ii++)
136 {
137 // j will be ij_sum - i; j >= 1.
138 const ordinal_type jj = ij_sum - ii; // jj >= 1
139 if ( (ii == i) && (jj == j))
140 {
141 // have reached the (i,j) we're looking for
142 return fieldOrdinal;
143 }
144 fieldOrdinal++;
145 }
146 }
147 return -1; // error: not found.
148 }
149
150 Hierarchical_HDIV_TET_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, int polyOrder)
151 : opType_(opType), output_(output), inputPoints_(inputPoints),
152 polyOrder_(polyOrder),
153 fad_size_output_(getScalarDimensionForView(output)),
154 numFaceFunctionsPerFace_(polyOrder * (polyOrder+1)/2), // p*(p+1)/2 functions per face
155 numFaceFunctions_(numFaceFunctionsPerFace_*numFaces), // 4 faces
156 numInteriorFunctionsPerFamily_((polyOrder-1)*polyOrder*(polyOrder+1)/6), // (p+1) choose 3
157 numInteriorFunctions_(numInteriorFunctionsPerFamily_ * numInteriorFamilies) // 3 families of interior functions
158 {
159 numFields_ = output.extent_int(0);
160 numPoints_ = output.extent_int(1);
161
162 const ordinal_type expectedCardinality = numFaceFunctions_ + numInteriorFunctions_;
163
164 // interior family I: computed in terms of face 012 (face ordinal 0), ordinal 0 in face family I.
165 // interior family II: computed in terms of face 123 (face ordinal 2), ordinal 2 in face family I.
166 // interior family III: computed in terms of face 230 (face ordinal 3), ordinal 3 in face family II.
167
168 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
169 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument, "output field size does not match basis cardinality");
170 }
171
172 KOKKOS_INLINE_FUNCTION
173 void computeFaceJacobi(OutputScratchView &P_2ip1,
174 const ordinal_type &zeroBasedFaceOrdinal,
175 const ordinal_type &i,
176 const PointScalar* lambda) const
177 {
178 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
179 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
180 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
181 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];
182
183 const auto & s0 = lambda[s0_index];
184 const auto & s1 = lambda[s1_index];
185 const auto & s2 = lambda[s2_index];
186 const PointScalar jacobiScaling = s0 + s1 + s2;
187
188 const double alpha = i*2.0 + 1;
189 Polynomials::shiftedScaledJacobiValues(P_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);
190 }
191
193 KOKKOS_INLINE_FUNCTION
194 void computeFaceJacobiForInterior(OutputScratchView &P_2ip1,
195 const ordinal_type &zeroBasedInteriorFamilyOrdinal,
196 const ordinal_type &i,
197 const PointScalar* lambda) const
198 {
199 const ordinal_type & relatedFaceOrdinal = faceOrdinalForInterior_[zeroBasedInteriorFamilyOrdinal];
200 const auto &s0_vertex_number = interior_face_family_start_ [zeroBasedInteriorFamilyOrdinal];
201 const auto &s1_vertex_number = interior_face_family_middle_[zeroBasedInteriorFamilyOrdinal];
202 const auto &s2_vertex_number = interior_face_family_end_ [zeroBasedInteriorFamilyOrdinal];
203
204 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
205 const auto &s0_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s0_vertex_number];
206 const auto &s1_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s1_vertex_number];
207 const auto &s2_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s2_vertex_number];
208
209 const auto & s0 = lambda[s0_index];
210 const auto & s1 = lambda[s1_index];
211 const auto & s2 = lambda[s2_index];
212 const PointScalar jacobiScaling = s0 + s1 + s2;
213
214 const double alpha = i*2.0 + 1;
215 Polynomials::shiftedScaledJacobiValues(P_2ip1, alpha, polyOrder_-1, s2, jacobiScaling);
216 }
217
218 KOKKOS_INLINE_FUNCTION
219 void computeFaceLegendre(OutputScratchView &P,
220 const ordinal_type &zeroBasedFaceOrdinal,
221 const PointScalar* lambda) const
222 {
223 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
224 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
225 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
226
227 const auto & s0 = lambda[s0_index];
228 const auto & s1 = lambda[s1_index];
229 const PointScalar legendreScaling = s0 + s1;
230
231 Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, s1, legendreScaling);
232 }
233
234 KOKKOS_INLINE_FUNCTION
235 void computeFaceLegendreForInterior(OutputScratchView &P,
236 const ordinal_type &zeroBasedInteriorFamilyOrdinal,
237 const PointScalar* lambda) const
238 {
239 const ordinal_type & relatedFaceOrdinal = faceOrdinalForInterior_[zeroBasedInteriorFamilyOrdinal];
240 const auto &s0_vertex_number = interior_face_family_start_ [zeroBasedInteriorFamilyOrdinal];
241 const auto &s1_vertex_number = interior_face_family_middle_[zeroBasedInteriorFamilyOrdinal];
242
243 // index into face_vertices with faceOrdinal * numVerticesPerFace + vertexNumber
244 const auto &s0_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s0_vertex_number];
245 const auto &s1_index = face_vertices[relatedFaceOrdinal * numVerticesPerFace + s1_vertex_number];
246
247 const auto & s0 = lambda[s0_index];
248 const auto & s1 = lambda[s1_index];
249 const PointScalar legendreScaling = s0 + s1;
250
251 Polynomials::shiftedScaledLegendreValues(P, polyOrder_-1, s1, legendreScaling);
252 }
253
254 KOKKOS_INLINE_FUNCTION
255 void computeFaceVectorWeight(OutputScalar &vectorWeight_x,
256 OutputScalar &vectorWeight_y,
257 OutputScalar &vectorWeight_z,
258 const ordinal_type &zeroBasedFaceOrdinal,
259 const PointScalar* lambda,
260 const PointScalar* lambda_dx,
261 const PointScalar* lambda_dy,
262 const PointScalar* lambda_dz) const
263 {
264 // compute s0 (grad s1 x grad s2) + s1 (grad s2 x grad s0) + s2 (grad s0 x grad s1)
265
266 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
267 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
268 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];
269
270 const auto & s0 = lambda [s0_index];
271 const auto & s0_dx = lambda_dx[s0_index];
272 const auto & s0_dy = lambda_dy[s0_index];
273 const auto & s0_dz = lambda_dz[s0_index];
274
275 const auto & s1 = lambda [s1_index];
276 const auto & s1_dx = lambda_dx[s1_index];
277 const auto & s1_dy = lambda_dy[s1_index];
278 const auto & s1_dz = lambda_dz[s1_index];
279
280 const auto & s2 = lambda [s2_index];
281 const auto & s2_dx = lambda_dx[s2_index];
282 const auto & s2_dy = lambda_dy[s2_index];
283 const auto & s2_dz = lambda_dz[s2_index];
284
285 vectorWeight_x = s0 * (s1_dy * s2_dz - s1_dz * s2_dy)
286 + s1 * (s2_dy * s0_dz - s2_dz * s0_dy)
287 + s2 * (s0_dy * s1_dz - s0_dz * s1_dy);
288
289 vectorWeight_y = s0 * (s1_dz * s2_dx - s1_dx * s2_dz)
290 + s1 * (s2_dz * s0_dx - s2_dx * s0_dz)
291 + s2 * (s0_dz * s1_dx - s0_dx * s1_dz);
292
293 vectorWeight_z = s0 * (s1_dx * s2_dy - s1_dy * s2_dx)
294 + s1 * (s2_dx * s0_dy - s2_dy * s0_dx)
295 + s2 * (s0_dx * s1_dy - s0_dy * s1_dx);
296 }
297
298 // This is the "Ancillary Operator" V^{tri}_{ij} on p. 433 of Fuentes et al.
299 KOKKOS_INLINE_FUNCTION
300 void faceFunctionValue(OutputScalar &value_x,
301 OutputScalar &value_y,
302 OutputScalar &value_z,
303 const ordinal_type &i, // i >= 0
304 const ordinal_type &j, // j >= 0
305 const OutputScratchView &P, // container in which shiftedScaledLegendreValues have been computed for the appropriate face
306 const OutputScratchView &P_2ip1, // container in which shiftedScaledJacobiValues have been computed for (2i+1) for the appropriate face
307 const OutputScalar &vectorWeight_x, // x component of s0 (grad s1 x grad s2) + s1 (grad s2 x grad s0) + s2 (grad s0 x grad s1)
308 const OutputScalar &vectorWeight_y, // y component
309 const OutputScalar &vectorWeight_z, // z component
310 const PointScalar* lambda) const
311 {
312 const auto &P_i = P(i);
313 const auto &P_2ip1_j = P_2ip1(j);
314
315 value_x = P_i * P_2ip1_j * vectorWeight_x;
316 value_y = P_i * P_2ip1_j * vectorWeight_y;
317 value_z = P_i * P_2ip1_j * vectorWeight_z;
318 }
319
320 KOKKOS_INLINE_FUNCTION
321 void computeFaceDivWeight(OutputScalar &divWeight,
322 const ordinal_type &zeroBasedFaceOrdinal,
323 const PointScalar* lambda_dx,
324 const PointScalar* lambda_dy,
325 const PointScalar* lambda_dz) const
326 {
327 // grad s0 \dot (grad s1 x grad s2)
328
329 const auto &s0_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 0];
330 const auto &s1_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 1];
331 const auto &s2_index = face_vertices[zeroBasedFaceOrdinal * numVerticesPerFace + 2];
332
333 const auto & s0_dx = lambda_dx[s0_index];
334 const auto & s0_dy = lambda_dy[s0_index];
335 const auto & s0_dz = lambda_dz[s0_index];
336
337 const auto & s1_dx = lambda_dx[s1_index];
338 const auto & s1_dy = lambda_dy[s1_index];
339 const auto & s1_dz = lambda_dz[s1_index];
340
341 const auto & s2_dx = lambda_dx[s2_index];
342 const auto & s2_dy = lambda_dy[s2_index];
343 const auto & s2_dz = lambda_dz[s2_index];
344
345 divWeight = s0_dx * (s1_dy * s2_dz - s1_dz * s2_dy)
346 + s0_dy * (s1_dz * s2_dx - s1_dx * s2_dz)
347 + s0_dz * (s1_dx * s2_dy - s1_dy * s2_dx);
348 }
349
350 KOKKOS_INLINE_FUNCTION
351 void computeInteriorIntegratedJacobi(OutputScratchView &L_2ipjp1,
352 const ordinal_type &i,
353 const ordinal_type &j,
354 const ordinal_type &zeroBasedFamilyOrdinal,
355 const PointScalar* lambda) const
356 {
357 const auto &lambda_m = lambda[interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal]];
358
359 const double alpha = 2 * (i + j + 1);
360
361 const PointScalar jacobiScaling = 1.0;
362 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ipjp1, alpha, polyOrder_-1, lambda_m, jacobiScaling);
363 }
364
365 KOKKOS_INLINE_FUNCTION
366 void computeInteriorJacobi(OutputScratchView &P_2ipjp1,
367 const ordinal_type &i,
368 const ordinal_type &j,
369 const ordinal_type &zeroBasedFamilyOrdinal,
370 const PointScalar* lambda) const
371 {
372 const auto &lambda_m = lambda[interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal]];
373
374 const double alpha = 2 * (i + j + 1);
375
376 const PointScalar jacobiScaling = 1.0;
377 Polynomials::shiftedScaledJacobiValues(P_2ipjp1, alpha, polyOrder_-1, lambda_m, jacobiScaling);
378 }
379
380 // Divergence of the "Ancillary Operator" V^{tri}_{ij} on p. 433 of Fuentes et al.
381 KOKKOS_INLINE_FUNCTION
382 void faceFunctionDiv(OutputScalar &divValue,
383 const ordinal_type &i, // i >= 0
384 const ordinal_type &j, // j >= 0
385 const OutputScratchView &P, // container in which shiftedScaledLegendreValues have been computed for the appropriate face
386 const OutputScratchView &P_2ip1, // container in which shiftedScaledJacobiValues have been computed for (2i+1) for the appropriate face
387 const OutputScalar &divWeight, // grad s0 \dot (grad s1 x grad s2)
388 const PointScalar* lambda) const
389 {
390 const auto &P_i = P(i);
391 const auto &P_2ip1_j = P_2ip1(j);
392
393 divValue = (i + j + 3.) * P_i * P_2ip1_j * divWeight;
394 }
395
396 // grad ([L^{2(i+j+1)}_k](1-lambda_m,lambda_m)), used in divergence of interior basis functions
397 KOKKOS_INLINE_FUNCTION
398 void gradInteriorIntegratedJacobi(OutputScalar &L_2ipjp1_dx,
399 OutputScalar &L_2ipjp1_dy,
400 OutputScalar &L_2ipjp1_dz,
401 const ordinal_type &zeroBasedFamilyOrdinal,
402 const ordinal_type &j,
403 const ordinal_type &k,
404 const OutputScratchView &P_2ipjp1, // container in which shiftedScaledJacobiValues have been computed for alpha=2(i+j+1), t0=1-lambda_m, t1=lambda_m
405 const PointScalar* lambda,
406 const PointScalar* lambda_dx,
407 const PointScalar* lambda_dy,
408 const PointScalar* lambda_dz) const
409 {
410 // grad [L^alpha_k](t0,t1) = [P^alpha_{k-1}](t0,t1) grad(t1) + [R^alpha_{k-1}](t0,t1) grad(t1+t0)
411 // here, t0 = 1-lambda_m, t1 = lambda_m ==> t1 + t0 = 1 ==> grad(t1+t0) = 0 ==> the R term vanishes.
412
413 const ordinal_type &m = interiorCoordinateOrdinal_[zeroBasedFamilyOrdinal];
414
415 L_2ipjp1_dx = P_2ipjp1(k-1) * lambda_dx[m];
416 L_2ipjp1_dy = P_2ipjp1(k-1) * lambda_dy[m];
417 L_2ipjp1_dz = P_2ipjp1(k-1) * lambda_dz[m];
418 }
419
420 KOKKOS_INLINE_FUNCTION
421 void interiorFunctionDiv(OutputScalar &outputDiv,
422 OutputScalar &L_2ipjp1_k,
423 OutputScalar &faceDiv,
424 OutputScalar &L_2ipjp1_k_dx,
425 OutputScalar &L_2ipjp1_k_dy,
426 OutputScalar &L_2ipjp1_k_dz,
427 OutputScalar &faceValue_x,
428 OutputScalar &faceValue_y,
429 OutputScalar &faceValue_z) const
430 {
431 outputDiv = L_2ipjp1_k * faceDiv + L_2ipjp1_k_dx * faceValue_x + L_2ipjp1_k_dy * faceValue_y + L_2ipjp1_k_dz * faceValue_z;
432 }
433
434 KOKKOS_INLINE_FUNCTION
435 void operator()(const TeamMember &teamMember) const {
436 auto pointOrdinal = teamMember.league_rank();
437 OutputScratchView scratch0, scratch1, scratch2, scratch3;
438 if (fad_size_output_ > 0) {
439 scratch0 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
440 scratch1 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
441 scratch2 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
442 scratch3 = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
443 }
444 else {
445 scratch0 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
446 scratch1 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
447 scratch2 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
448 scratch3 = OutputScratchView(teamMember.team_shmem(), polyOrder_);
449 }
450
451 const auto & x = inputPoints_(pointOrdinal,0);
452 const auto & y = inputPoints_(pointOrdinal,1);
453 const auto & z = inputPoints_(pointOrdinal,2);
454
455 // write as barycentric coordinates:
456 const PointScalar lambda[4] = {1. - x - y - z, x, y, z};
457 const PointScalar lambda_dx[4] = {-1., 1., 0., 0.};
458 const PointScalar lambda_dy[4] = {-1., 0., 1., 0.};
459 const PointScalar lambda_dz[4] = {-1., 0., 0., 1.};
460
461 switch (opType_)
462 {
463 case OPERATOR_VALUE:
464 {
465 // face functions
466 {
467 // relabel scratch views
468 auto &scratchP = scratch0;
469 auto &scratchP_2ip1 = scratch1;
470
471 const ordinal_type max_ij_sum = polyOrder_ - 1;
472
473 for (ordinal_type faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
474 {
475 OutputScalar divWeight;
476 computeFaceDivWeight(divWeight, faceOrdinal, lambda_dx, lambda_dy, lambda_dz);
477
478 OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;
479 computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, faceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);
480
481 ordinal_type fieldOrdinal = faceOrdinal * numFaceFunctionsPerFace_;
482 computeFaceLegendre(scratchP, faceOrdinal, lambda);
483
484 for (int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)
485 {
486 for (int i=0; i<=ij_sum; i++)
487 {
488 computeFaceJacobi(scratchP_2ip1, faceOrdinal, i, lambda);
489
490 const int j = ij_sum - i; // j >= 1
491
492 auto & output_x = output_(fieldOrdinal,pointOrdinal,0);
493 auto & output_y = output_(fieldOrdinal,pointOrdinal,1);
494 auto & output_z = output_(fieldOrdinal,pointOrdinal,2);
495
496 faceFunctionValue(output_x, output_y, output_z, i, j,
497 scratchP, scratchP_2ip1, vectorWeight_x,
498 vectorWeight_y, vectorWeight_z, lambda);
499
500 fieldOrdinal++;
501 } // i
502 } // ij_sum
503 } // faceOrdinal
504 } // face functions block
505
506 // interior functions
507 {
508 // relabel scratch views
509 auto &scratchP = scratch0;
510 auto &scratchP_2ip1 = scratch1;
511 auto &scratchL_2ipjp1 = scratch2; // L^{2(i+j+1)}, integrated Jacobi
512
513 const ordinal_type min_ijk_sum = 1;
514 const ordinal_type max_ijk_sum = polyOrder_-1;
515 const ordinal_type min_ij_sum = 0;
516 const ordinal_type min_k = 1;
517 const ordinal_type min_j = 0;
518 const ordinal_type min_i = 0;
519
520 OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;
521
522 for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
523 {
524 // following ESEAS, we interleave the interior families. This groups all the interior dofs of a given degree together.
525
526 ordinal_type interiorFamilyFieldOrdinal =
527 numFaceFunctions_ + interiorFamilyOrdinal - 1;
528
529 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
530
531 computeFaceLegendreForInterior(scratchP,
532 interiorFamilyOrdinal - 1, lambda);
533 computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, relatedFaceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);
534
535 for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
536 {
537 for (int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)
538 {
539 for (int i=min_i; i<=ij_sum-min_j; i++)
540 {
541 const ordinal_type j = ij_sum-i;
542 const ordinal_type k = ijk_sum - ij_sum;
543
545 scratchP_2ip1, interiorFamilyOrdinal - 1, i, lambda);
546 computeInteriorIntegratedJacobi(scratchL_2ipjp1, i, j,
547 interiorFamilyOrdinal - 1,
548 lambda);
549
550 OutputScalar V_x, V_y, V_z;
551
552 faceFunctionValue(V_x, V_y, V_z, i, j, scratchP,
553 scratchP_2ip1, vectorWeight_x,
554 vectorWeight_y, vectorWeight_z, lambda);
555
556 auto &output_x =
557 output_(interiorFamilyFieldOrdinal, pointOrdinal, 0);
558 auto &output_y =
559 output_(interiorFamilyFieldOrdinal, pointOrdinal, 1);
560 auto &output_z =
561 output_(interiorFamilyFieldOrdinal, pointOrdinal, 2);
562
563 output_x = V_x * scratchL_2ipjp1(k);
564 output_y = V_y * scratchL_2ipjp1(k);
565 output_z = V_z * scratchL_2ipjp1(k);
566
567 interiorFamilyFieldOrdinal +=
568 numInteriorFamilies; // increment due to the
569 // interleaving.
570 }
571 }
572 }
573 }
574 } // interior functions block
575
576 } // end OPERATOR_VALUE
577 break;
578 case OPERATOR_DIV:
579 {
580 // rename the scratch memory to match our usage here:
581 auto &scratchP = scratch0;
582 auto &scratchP_2ip1 = scratch1;
583
584 // following ESEAS, we interleave the face families. This groups all the face dofs of a given degree together.
585 ordinal_type fieldOrdinal = 0;
586 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
587 {
588 const int max_ij_sum = polyOrder_ - 1;
589 computeFaceLegendre(scratchP, faceOrdinal, lambda);
590 OutputScalar divWeight;
591 computeFaceDivWeight(divWeight, faceOrdinal, lambda_dx, lambda_dy, lambda_dz);
592 for (int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)
593 {
594 for (int i=0; i<=ij_sum; i++)
595 {
596 const int j = ij_sum - i; // j >= 0
597
598 computeFaceJacobi(scratchP_2ip1, faceOrdinal, i, lambda);
599 auto &outputValue = output_(fieldOrdinal,pointOrdinal);
600 faceFunctionDiv(outputValue, i, j, scratchP, scratchP_2ip1,
601 divWeight, lambda);
602
603 fieldOrdinal++;
604 } // i
605 } // ij_sum
606 } // faceOrdinal
607
608 // interior functions
609 {
610 // rename the scratch memory to match our usage here:
611 auto &scratchL_2ipjp1 = scratch2;
612 auto &scratchP_2ipjp1 = scratch3;
613
614 const int interiorFieldOrdinalOffset = numFaceFunctions_;
615 for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
616 {
617 // following ESEAS, we interleave the interior families. This groups all the interior dofs of a given degree together.
618
619 const ordinal_type relatedFaceOrdinal = faceOrdinalForInterior_[interiorFamilyOrdinal-1];
620
621 computeFaceLegendreForInterior(scratchP,
622 interiorFamilyOrdinal - 1, lambda);
623 OutputScalar divWeight;
624 computeFaceDivWeight(divWeight, relatedFaceOrdinal, lambda_dx, lambda_dy, lambda_dz);
625
626 OutputScalar vectorWeight_x, vectorWeight_y, vectorWeight_z;
627 computeFaceVectorWeight(vectorWeight_x, vectorWeight_y, vectorWeight_z, relatedFaceOrdinal, lambda, lambda_dx, lambda_dy, lambda_dz);
628
629 ordinal_type interiorFieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
630
631 const ordinal_type min_ijk_sum = 1;
632 const ordinal_type max_ijk_sum = polyOrder_-1;
633 const ordinal_type min_ij_sum = 0;
634 const ordinal_type min_k = 1;
635 const ordinal_type min_j = 0;
636 const ordinal_type min_i = 0;
637 for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
638 {
639 for (int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)
640 {
641 for (int i=min_i; i<=ij_sum-min_j; i++)
642 {
643 const ordinal_type j = ij_sum-i;
644 const ordinal_type k = ijk_sum - ij_sum;
646 scratchP_2ip1, interiorFamilyOrdinal - 1, i, lambda);
647
648 OutputScalar faceDiv;
649 faceFunctionDiv(faceDiv, i, j, scratchP, scratchP_2ip1,
650 divWeight, lambda);
651
652 OutputScalar faceValue_x, faceValue_y, faceValue_z;
653
654 faceFunctionValue(faceValue_x, faceValue_y, faceValue_z, i,
655 j, scratchP, scratchP_2ip1,
656 vectorWeight_x, vectorWeight_y,
657 vectorWeight_z, lambda);
658 computeInteriorJacobi(scratchP_2ipjp1, i, j,
659 interiorFamilyOrdinal - 1, lambda);
660
661 computeInteriorIntegratedJacobi(scratchL_2ipjp1, i, j,
662 interiorFamilyOrdinal - 1,
663 lambda);
664
665 OutputScalar L_2ipjp1_k_dx, L_2ipjp1_k_dy, L_2ipjp1_k_dz;
666 gradInteriorIntegratedJacobi(
667 L_2ipjp1_k_dx, L_2ipjp1_k_dy, L_2ipjp1_k_dz,
668 interiorFamilyOrdinal - 1, j, k, scratchP_2ipjp1,
669 lambda, lambda_dx, lambda_dy, lambda_dz);
670
671 auto &outputDiv = output_(interiorFieldOrdinal, pointOrdinal);
672 interiorFunctionDiv(outputDiv, scratchL_2ipjp1(k), faceDiv,
673 L_2ipjp1_k_dx, L_2ipjp1_k_dy,
674 L_2ipjp1_k_dz, faceValue_x, faceValue_y,
675 faceValue_z);
676
677 interiorFieldOrdinal += numInteriorFamilies; // increment due to the interleaving.
678 }
679 }
680 }
681 }
682 } // interior functions block
683 } // OPERATOR_DIV
684 break;
685 case OPERATOR_GRAD:
686 case OPERATOR_D1:
687 case OPERATOR_D2:
688 case OPERATOR_D3:
689 case OPERATOR_D4:
690 case OPERATOR_D5:
691 case OPERATOR_D6:
692 case OPERATOR_D7:
693 case OPERATOR_D8:
694 case OPERATOR_D9:
695 case OPERATOR_D10:
696 INTREPID2_TEST_FOR_ABORT( true,
697 ">>> ERROR: (Intrepid2::Hierarchical_HDIV_TET_Functor) Unsupported differential operator");
698 default:
699 // unsupported operator type
700 device_assert(false);
701 }
702 }
703
704 // Provide the shared memory capacity.
705 // This function takes the team_size as an argument,
706 // which allows team_size-dependent allocations.
707 size_t team_shmem_size (int team_size) const
708 {
709 // we will use shared memory to create a fast buffer for basis computations
710 size_t shmem_size = 0;
711 if (fad_size_output_ > 0)
712 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
713 else
714 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
715
716 return shmem_size;
717 }
718 };
719
728 template<typename DeviceType,
729 typename OutputScalar = double,
730 typename PointScalar = double,
731 bool useCGBasis = true> // if useCGBasis is false, all basis functions will be associated with the interior
733 : public Basis<DeviceType,OutputScalar,PointScalar>
734 {
735 public:
737
740
741 using typename BasisBase::OutputViewType;
742 using typename BasisBase::PointViewType;
743 using typename BasisBase::ScalarViewType;
744
745 using typename BasisBase::ExecutionSpace;
746
747 protected:
748 int polyOrder_; // the maximum order of the polynomial
749 public:
754 HierarchicalBasis_HDIV_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
755 :
756 polyOrder_(polyOrder)
757 {
758
759 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
760 const int numFaces = this->basisCellTopology_.getFaceCount();
761
762 const int numVertexFunctions = 0;
763 const int numEdgeFunctions = 0;
764 const int numFaceFunctions = numFaces * polyOrder * (polyOrder+1) / 2; // 4 faces; 2 families, each with p*(p-1)/2 functions per face
765 const int numInteriorFunctionsPerFamily = (polyOrder-1)*polyOrder*(polyOrder+1)/6;
766 const int numInteriorFunctions = numInteriorFunctionsPerFamily * 3; // 3 families of interior functions
767 this->basisCardinality_ = numVertexFunctions + numEdgeFunctions + numFaceFunctions + numInteriorFunctions;
768 this->basisDegree_ = polyOrder;
769
770 this->basisType_ = BASIS_FEM_HIERARCHICAL;
771 this->basisCoordinates_ = COORDINATES_CARTESIAN;
772 this->functionSpace_ = FUNCTION_SPACE_HDIV;
773
774 const int degreeLength = 1;
775 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Hierarchical H(div) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
776
777 // **** vertex functions **** //
778 // no vertex functions in H(div)
779
780 // **** edge functions **** //
781 // no edge functions in H(div)
782
783 // **** face functions **** //
784 const int max_ij_sum = polyOrder-1;
785 ordinal_type fieldOrdinal = 0;
786 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
787 {
788 for (int ij_sum=0; ij_sum <= max_ij_sum; ij_sum++)
789 {
790 for (int i=0; i<=ij_sum; i++)
791 {
792 this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ij_sum+1;
793 fieldOrdinal++;
794 }
795 }
796 }
797 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal != numEdgeFunctions + numFaceFunctions, std::invalid_argument, "Internal error: basis enumeration is incorrect");
798
799 const int numInteriorFamilies = 3;
800 const int interiorFieldOrdinalOffset = fieldOrdinal;
801 const ordinal_type min_ijk_sum = 1;
802 const ordinal_type max_ijk_sum = polyOrder_-1;
803 const ordinal_type min_ij_sum = 0;
804 const ordinal_type min_k = 1;
805 const ordinal_type min_j = 0;
806 const ordinal_type min_i = 0;
807 for (int interiorFamilyOrdinal=1; interiorFamilyOrdinal<=numInteriorFamilies; interiorFamilyOrdinal++)
808 {
809 // following ESEAS, we interleave the interior families. This groups all the interior dofs of a given degree together.
810 fieldOrdinal = interiorFieldOrdinalOffset + interiorFamilyOrdinal - 1;
811 for (int ijk_sum=min_ijk_sum; ijk_sum <= max_ijk_sum; ijk_sum++)
812 {
813 for (int ij_sum=min_ij_sum; ij_sum<=ijk_sum-min_k; ij_sum++)
814 {
815 for (int i=min_i; i<=ij_sum-min_j; i++)
816 {
817 this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = ijk_sum+1;
818 fieldOrdinal += numInteriorFamilies; // increment due to the interleaving.
819 }
820 }
821 }
822 fieldOrdinal = fieldOrdinal - numInteriorFamilies + 1; // due to the interleaving increment, we've gone numInteriorFamilies past the last interior ordinal. Set fieldOrdinal to be one past.
823 }
824
825 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
826
827 // initialize tags
828 {
829 // ESEAS numbers tetrahedron faces differently from Intrepid2
830 // ESEAS: 012, 013, 123, 023
831 // Intrepid2: 013, 123, 032, 021
832 const int intrepid2FaceOrdinals[4] {3,0,1,2}; // index is the ESEAS face ordinal; value is the intrepid2 ordinal
833
834 const auto & cardinality = this->basisCardinality_;
835
836 // Basis-dependent initializations
837 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
838 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
839 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
840 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
841
842 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
843 const ordinal_type faceDim = 2, volumeDim = 3;
844
845 if (useCGBasis) {
846 {
847 int tagNumber = 0;
848 const int numFunctionsPerFace = numFaceFunctions / numFaces;
849 for (int faceOrdinalESEAS=0; faceOrdinalESEAS<numFaces; faceOrdinalESEAS++)
850 {
851 int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
852 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
853 {
854 tagView(tagNumber*tagSize+0) = faceDim; // face dimension
855 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2; // face id
856 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
857 tagView(tagNumber*tagSize+3) = numFunctionsPerFace; // total number of dofs on this face
858 tagNumber++;
859 }
860 }
861
862 // interior
863 for (int functionOrdinal=0; functionOrdinal<numInteriorFunctions; functionOrdinal++)
864 {
865 tagView(tagNumber*tagSize+0) = volumeDim; // interior dimension
866 tagView(tagNumber*tagSize+1) = 0; // volume id
867 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
868 tagView(tagNumber*tagSize+3) = numInteriorFunctions; // total number of interior dofs
869 tagNumber++;
870 }
871 }
872 }
873 else
874 {
875 // DG basis: all functions are associated with interior
876 for (ordinal_type i=0;i<cardinality;++i) {
877 tagView(i*tagSize+0) = volumeDim; // face dimension
878 tagView(i*tagSize+1) = 0; // interior/volume id
879 tagView(i*tagSize+2) = i; // local dof id
880 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
881 }
882 }
883
884 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
885 // tags are constructed on host
887 this->ordinalToTag_,
888 tagView,
889 this->basisCardinality_,
890 tagSize,
891 posScDim,
892 posScOrd,
893 posDfOrd);
894 }
895 }
896
901 const char* getName() const override {
902 return "Intrepid2_HierarchicalBasis_HDIV_TET";
903 }
904
907 virtual bool requireOrientation() const override {
908 return (this->getDegree() > 2);
909 }
910
911 // since the getValues() below only overrides the FEM variant, we specify that
912 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
913 // (It's an error to use the FVD variant on this basis.)
915
934 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
935 const EOperator operatorType = OPERATOR_VALUE ) const override
936 {
937 auto numPoints = inputPoints.extent_int(0);
938
940
941 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
942
943 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
944 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
945 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
946 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
947
948 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
949 Kokkos::parallel_for("Hierarchical_HDIV_TET_Functor", policy , functor);
950 }
951
961 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
962 {
964 if (subCellDim == 2)
965 {
966 return Teuchos::rcp(new HVOL_Tri(this->basisDegree_-1));
967 }
968 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
969 }
970
976 getHostBasis() const override {
977 using HostDeviceType = typename Kokkos::HostSpace::device_type;
979 return Teuchos::rcp( new HostBasisType(polyOrder_) );
980 }
981 };
982} // end namespace Intrepid2
983
984#endif /* Intrepid2_HierarchicalBasis_HDIV_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(vol) basis on the triangle based on integrated 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.
const char * getName() const override
Returns basis name.
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...
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.
virtual bool requireOrientation() const override
True if orientation is required.
HierarchicalBasis_HDIV_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.
Basis defining Legendre basis on the line, a polynomial subspace of H(vol) on the line: extension to ...
Functor for computing values for the HierarchicalBasis_HDIV_TET class.
KOKKOS_INLINE_FUNCTION void computeFaceJacobiForInterior(OutputScratchView &P_2ip1, const ordinal_type &zeroBasedInteriorFamilyOrdinal, const ordinal_type &i, const PointScalar *lambda) const
The face functions we compute for interior blending can have a different orientation than the ones us...