Intrepid2
Intrepid2_IntegratedLegendreBasis_HGRAD_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 Kyungjoo Kim (kyukim@sandia.gov),
38// Mauro Perego (mperego@sandia.gov), or
39// Nate Roberts (nvrober@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h
50#define Intrepid2_IntegratedLegendreBasis_HGRAD_TET_h
51
52#include <Kokkos_DynRankView.hpp>
53
54#include <Intrepid2_config.h>
55
56#include "Intrepid2_Basis.hpp"
60#include "Intrepid2_Utils.hpp"
61
62namespace Intrepid2
63{
69 template<class DeviceType, class OutputScalar, class PointScalar,
70 class OutputFieldType, class InputPointsType>
72 {
73 using ExecutionSpace = typename DeviceType::execution_space;
74 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
75 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
76 using OutputScratchView2D = Kokkos::View<OutputScalar**,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
77 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
78
79 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
80 using TeamMember = typename TeamPolicy::member_type;
81
82 EOperator opType_;
83
84 OutputFieldType output_; // F,P
85 InputPointsType inputPoints_; // P,D
86
87 int polyOrder_;
88 bool defineVertexFunctions_;
89 int numFields_, numPoints_;
90
91 size_t fad_size_output_;
92
93 static const int numVertices = 4;
94 static const int numEdges = 6;
95 // the following ordering of the edges matches that used by ESEAS
96 const int edge_start_[numEdges] = {0,1,0,0,1,2}; // edge i is from edge_start_[i] to edge_end_[i]
97 const int edge_end_[numEdges] = {1,2,2,3,3,3}; // edge i is from edge_start_[i] to edge_end_[i]
98
99 static const int numFaces = 4;
100 const int face_vertex_0[numFaces] = {0,0,1,0}; // faces are abc where 0 ≤ a < b < c ≤ 3
101 const int face_vertex_1[numFaces] = {1,1,2,2};
102 const int face_vertex_2[numFaces] = {2,3,3,3};
103
104 // this allows us to look up the edge ordinal of the first edge of a face
105 // this is useful because face functions are defined using edge basis functions of the first edge of the face
106 const int face_ordinal_of_first_edge[numFaces] = {0,0,1,2};
107
108 Hierarchical_HGRAD_TET_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints,
109 int polyOrder, bool defineVertexFunctions)
110 : opType_(opType), output_(output), inputPoints_(inputPoints),
111 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
112 fad_size_output_(getScalarDimensionForView(output))
113 {
114 numFields_ = output.extent_int(0);
115 numPoints_ = output.extent_int(1);
116 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
117 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)*(polyOrder_+3)/6, std::invalid_argument, "output field size does not match basis cardinality");
118 }
119
120 KOKKOS_INLINE_FUNCTION
121 void operator()( const TeamMember & teamMember ) const
122 {
123 const int numFaceBasisFunctionsPerFace = (polyOrder_-2) * (polyOrder_-1) / 2;
124 const int numInteriorBasisFunctions = (polyOrder_-3) * (polyOrder_-2) * (polyOrder_-1) / 6;
125
126 auto pointOrdinal = teamMember.league_rank();
127 OutputScratchView legendre_values1_at_point, legendre_values2_at_point;
128 OutputScratchView2D jacobi_values1_at_point, jacobi_values2_at_point, jacobi_values3_at_point;
129 const int numAlphaValues = (polyOrder_-1 > 1) ? (polyOrder_-1) : 1; // make numAlphaValues at least 1 so we can avoid zero-extent allocations…
130 if (fad_size_output_ > 0) {
131 legendre_values1_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
132 legendre_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
133 jacobi_values1_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
134 jacobi_values2_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
135 jacobi_values3_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1, fad_size_output_);
136 }
137 else {
138 legendre_values1_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
139 legendre_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
140 jacobi_values1_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
141 jacobi_values2_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
142 jacobi_values3_at_point = OutputScratchView2D(teamMember.team_shmem(), numAlphaValues, polyOrder_ + 1);
143 }
144
145 const auto & x = inputPoints_(pointOrdinal,0);
146 const auto & y = inputPoints_(pointOrdinal,1);
147 const auto & z = inputPoints_(pointOrdinal,2);
148
149 // write as barycentric coordinates:
150 const PointScalar lambda[numVertices] = {1. - x - y - z, x, y, z};
151 const PointScalar lambda_dx[numVertices] = {-1., 1., 0., 0.};
152 const PointScalar lambda_dy[numVertices] = {-1., 0., 1., 0.};
153 const PointScalar lambda_dz[numVertices] = {-1., 0., 0., 1.};
154
155 const int num1DEdgeFunctions = polyOrder_ - 1;
156
157 switch (opType_)
158 {
159 case OPERATOR_VALUE:
160 {
161 // vertex functions come first, according to vertex ordering: (0,0,0), (1,0,0), (0,1,0), (0,0,1)
162 for (int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
163 {
164 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
165 }
166 if (!defineVertexFunctions_)
167 {
168 // "DG" basis case
169 // here, we overwrite the first vertex function with 1:
170 output_(0,pointOrdinal) = 1.0;
171 }
172
173 // edge functions
174 int fieldOrdinalOffset = numVertices;
175 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
176 {
177 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
178 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
179
180 Polynomials::shiftedScaledIntegratedLegendreValues(legendre_values1_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
181 for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
182 {
183 // the first two integrated legendre functions are essentially the vertex functions; hence the +2 on on the RHS here:
184 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = legendre_values1_at_point(edgeFunctionOrdinal+2);
185 }
186 fieldOrdinalOffset += num1DEdgeFunctions;
187 }
188 /*
189 Face functions for face abc are the product of edge functions on their ab edge
190 and a Jacobi polynomial [L^2i_j](s0+s1,s2) = L^2i_j(s2;s0+s1+s2)
191 */
192 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
193 {
194 const auto & s0 = lambda[face_vertex_0[faceOrdinal]];
195 const auto & s1 = lambda[face_vertex_1[faceOrdinal]];
196 const auto & s2 = lambda[face_vertex_2[faceOrdinal]];
197 const PointScalar jacobiScaling = s0 + s1 + s2;
198
199 // compute integrated Jacobi values for each desired value of alpha
200 for (int n=2; n<=polyOrder_; n++)
201 {
202 const double alpha = n*2;
203 const int alphaOrdinal = n-2;
204 using Kokkos::subview;
205 using Kokkos::ALL;
206 auto jacobi_alpha = subview(jacobi_values1_at_point, alphaOrdinal, ALL);
207 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
208 }
209
210 const int edgeOrdinal = face_ordinal_of_first_edge[faceOrdinal];
211 int localFaceBasisOrdinal = 0;
212 for (int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
213 {
214 for (int i=2; i<totalPolyOrder; i++)
215 {
216 const int edgeBasisOrdinal = edgeOrdinal*num1DEdgeFunctions + i-2 + numVertices;
217 const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
218 const int alphaOrdinal = i-2;
219
220 const int j = totalPolyOrder - i;
221 const auto & jacobiValue = jacobi_values1_at_point(alphaOrdinal,j);
222 const int fieldOrdinal = fieldOrdinalOffset + localFaceBasisOrdinal;
223 output_(fieldOrdinal,pointOrdinal) = edgeValue * jacobiValue;
224
225 localFaceBasisOrdinal++;
226 }
227 }
228 fieldOrdinalOffset += numFaceBasisFunctionsPerFace;
229 }
230 // interior functions
231 // compute integrated Jacobi values for each desired value of alpha
232 for (int n=3; n<=polyOrder_; n++)
233 {
234 const double alpha = n*2;
235 const double jacobiScaling = 1.0;
236 const int alphaOrdinal = n-3;
237 using Kokkos::subview;
238 using Kokkos::ALL;
239 auto jacobi_alpha = subview(jacobi_values1_at_point, alphaOrdinal, ALL);
240 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-3, lambda[3], jacobiScaling);
241 }
242 const int min_i = 2;
243 const int min_j = 1;
244 const int min_k = 1;
245 const int min_ij = min_i + min_j;
246 const int min_ijk = min_ij + min_k;
247 int localInteriorBasisOrdinal = 0;
248 for (int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= polyOrder_; totalPolyOrder_ijk++)
249 {
250 int localFaceBasisOrdinal = 0;
251 for (int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
252 {
253 for (int i=2; i <= totalPolyOrder_ij-min_j; i++)
254 {
255 const int j = totalPolyOrder_ij - i;
256 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
257 const int faceBasisOrdinal = numEdges*num1DEdgeFunctions + numVertices + localFaceBasisOrdinal;
258 const auto & faceValue = output_(faceBasisOrdinal,pointOrdinal);
259 const int alphaOrdinal = (i+j)-3;
260 localFaceBasisOrdinal++;
261
262 const int fieldOrdinal = fieldOrdinalOffset + localInteriorBasisOrdinal;
263 const auto & jacobiValue = jacobi_values1_at_point(alphaOrdinal,k);
264 output_(fieldOrdinal,pointOrdinal) = faceValue * jacobiValue;
265 localInteriorBasisOrdinal++;
266 } // end i loop
267 } // end totalPolyOrder_ij loop
268 } // end totalPolyOrder_ijk loop
269 fieldOrdinalOffset += numInteriorBasisFunctions;
270 } // end OPERATOR_VALUE
271 break;
272 case OPERATOR_GRAD:
273 case OPERATOR_D1:
274 {
275 // vertex functions
276 if (defineVertexFunctions_)
277 {
278 // standard, "CG" basis case
279 // first vertex function is 1-x-y-z
280 output_(0,pointOrdinal,0) = -1.0;
281 output_(0,pointOrdinal,1) = -1.0;
282 output_(0,pointOrdinal,2) = -1.0;
283 }
284 else
285 {
286 // "DG" basis case
287 // here, the first "vertex" function is 1, so the derivative is 0:
288 output_(0,pointOrdinal,0) = 0.0;
289 output_(0,pointOrdinal,1) = 0.0;
290 output_(0,pointOrdinal,2) = 0.0;
291 }
292 // second vertex function is x
293 output_(1,pointOrdinal,0) = 1.0;
294 output_(1,pointOrdinal,1) = 0.0;
295 output_(1,pointOrdinal,2) = 0.0;
296 // third vertex function is y
297 output_(2,pointOrdinal,0) = 0.0;
298 output_(2,pointOrdinal,1) = 1.0;
299 output_(2,pointOrdinal,2) = 0.0;
300 // fourth vertex function is z
301 output_(3,pointOrdinal,0) = 0.0;
302 output_(3,pointOrdinal,1) = 0.0;
303 output_(3,pointOrdinal,2) = 1.0;
304
305 // edge functions
306 int fieldOrdinalOffset = numVertices;
307 /*
308 Per Fuentes et al. (see Appendix E.1, E.2), the edge functions, defined for i ≥ 2, are
309 [L_i](s0,s1) = L_i(s1; s0+s1)
310 and have gradients:
311 grad [L_i](s0,s1) = [P_{i-1}](s0,s1) grad s1 + [R_{i-1}](s0,s1) grad (s0 + s1)
312 where
313 [R_{i-1}](s0,s1) = R_{i-1}(s1; s0+s1) = d/dt L_{i}(s0; s0+s1)
314 The P_i we have implemented in shiftedScaledLegendreValues, while d/dt L_{i+1} is
315 implemented in shiftedScaledIntegratedLegendreValues_dt.
316 */
317 // rename the scratch memory to match our usage here:
318 auto & P_i_minus_1 = legendre_values1_at_point;
319 auto & L_i_dt = legendre_values2_at_point;
320 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
321 {
322 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
323 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
324
325 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
326 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
327 const auto & s0_dz = lambda_dz[edge_start_[edgeOrdinal]];
328 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
329 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
330 const auto & s1_dz = lambda_dz[ edge_end_[edgeOrdinal]];
331
332 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
333 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
334 for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
335 {
336 // the first two (integrated) Legendre functions are essentially the vertex functions; hence the +2 here:
337 const int i = edgeFunctionOrdinal+2;
338 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
339 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
340 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,2) = P_i_minus_1(i-1) * s1_dz + L_i_dt(i) * (s1_dz + s0_dz);
341 }
342 fieldOrdinalOffset += num1DEdgeFunctions;
343 }
344
345 /*
346 Fuentes et al give the face functions as phi_{ij}, with gradient:
347 grad phi_{ij}(s0,s1,s2) = [L^{2i}_j](s0+s1,s2) grad [L_i](s0,s1) + [L_i](s0,s1) grad [L^{2i}_j](s0+s1,s2)
348 where:
349 - grad [L_i](s0,s1) is the edge function gradient we computed above
350 - [L_i](s0,s1) is the edge function which we have implemented above (in OPERATOR_VALUE)
351 - L^{2i}_j is a Jacobi polynomial with:
352 [L^{2i}_j](s0,s1) = L^{2i}_j(s1;s0+s1)
353 and the gradient for j ≥ 1 is
354 grad [L^{2i}_j](s0,s1) = [P^{2i}_{j-1}](s0,s1) grad s1 + [R^{2i}_{j-1}(s0,s1)] grad (s0 + s1)
355 Here,
356 [P^{2i}_{j-1}](s0,s1) = P^{2i}_{j-1}(s1,s0+s1)
357 and
358 [R^{2i}_{j-1}(s0,s1)] = d/dt L^{2i}_j(s1,s0+s1)
359 We have implemented P^{alpha}_{j} as shiftedScaledJacobiValues,
360 and d/dt L^{alpha}_{j} as shiftedScaledIntegratedJacobiValues_dt.
361 */
362 // rename the scratch memory to match our usage here:
363 auto & L_i = legendre_values2_at_point;
364 auto & L_2i_j_dt = jacobi_values1_at_point;
365 auto & L_2i_j = jacobi_values2_at_point;
366 auto & P_2i_j_minus_1 = jacobi_values3_at_point;
367
368 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
369 {
370 const auto & s0 = lambda[face_vertex_0[faceOrdinal]];
371 const auto & s1 = lambda[face_vertex_1[faceOrdinal]];
372 const auto & s2 = lambda[face_vertex_2[faceOrdinal]];
373 Polynomials::shiftedScaledIntegratedLegendreValues(L_i, polyOrder_, s1, s0+s1);
374
375 const PointScalar jacobiScaling = s0 + s1 + s2;
376
377 // compute integrated Jacobi values for each desired value of alpha
378 for (int n=2; n<=polyOrder_; n++)
379 {
380 const double alpha = n*2;
381 const int alphaOrdinal = n-2;
382 using Kokkos::subview;
383 using Kokkos::ALL;
384 auto L_2i_j_dt_alpha = subview(L_2i_j_dt, alphaOrdinal, ALL);
385 auto L_2i_j_alpha = subview(L_2i_j, alphaOrdinal, ALL);
386 auto P_2i_j_minus_1_alpha = subview(P_2i_j_minus_1, alphaOrdinal, ALL);
387 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2i_j_dt_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
388 Polynomials::shiftedScaledIntegratedJacobiValues (L_2i_j_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
389 Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1_alpha, alpha, polyOrder_-1, s2, jacobiScaling);
390 }
391
392 const int edgeOrdinal = face_ordinal_of_first_edge[faceOrdinal];
393 int localFaceOrdinal = 0;
394 for (int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
395 {
396 for (int i=2; i<totalPolyOrder; i++)
397 {
398 const int edgeBasisOrdinal = edgeOrdinal*num1DEdgeFunctions + i-2 + numVertices;
399 const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
400 const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
401 const auto & grad_L_i_dz = output_(edgeBasisOrdinal,pointOrdinal,2);
402
403 const int alphaOrdinal = i-2;
404
405 const auto & s0_dx = lambda_dx[face_vertex_0[faceOrdinal]];
406 const auto & s0_dy = lambda_dy[face_vertex_0[faceOrdinal]];
407 const auto & s0_dz = lambda_dz[face_vertex_0[faceOrdinal]];
408 const auto & s1_dx = lambda_dx[face_vertex_1[faceOrdinal]];
409 const auto & s1_dy = lambda_dy[face_vertex_1[faceOrdinal]];
410 const auto & s1_dz = lambda_dz[face_vertex_1[faceOrdinal]];
411 const auto & s2_dx = lambda_dx[face_vertex_2[faceOrdinal]];
412 const auto & s2_dy = lambda_dy[face_vertex_2[faceOrdinal]];
413 const auto & s2_dz = lambda_dz[face_vertex_2[faceOrdinal]];
414
415 int j = totalPolyOrder - i;
416
417 // put references to the entries of interest in like-named variables with lowercase first letters
418 auto & l_2i_j = L_2i_j(alphaOrdinal,j);
419 auto & l_i = L_i(i);
420 auto & l_2i_j_dt = L_2i_j_dt(alphaOrdinal,j);
421 auto & p_2i_j_minus_1 = P_2i_j_minus_1(alphaOrdinal,j-1);
422
423 const OutputScalar basisValue_dx = l_2i_j * grad_L_i_dx + l_i * (p_2i_j_minus_1 * s2_dx + l_2i_j_dt * (s0_dx + s1_dx + s2_dx));
424 const OutputScalar basisValue_dy = l_2i_j * grad_L_i_dy + l_i * (p_2i_j_minus_1 * s2_dy + l_2i_j_dt * (s0_dy + s1_dy + s2_dy));
425 const OutputScalar basisValue_dz = l_2i_j * grad_L_i_dz + l_i * (p_2i_j_minus_1 * s2_dz + l_2i_j_dt * (s0_dz + s1_dz + s2_dz));
426
427 const int fieldOrdinal = fieldOrdinalOffset + localFaceOrdinal;
428
429 output_(fieldOrdinal,pointOrdinal,0) = basisValue_dx;
430 output_(fieldOrdinal,pointOrdinal,1) = basisValue_dy;
431 output_(fieldOrdinal,pointOrdinal,2) = basisValue_dz;
432
433 localFaceOrdinal++;
434 }
435 }
436 fieldOrdinalOffset += numFaceBasisFunctionsPerFace;
437 }
438 // interior functions
439 /*
440 Per Fuentes et al. (see Appendix E.1, E.2), the interior functions, defined for i ≥ 2, are
441 phi_ij(lambda_012) [L^{2(i+j)}_k](1-lambda_3,lambda_3) = phi_ij(lambda_012) L^{2(i+j)}_k (lambda_3; 1)
442 and have gradients:
443 L^{2(i+j)}_k (lambda_3; 1) grad (phi_ij(lambda_012)) + phi_ij(lambda_012) grad (L^{2(i+j)}_k (lambda_3; 1))
444 where:
445 - phi_ij(lambda_012) is the (i,j) basis function on face 012,
446 - L^{alpha}_j(t0; t1) is the jth integrated Jacobi polynomial
447 and the gradient of the integrated Jacobi polynomial can be computed:
448 - grad L^{alpha}_j(t0; t1) = P^{alpha}_{j-1} (t0;t1) grad t0 + R^{alpha}_{j-1}(t0,t1) grad t1
449 Here, t1=1, so this simplifies to:
450 - grad L^{alpha}_j(t0; t1) = P^{alpha}_{j-1} (t0;t1) grad t0
451
452 The P_i we have implemented in shiftedScaledJacobiValues.
453 */
454 // rename the scratch memory to match our usage here:
455 auto & L_alpha = jacobi_values1_at_point;
456 auto & P_alpha = jacobi_values2_at_point;
457
458 // precompute values used in face ordinal 0:
459 {
460 const auto & s0 = lambda[0];
461 const auto & s1 = lambda[1];
462 const auto & s2 = lambda[2];
463 // Legendre:
464 Polynomials::shiftedScaledIntegratedLegendreValues(legendre_values1_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
465
466 // Jacobi for each desired alpha value:
467 const PointScalar jacobiScaling = s0 + s1 + s2;
468 for (int n=2; n<=polyOrder_; n++)
469 {
470 const double alpha = n*2;
471 const int alphaOrdinal = n-2;
472 using Kokkos::subview;
473 using Kokkos::ALL;
474 auto jacobi_alpha = subview(jacobi_values3_at_point, alphaOrdinal, ALL);
475 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_alpha, alpha, polyOrder_-2, s2, jacobiScaling);
476 }
477 }
478
479 // interior
480 for (int n=3; n<=polyOrder_; n++)
481 {
482 const double alpha = n*2;
483 const double jacobiScaling = 1.0;
484 const int alphaOrdinal = n-3;
485 using Kokkos::subview;
486 using Kokkos::ALL;
487
488 // values for interior functions:
489 auto L = subview(L_alpha, alphaOrdinal, ALL);
490 auto P = subview(P_alpha, alphaOrdinal, ALL);
491 Polynomials::shiftedScaledIntegratedJacobiValues(L, alpha, polyOrder_-3, lambda[3], jacobiScaling);
492 Polynomials::shiftedScaledJacobiValues (P, alpha, polyOrder_-3, lambda[3], jacobiScaling);
493 }
494
495 const int min_i = 2;
496 const int min_j = 1;
497 const int min_k = 1;
498 const int min_ij = min_i + min_j;
499 const int min_ijk = min_ij + min_k;
500 int localInteriorBasisOrdinal = 0;
501 for (int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= polyOrder_; totalPolyOrder_ijk++)
502 {
503 int localFaceBasisOrdinal = 0;
504 for (int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
505 {
506 for (int i=2; i <= totalPolyOrder_ij-min_j; i++)
507 {
508 const int j = totalPolyOrder_ij - i;
509 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
510 // interior functions use basis values belonging to the first face, 012
511 const int faceBasisOrdinal = numEdges*num1DEdgeFunctions + numVertices + localFaceBasisOrdinal;
512
513 const auto & faceValue_dx = output_(faceBasisOrdinal,pointOrdinal,0);
514 const auto & faceValue_dy = output_(faceBasisOrdinal,pointOrdinal,1);
515 const auto & faceValue_dz = output_(faceBasisOrdinal,pointOrdinal,2);
516
517 // determine faceValue (on face 0)
518 OutputScalar faceValue;
519 {
520 const auto & edgeValue = legendre_values1_at_point(i);
521 const int alphaOrdinal = i-2;
522 const auto & jacobiValue = jacobi_values3_at_point(alphaOrdinal,j);
523 faceValue = edgeValue * jacobiValue;
524 }
525 localFaceBasisOrdinal++;
526
527 const int alphaOrdinal = (i+j)-3;
528
529 const int fieldOrdinal = fieldOrdinalOffset + localInteriorBasisOrdinal;
530 const auto & integratedJacobiValue = L_alpha(alphaOrdinal,k);
531 const auto & jacobiValue = P_alpha(alphaOrdinal,k-1);
532 output_(fieldOrdinal,pointOrdinal,0) = integratedJacobiValue * faceValue_dx + faceValue * jacobiValue * lambda_dx[3];
533 output_(fieldOrdinal,pointOrdinal,1) = integratedJacobiValue * faceValue_dy + faceValue * jacobiValue * lambda_dy[3];
534 output_(fieldOrdinal,pointOrdinal,2) = integratedJacobiValue * faceValue_dz + faceValue * jacobiValue * lambda_dz[3];
535
536 localInteriorBasisOrdinal++;
537 } // end i loop
538 } // end totalPolyOrder_ij loop
539 } // end totalPolyOrder_ijk loop
540 fieldOrdinalOffset += numInteriorBasisFunctions;
541 }
542 break;
543 case OPERATOR_D2:
544 case OPERATOR_D3:
545 case OPERATOR_D4:
546 case OPERATOR_D5:
547 case OPERATOR_D6:
548 case OPERATOR_D7:
549 case OPERATOR_D8:
550 case OPERATOR_D9:
551 case OPERATOR_D10:
552 INTREPID2_TEST_FOR_ABORT( true,
553 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
554 default:
555 // unsupported operator type
556 device_assert(false);
557 }
558 }
559
560 // Provide the shared memory capacity.
561 // This function takes the team_size as an argument,
562 // which allows team_size-dependent allocations.
563 size_t team_shmem_size (int team_size) const
564 {
565 // we will use shared memory to create a fast buffer for basis computations
566 // for the (integrated) Legendre computations, we just need p+1 values stored
567 // for the (integrated) Jacobi computations, though, we want (p+1)*(# alpha values)
568 // alpha is either 2i or 2(i+j), where i=2,…,p or i+j=3,…,p. So there are at most (p-1) alpha values needed.
569 // We can have up to 3 of the (integrated) Jacobi values needed at once.
570 const int numAlphaValues = std::max(polyOrder_-1, 1); // make it at least 1 so we can avoid zero-extent ranks…
571 size_t shmem_size = 0;
572 if (fad_size_output_ > 0)
573 {
574 // Legendre:
575 shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
576 // Jacobi:
577 shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1, fad_size_output_);
578 }
579 else
580 {
581 // Legendre:
582 shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1);
583 // Jacobi:
584 shmem_size += 3 * OutputScratchView2D::shmem_size(numAlphaValues, polyOrder_ + 1);
585 }
586
587 return shmem_size;
588 }
589 };
590
608 template<typename DeviceType,
609 typename OutputScalar = double,
610 typename PointScalar = double,
611 bool defineVertexFunctions = true> // if defineVertexFunctions is true, first four basis functions are 1-x-y-z, x, y, and z. Otherwise, they are 1, x, y, and z.
613 : public Basis<DeviceType,OutputScalar,PointScalar>
614 {
615 public:
617
620
621 using typename BasisBase::OutputViewType;
622 using typename BasisBase::PointViewType;
623 using typename BasisBase::ScalarViewType;
624
625 using typename BasisBase::ExecutionSpace;
626
627 protected:
628 int polyOrder_; // the maximum order of the polynomial
629 EPointType pointType_;
630 public:
641 IntegratedLegendreBasis_HGRAD_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
642 :
643 polyOrder_(polyOrder),
644 pointType_(pointType)
645 {
646 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
647 this->basisCardinality_ = ((polyOrder+1) * (polyOrder+2) * (polyOrder+3)) / 6;
648 this->basisDegree_ = polyOrder;
649 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
650 this->basisType_ = BASIS_FEM_HIERARCHICAL;
651 this->basisCoordinates_ = COORDINATES_CARTESIAN;
652 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
653
654 const int degreeLength = 1;
655 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) tetrahedron polynomial degree lookup", this->basisCardinality_, degreeLength);
656 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) tetrahedron polynomial H^1 degree lookup", this->basisCardinality_, degreeLength);
657
658 int fieldOrdinalOffset = 0;
659 // **** vertex functions **** //
660 const int numVertices = this->basisCellTopology_.getVertexCount();
661 const int numFunctionsPerVertex = 1;
662 const int numVertexFunctions = numVertices * numFunctionsPerVertex;
663 for (int i=0; i<numVertexFunctions; i++)
664 {
665 // for H(grad) on tetrahedron, if defineVertexFunctions is false, first four basis members are linear
666 // if not, then the only difference is that the first member is constant
667 this->fieldOrdinalPolynomialDegree_ (i,0) = 1;
668 this->fieldOrdinalH1PolynomialDegree_(i,0) = 1;
669 }
670 if (!defineVertexFunctions)
671 {
672 this->fieldOrdinalPolynomialDegree_ (0,0) = 0;
673 this->fieldOrdinalH1PolynomialDegree_(0,0) = 0;
674 }
675 fieldOrdinalOffset += numVertexFunctions;
676
677 // **** edge functions **** //
678 const int numFunctionsPerEdge = polyOrder - 1; // bubble functions: all but the vertices
679 const int numEdges = this->basisCellTopology_.getEdgeCount();
680 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
681 {
682 for (int i=0; i<numFunctionsPerEdge; i++)
683 {
684 this->fieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2
685 this->fieldOrdinalH1PolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2
686 }
687 fieldOrdinalOffset += numFunctionsPerEdge;
688 }
689
690 // **** face functions **** //
691 const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
692 const int numFaces = 4;
693 for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
694 {
695 for (int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
696 {
697 const int totalFaceDofs = (totalPolyOrder-2)*(totalPolyOrder-1)/2;
698 const int totalFaceDofsPrevious = (totalPolyOrder-3)*(totalPolyOrder-2)/2;
699 const int faceDofsForPolyOrder = totalFaceDofs - totalFaceDofsPrevious;
700 for (int i=0; i<faceDofsForPolyOrder; i++)
701 {
702 this->fieldOrdinalPolynomialDegree_ (fieldOrdinalOffset,0) = totalPolyOrder;
703 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinalOffset,0) = totalPolyOrder;
704 fieldOrdinalOffset++;
705 }
706 }
707 }
708
709 // **** interior functions **** //
710 const int numFunctionsPerVolume = ((polyOrder-1)*(polyOrder-2)*(polyOrder-3))/6;
711 const int numVolumes = 1; // interior
712 for (int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
713 {
714 for (int totalPolyOrder=4; totalPolyOrder<=polyOrder_; totalPolyOrder++)
715 {
716 const int totalInteriorDofs = (totalPolyOrder-3)*(totalPolyOrder-2)*(totalPolyOrder-1)/6;
717 const int totalInteriorDofsPrevious = (totalPolyOrder-4)*(totalPolyOrder-3)*(totalPolyOrder-2)/6;
718 const int interiorDofsForPolyOrder = totalInteriorDofs - totalInteriorDofsPrevious;
719
720 for (int i=0; i<interiorDofsForPolyOrder; i++)
721 {
722 this->fieldOrdinalPolynomialDegree_ (fieldOrdinalOffset,0) = totalPolyOrder;
723 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinalOffset,0) = totalPolyOrder;
724 fieldOrdinalOffset++;
725 }
726 }
727 }
728
729 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
730
731 // initialize tags
732 {
733 // ESEAS numbers tetrahedron faces differently from Intrepid2
734 // ESEAS: 012, 013, 123, 023
735 // Intrepid2: 013, 123, 032, 021
736 const int intrepid2FaceOrdinals[4] {3,0,1,2}; // index is the ESEAS face ordinal; value is the intrepid2 ordinal
737
738 const auto & cardinality = this->basisCardinality_;
739
740 // Basis-dependent initializations
741 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
742 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
743 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
744 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
745
746 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
747 const int vertexDim = 0, edgeDim = 1, faceDim = 2, volumeDim = 3;
748
749 if (defineVertexFunctions) {
750 {
751 int tagNumber = 0;
752 for (int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
753 {
754 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerVertex; functionOrdinal++)
755 {
756 tagView(tagNumber*tagSize+0) = vertexDim; // vertex dimension
757 tagView(tagNumber*tagSize+1) = vertexOrdinal; // vertex id
758 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
759 tagView(tagNumber*tagSize+3) = numFunctionsPerVertex; // total number of dofs in this vertex
760 tagNumber++;
761 }
762 }
763 for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
764 {
765 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
766 {
767 tagView(tagNumber*tagSize+0) = edgeDim; // edge dimension
768 tagView(tagNumber*tagSize+1) = edgeOrdinal; // edge id
769 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
770 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge; // total number of dofs on this edge
771 tagNumber++;
772 }
773 }
774 for (int faceOrdinalESEAS=0; faceOrdinalESEAS<numFaces; faceOrdinalESEAS++)
775 {
776 int faceOrdinalIntrepid2 = intrepid2FaceOrdinals[faceOrdinalESEAS];
777 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
778 {
779 tagView(tagNumber*tagSize+0) = faceDim; // face dimension
780 tagView(tagNumber*tagSize+1) = faceOrdinalIntrepid2; // face id
781 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
782 tagView(tagNumber*tagSize+3) = numFunctionsPerFace; // total number of dofs on this face
783 tagNumber++;
784 }
785 }
786 for (int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
787 {
788 for (int functionOrdinal=0; functionOrdinal<numFunctionsPerVolume; functionOrdinal++)
789 {
790 tagView(tagNumber*tagSize+0) = volumeDim; // volume dimension
791 tagView(tagNumber*tagSize+1) = volumeOrdinal; // volume id
792 tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
793 tagView(tagNumber*tagSize+3) = numFunctionsPerVolume; // total number of dofs in this volume
794 tagNumber++;
795 }
796 }
797 INTREPID2_TEST_FOR_EXCEPTION(tagNumber != this->basisCardinality_, std::invalid_argument, "Internal error: basis tag enumeration is incorrect");
798 }
799 } else {
800 for (ordinal_type i=0;i<cardinality;++i) {
801 tagView(i*tagSize+0) = volumeDim; // volume dimension
802 tagView(i*tagSize+1) = 0; // volume ordinal
803 tagView(i*tagSize+2) = i; // local dof id
804 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
805 }
806 }
807
808 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
809 // tags are constructed on host
811 this->ordinalToTag_,
812 tagView,
813 this->basisCardinality_,
814 tagSize,
815 posScDim,
816 posScOrd,
817 posDfOrd);
818 }
819 }
820
825 const char* getName() const override {
826 return "Intrepid2_IntegratedLegendreBasis_HGRAD_TET";
827 }
828
831 virtual bool requireOrientation() const override {
832 return (this->getDegree() > 2);
833 }
834
835 // since the getValues() below only overrides the FEM variant, we specify that
836 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
837 // (It's an error to use the FVD variant on this basis.)
839
858 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
859 const EOperator operatorType = OPERATOR_VALUE ) const override
860 {
861 auto numPoints = inputPoints.extent_int(0);
862
864
865 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
866
867 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
868 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
869 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
870 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
871
872 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
873 Kokkos::parallel_for("Hierarchical_HGRAD_TET_Functor", policy, functor);
874 }
875
885 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
886 if(subCellDim == 1) {
887 return Teuchos::rcp(new
889 (this->basisDegree_));
890 } else if(subCellDim == 2) {
891 return Teuchos::rcp(new
893 (this->basisDegree_));
894 }
895 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
896 }
897
903 getHostBasis() const override {
904 using HostDeviceType = typename Kokkos::HostSpace::device_type;
906 return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
907 }
908 };
909} // end namespace Intrepid2
910
911#endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_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(grad) basis on the line based on integrated Legendre polynomials.
H(grad) 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
OrdinalTypeArray2DHost fieldOrdinalH1PolynomialDegree_
H^1 polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
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.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line.
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 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 void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
virtual bool requireOrientation() const override
True if orientation is required.
IntegratedLegendreBasis_HGRAD_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line: e...
Functor for computing values for the IntegratedLegendreBasis_HGRAD_TET class.