Intrepid2
Intrepid2_IntegratedLegendreBasis_HGRAD_TRI.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_TRI_h
50 #define Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_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 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  int polyOrder_;
85  bool defineVertexFunctions_;
86  int numFields_, numPoints_;
87 
88  size_t fad_size_output_;
89 
90  static const int numVertices = 3;
91  static const int numEdges = 3;
92  const int edge_start_[numEdges] = {0,1,0}; // edge i is from edge_start_[i] to edge_end_[i]
93  const int edge_end_[numEdges] = {1,2,2}; // edge i is from edge_start_[i] to edge_end_[i]
94 
95  Hierarchical_HGRAD_TRI_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints,
96  int polyOrder, bool defineVertexFunctions)
97  : opType_(opType), output_(output), inputPoints_(inputPoints),
98  polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
99  fad_size_output_(getScalarDimensionForView(output))
100  {
101  numFields_ = output.extent_int(0);
102  numPoints_ = output.extent_int(1);
103  INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
104  INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument, "output field size does not match basis cardinality");
105  }
106 
107  KOKKOS_INLINE_FUNCTION
108  void operator()( const TeamMember & teamMember ) const
109  {
110  auto pointOrdinal = teamMember.league_rank();
111  OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
112  if (fad_size_output_ > 0) {
113  edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
114  jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
115  other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
116  other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
117  }
118  else {
119  edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
120  jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
121  other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
122  other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
123  }
124 
125  const auto & x = inputPoints_(pointOrdinal,0);
126  const auto & y = inputPoints_(pointOrdinal,1);
127 
128  // write as barycentric coordinates:
129  const PointScalar lambda[3] = {1. - x - y, x, y};
130  const PointScalar lambda_dx[3] = {-1., 1., 0.};
131  const PointScalar lambda_dy[3] = {-1., 0., 1.};
132 
133  const int num1DEdgeFunctions = polyOrder_ - 1;
134 
135  switch (opType_)
136  {
137  case OPERATOR_VALUE:
138  {
139  // vertex functions come first, according to vertex ordering: (0,0), (1,0), (0,1)
140  for (int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
141  {
142  output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
143  }
144  if (!defineVertexFunctions_)
145  {
146  // "DG" basis case
147  // here, we overwrite the first vertex function with 1:
148  output_(0,pointOrdinal) = 1.0;
149  }
150 
151  // edge functions
152  int fieldOrdinalOffset = 3;
153  for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
154  {
155  const auto & s0 = lambda[edge_start_[edgeOrdinal]];
156  const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
157 
158  Polynomials::shiftedScaledIntegratedLegendreValues(edge_field_values_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
159  for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
160  {
161  // the first two integrated legendre functions are essentially the vertex functions; hence the +2 on on the RHS here:
162  output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = edge_field_values_at_point(edgeFunctionOrdinal+2);
163  }
164  fieldOrdinalOffset += num1DEdgeFunctions;
165  }
166 
167  // face functions
168  {
169  // these functions multiply the edge functions from the 01 edge by integrated Jacobi functions, appropriately scaled
170  const double jacobiScaling = 1.0; // s0 + s1 + s2
171 
172  for (int i=2; i<polyOrder_; i++)
173  {
174  const int edgeBasisOrdinal = i+numVertices-2; // i+1: where the value of the edge function is stored in output_
175  const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
176  const double alpha = i*2.0;
177 
178  Polynomials::integratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-2, lambda[2], jacobiScaling);
179  for (int j=1; i+j <= polyOrder_; j++)
180  {
181  const auto & jacobiValue = jacobi_values_at_point(j);
182  output_(fieldOrdinalOffset,pointOrdinal) = edgeValue * jacobiValue;
183  fieldOrdinalOffset++;
184  }
185  }
186  }
187  } // end OPERATOR_VALUE
188  break;
189  case OPERATOR_GRAD:
190  case OPERATOR_D1:
191  {
192  // vertex functions
193  if (defineVertexFunctions_)
194  {
195  // standard, "CG" basis case
196  // first vertex function is 1-x-y
197  output_(0,pointOrdinal,0) = -1.0;
198  output_(0,pointOrdinal,1) = -1.0;
199  }
200  else
201  {
202  // "DG" basis case
203  // here, the first "vertex" function is 1, so the derivative is 0:
204  output_(0,pointOrdinal,0) = 0.0;
205  output_(0,pointOrdinal,1) = 0.0;
206  }
207  // second vertex function is x
208  output_(1,pointOrdinal,0) = 1.0;
209  output_(1,pointOrdinal,1) = 0.0;
210  // third vertex function is y
211  output_(2,pointOrdinal,0) = 0.0;
212  output_(2,pointOrdinal,1) = 1.0;
213 
214  // edge functions
215  int fieldOrdinalOffset = 3;
216  /*
217  Per Fuentes et al. (see Appendix E.1, E.2), the edge functions, defined for i ≥ 2, are
218  [L_i](s0,s1) = L_i(s1; s0+s1)
219  and have gradients:
220  grad [L_i](s0,s1) = [P_{i-1}](s0,s1) grad s1 + [R_{i-1}](s0,s1) grad (s0 + s1)
221  where
222  [R_{i-1}](s0,s1) = R_{i-1}(s1; s0+s1) = d/dt L_{i}(s0; s0+s1)
223  The P_i we have implemented in shiftedScaledLegendreValues, while d/dt L_{i+1} is
224  implemented in shiftedScaledIntegratedLegendreValues_dt.
225  */
226  // rename the scratch memory to match our usage here:
227  auto & P_i_minus_1 = edge_field_values_at_point;
228  auto & L_i_dt = jacobi_values_at_point;
229  for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
230  {
231  const auto & s0 = lambda[edge_start_[edgeOrdinal]];
232  const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
233 
234  const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
235  const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
236  const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
237  const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
238 
239  Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
240  Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
241  for (int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
242  {
243  // the first two (integrated) Legendre functions are essentially the vertex functions; hence the +2 here:
244  const int i = edgeFunctionOrdinal+2;
245  output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
246  output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
247  }
248  fieldOrdinalOffset += num1DEdgeFunctions;
249  }
250 
251  /*
252  Fuentes et al give the face functions as phi_{ij}, with gradient:
253  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)
254  where:
255  - grad [L_i](s0,s1) is the edge function gradient we computed above
256  - [L_i](s0,s1) is the edge function which we have implemented above (in OPERATOR_VALUE)
257  - L^{2i}_j is a Jacobi polynomial with:
258  [L^{2i}_j](s0,s1) = L^{2i}_j(s1;s0+s1)
259  and the gradient for j ≥ 1 is
260  grad [L^{2i}_j](s0,s1) = [P^{2i}_{j-1}](s0,s1) grad s1 + [R^{2i}_{j-1}(s0,s1)] grad (s0 + s1)
261  Here,
262  [P^{2i}_{j-1}](s0,s1) = P^{2i}_{j-1}(s1,s0+s1)
263  and
264  [R^{2i}_{j-1}(s0,s1)] = d/dt L^{2i}_j(s1,s0+s1)
265  We have implemented P^{alpha}_{j} as shiftedScaledJacobiValues,
266  and d/dt L^{alpha}_{j} as integratedJacobiValues_dt.
267  */
268  // rename the scratch memory to match our usage here:
269  auto & P_2i_j_minus_1 = edge_field_values_at_point;
270  auto & L_2i_j_dt = jacobi_values_at_point;
271  auto & L_i = other_values_at_point;
272  auto & L_2i_j = other_values2_at_point;
273  {
274  // face functions multiply the edge functions from the 01 edge by integrated Jacobi functions, appropriately scaled
275  const double jacobiScaling = 1.0; // s0 + s1 + s2
276 
277  for (int i=2; i<polyOrder_; i++)
278  {
279  // the edge function here is for edge 01, in the first set of edge functions.
280  const int edgeBasisOrdinal = i+numVertices-2; // i+1: where the value of the edge function is stored in output_
281  const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
282  const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
283 
284  const double alpha = i*2.0;
285 
286  Polynomials::shiftedScaledIntegratedLegendreValues(L_i, polyOrder_, lambda[1], lambda[0]+lambda[1]);
287  Polynomials::integratedJacobiValues_dt( L_2i_j_dt, alpha, polyOrder_, lambda[2], jacobiScaling);
288  Polynomials::integratedJacobiValues ( L_2i_j, alpha, polyOrder_, lambda[2], jacobiScaling);
289  Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1, alpha, polyOrder_-1, lambda[2], jacobiScaling);
290 
291  const auto & s0_dx = lambda_dx[0];
292  const auto & s0_dy = lambda_dy[0];
293  const auto & s1_dx = lambda_dx[1];
294  const auto & s1_dy = lambda_dy[1];
295  const auto & s2_dx = lambda_dx[2];
296  const auto & s2_dy = lambda_dy[2];
297 
298  for (int j=1; i+j <= polyOrder_; j++)
299  {
300  const OutputScalar basisValue_dx = L_2i_j(j) * grad_L_i_dx + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dx + L_2i_j_dt(j) * (s0_dx + s1_dx + s2_dx));
301  const OutputScalar basisValue_dy = L_2i_j(j) * grad_L_i_dy + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dy + L_2i_j_dt(j) * (s0_dy + s1_dy + s2_dy));
302 
303  output_(fieldOrdinalOffset,pointOrdinal,0) = basisValue_dx;
304  output_(fieldOrdinalOffset,pointOrdinal,1) = basisValue_dy;
305  fieldOrdinalOffset++;
306  }
307  }
308  }
309  }
310  break;
311  case OPERATOR_D2:
312  case OPERATOR_D3:
313  case OPERATOR_D4:
314  case OPERATOR_D5:
315  case OPERATOR_D6:
316  case OPERATOR_D7:
317  case OPERATOR_D8:
318  case OPERATOR_D9:
319  case OPERATOR_D10:
320  INTREPID2_TEST_FOR_ABORT( true,
321  ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
322  default:
323  // unsupported operator type
324  device_assert(false);
325  }
326  }
327 
328  // Provide the shared memory capacity.
329  // This function takes the team_size as an argument,
330  // which allows team_size-dependent allocations.
331  size_t team_shmem_size (int team_size) const
332  {
333  // we will use shared memory to create a fast buffer for basis computations
334  size_t shmem_size = 0;
335  if (fad_size_output_ > 0)
336  shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
337  else
338  shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
339 
340  return shmem_size;
341  }
342  };
343 
361  template<typename DeviceType,
362  typename OutputScalar = double,
363  typename PointScalar = double,
364  bool defineVertexFunctions = true> // if defineVertexFunctions is true, first three basis functions are 1-x-y, x, and y. Otherwise, they are 1, x, and y.
366  : public Basis<DeviceType,OutputScalar,PointScalar>
367  {
368  public:
370 
371  using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost;
372  using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost;
373 
374  using OutputViewType = typename BasisBase::OutputViewType;
375  using PointViewType = typename BasisBase::PointViewType ;
376  using ScalarViewType = typename BasisBase::ScalarViewType;
377  protected:
378  int polyOrder_; // the maximum order of the polynomial
379  EPointType pointType_;
380  public:
391  IntegratedLegendreBasis_HGRAD_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
392  :
393  polyOrder_(polyOrder),
394  pointType_(pointType)
395  {
396  INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
397 
398  this->basisCardinality_ = ((polyOrder+2) * (polyOrder+1)) / 2;
399  this->basisDegree_ = polyOrder;
400  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
401  this->basisType_ = BASIS_FEM_HIERARCHICAL;
402  this->basisCoordinates_ = COORDINATES_CARTESIAN;
403  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
404 
405  const int degreeLength = 1;
406  this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
407 
408  int fieldOrdinalOffset = 0;
409  // **** vertex functions **** //
410  const int numVertices = this->basisCellTopology_.getVertexCount();
411  const int numFunctionsPerVertex = 1;
412  const int numVertexFunctions = numVertices * numFunctionsPerVertex;
413  for (int i=0; i<numVertexFunctions; i++)
414  {
415  // for H(grad) on triangle, if defineVertexFunctions is false, first three basis members are linear
416  // if not, then the only difference is that the first member is constant
417  this->fieldOrdinalPolynomialDegree_(i,0) = 1;
418  }
419  if (!defineVertexFunctions)
420  {
421  this->fieldOrdinalPolynomialDegree_(0,0) = 0;
422  }
423  fieldOrdinalOffset += numVertexFunctions;
424 
425  // **** edge functions **** //
426  const int numFunctionsPerEdge = polyOrder - 1; // bubble functions: all but the vertices
427  const int numEdges = this->basisCellTopology_.getEdgeCount();
428  for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
429  {
430  for (int i=0; i<numFunctionsPerEdge; i++)
431  {
432  this->fieldOrdinalPolynomialDegree_(i+fieldOrdinalOffset,0) = i+2; // vertex functions are 1st order; edge functions start at order 2
433  }
434  fieldOrdinalOffset += numFunctionsPerEdge;
435  }
436 
437  // **** face functions **** //
438  const int max_ij_sum = polyOrder;
439  for (int i=2; i<max_ij_sum; i++)
440  {
441  for (int j=1; i+j<=max_ij_sum; j++)
442  {
443  this->fieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = i+j;
444  fieldOrdinalOffset++;
445  }
446  }
447  const int numFaces = 1;
448  const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
449  INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
450 
451  // initialize tags
452  {
453  const auto & cardinality = this->basisCardinality_;
454 
455  // Basis-dependent initializations
456  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
457  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
458  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
459  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
460 
461  OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
462  const int vertexDim = 0, edgeDim = 1, faceDim = 2;
463 
464  if (defineVertexFunctions) {
465  {
466  int tagNumber = 0;
467  for (int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
468  {
469  for (int functionOrdinal=0; functionOrdinal<numFunctionsPerVertex; functionOrdinal++)
470  {
471  tagView(tagNumber*tagSize+0) = vertexDim; // vertex dimension
472  tagView(tagNumber*tagSize+1) = vertexOrdinal; // vertex id
473  tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
474  tagView(tagNumber*tagSize+3) = numFunctionsPerVertex; // total number of dofs in this vertex
475  tagNumber++;
476  }
477  }
478  for (int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
479  {
480  for (int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
481  {
482  tagView(tagNumber*tagSize+0) = edgeDim; // edge dimension
483  tagView(tagNumber*tagSize+1) = edgeOrdinal; // edge id
484  tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
485  tagView(tagNumber*tagSize+3) = numFunctionsPerEdge; // total number of dofs on this edge
486  tagNumber++;
487  }
488  }
489  for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
490  {
491  for (int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
492  {
493  tagView(tagNumber*tagSize+0) = faceDim; // face dimension
494  tagView(tagNumber*tagSize+1) = faceOrdinal; // face id
495  tagView(tagNumber*tagSize+2) = functionOrdinal; // local dof id
496  tagView(tagNumber*tagSize+3) = numFunctionsPerFace; // total number of dofs on this face
497  tagNumber++;
498  }
499  }
500  }
501  } else {
502  for (ordinal_type i=0;i<cardinality;++i) {
503  tagView(i*tagSize+0) = faceDim; // face dimension
504  tagView(i*tagSize+1) = 0; // face id
505  tagView(i*tagSize+2) = i; // local dof id
506  tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
507  }
508  }
509 
510  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
511  // tags are constructed on host
512  this->setOrdinalTagData(this->tagToOrdinal_,
513  this->ordinalToTag_,
514  tagView,
515  this->basisCardinality_,
516  tagSize,
517  posScDim,
518  posScOrd,
519  posDfOrd);
520  }
521  }
522 
527  const char* getName() const override {
528  return "Intrepid2_IntegratedLegendreBasis_HGRAD_TRI";
529  }
530 
533  virtual bool requireOrientation() const override {
534  return (this->getDegree() > 2);
535  }
536 
537  // since the getValues() below only overrides the FEM variant, we specify that
538  // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
539  // (It's an error to use the FVD variant on this basis.)
540  using BasisBase::getValues;
541 
560  virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
561  const EOperator operatorType = OPERATOR_VALUE ) const override
562  {
563  auto numPoints = inputPoints.extent_int(0);
564 
566 
567  FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
568 
569  const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
570  const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
571  const int vectorSize = std::max(outputVectorSize,pointVectorSize);
572  const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
573 
574  using ExecutionSpace = typename BasisBase::ExecutionSpace;
575 
576  auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
577  Kokkos::parallel_for( policy , functor, "Hierarchical_HGRAD_TRI_Functor");
578  }
579 
589  getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
590  if(subCellDim == 1) {
591  return Teuchos::rcp(new
593  (this->basisDegree_));
594  }
595  INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
596  }
597 
603  getHostBasis() const override {
604  using HostDeviceType = typename Kokkos::HostSpace::device_type;
606  return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
607  }
608  };
609 } // end namespace Intrepid2
610 
611 #endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_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
virtual bool requireOrientation() const override
True if orientation is required.
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.
IntegratedLegendreBasis_HGRAD_TRI(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 (...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
EFunctionSpace functionSpace_
The function space in which the basis is defined.
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.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
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...
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...
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.
Functor for computing values for the IntegratedLegendreBasis_HGRAD_TRI class.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
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...