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