Intrepid2
Intrepid2_IntegratedLegendreBasis_HGRAD_LINE.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_LINE_h
50 #define Intrepid2_IntegratedLegendreBasis_HGRAD_LINE_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_; // OPERATOR_VALUE or OPERATOR_GRAD
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  Hierarchical_HGRAD_LINE_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints,
91  int polyOrder, bool defineVertexFunctions)
92  : opType_(opType), output_(output), inputPoints_(inputPoints),
93  polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
94  fad_size_output_(getScalarDimensionForView(output))
95  {
96  numFields_ = output.extent_int(0);
97  numPoints_ = output.extent_int(1);
98  INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
99  INTREPID2_TEST_FOR_EXCEPTION(numFields_ != polyOrder_+1, std::invalid_argument, "output field size does not match basis cardinality");
100  }
101 
102  KOKKOS_INLINE_FUNCTION
103  void operator()( const TeamMember & teamMember ) const
104  {
105  auto pointOrdinal = teamMember.league_rank();
106  OutputScratchView field_values_at_point;
107  if (fad_size_output_ > 0) {
108  field_values_at_point = OutputScratchView(teamMember.team_shmem(), numFields_, fad_size_output_);
109  }
110  else {
111  field_values_at_point = OutputScratchView(teamMember.team_shmem(), numFields_);
112  }
113 
114  const auto & input_x = inputPoints_(pointOrdinal,0);
115  const bool taking_derivative = (opType_ != OPERATOR_VALUE);
116  const bool callingShiftedScaledLegendre = (opType_ == OPERATOR_VALUE) || (opType_ == OPERATOR_GRAD) || (opType_ == OPERATOR_D1);
117 
118  // shiftedScaledIntegratedLegendreValues{_dx} expects x in [0,1]
119  const PointScalar x = callingShiftedScaledLegendre ? PointScalar((input_x + 1.0)/2.0) : PointScalar(input_x);
120  const double legendreScaling = 1.0;
121  const double outputScaling = taking_derivative ? 0.5 : 1.0; // output scaling -- 0.5 if we take derivatives, 1.0 otherwise
122 
123  switch (opType_)
124  {
125  case OPERATOR_VALUE:
126  // field values are integrated Legendre polynomials, except for the first and second field,
127  // which may be 1 and x or x and 1-x, depending on whether the vertex compatibility flag is set.
128  Polynomials::shiftedScaledIntegratedLegendreValues(field_values_at_point, polyOrder_, x, legendreScaling);
129 
130  // note that because shiftedScaledIntegratedLegendreValues determines field values recursively, there is not much
131  // opportunity at that level for further parallelism
132 
133  if (defineVertexFunctions_)
134  {
135  field_values_at_point(0) = 1. - x;
136  field_values_at_point(1) = x;
137  }
138  break;
139  case OPERATOR_GRAD:
140  case OPERATOR_D1:
141  // field values are Legendre polynomials, except for the first and second field,
142  // which may be 0 and 1 or -1 and 1, depending on whether the vertex compatibility flag is set.
143  Polynomials::shiftedScaledIntegratedLegendreValues_dx(field_values_at_point, polyOrder_, x, legendreScaling);
144 
145  // note that because shiftedScaledIntegratedLegendreValues_dx determines field values recursively, there is not much
146  // opportunity at that level for further parallelism
147 
148  if (defineVertexFunctions_)
149  {
150  field_values_at_point(0) = -1.0; // derivative of 1-x
151  field_values_at_point(1) = 1.0; // derivative of x
152  }
153  break;
154  case OPERATOR_D2:
155  case OPERATOR_D3:
156  case OPERATOR_D4:
157  case OPERATOR_D5:
158  case OPERATOR_D6:
159  case OPERATOR_D7:
160  case OPERATOR_D8:
161  case OPERATOR_D9:
162  case OPERATOR_D10:
163  {
164  auto derivativeOrder = getOperatorOrder(opType_) - 1;
165  Polynomials::legendreDerivativeValues(field_values_at_point, polyOrder_, x, derivativeOrder);
166 
167  // L_i is defined in terms of an integral of P_(i-1), so we need to shift the values by 1
168  if (numFields_ >= 3)
169  {
170  OutputScalar Pn_minus_one = field_values_at_point(1);
171  for (int fieldOrdinal=2; fieldOrdinal<numFields_; fieldOrdinal++)
172  {
173  OutputScalar Pn = field_values_at_point(fieldOrdinal);
174  field_values_at_point(fieldOrdinal) = Pn_minus_one;
175  Pn_minus_one = Pn;
176  }
177  }
178  if (numFields_ >= 1) field_values_at_point(0) = 0.0;
179  if (numFields_ >= 2) field_values_at_point(1) = 0.0;
180  // legendreDerivativeValues works on [-1,1], so no per-derivative scaling is necessary
181  // however, there is a factor of 0.5 that comes from the scaling of the Legendre polynomials prior to integration
182  // in the shiftedScaledIntegratedLegendreValues -- the first derivative of our integrated polynomials is 0.5 times the Legendre polynomial
183  break;
184  }
185  default:
186  // unsupported operator type
187  device_assert(false);
188  }
189 
190  // copy the values into the output container
191  for (int fieldOrdinal=0; fieldOrdinal<numFields_; fieldOrdinal++)
192  {
193  // access() allows us to write one line that applies both to gradient (for which outputValues has rank 3, but third rank has only one entry) and to value (rank 2)
194  output_.access(fieldOrdinal,pointOrdinal,0) = outputScaling * field_values_at_point(fieldOrdinal);
195  }
196  }
197 
198  // Provide the shared memory capacity.
199  // This function takes the team_size as an argument,
200  // which allows team_size-dependent allocations.
201  size_t team_shmem_size (int team_size) const
202  {
203  // we want to use shared memory to create a fast buffer that we can use for basis computations
204  size_t shmem_size = 0;
205  if (fad_size_output_ > 0)
206  shmem_size += OutputScratchView::shmem_size(numFields_, fad_size_output_);
207  else
208  shmem_size += OutputScratchView::shmem_size(numFields_);
209 
210  return shmem_size;
211  }
212  };
213 
231  template<typename DeviceType,
232  typename OutputScalar = double,
233  typename PointScalar = double,
234  bool defineVertexFunctions = true, // if defineVertexFunctions is true, first and second basis functions are x and 1-x. Otherwise, they are 1 and x.
235  bool useMinusOneToOneReferenceElement = true> // if useMinusOneToOneReferenceElement is true, basis is define on [-1,1]. Otherwise, [0,1].
237  : public Basis<DeviceType,OutputScalar,PointScalar>
238  {
239  public:
242 
243  using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost;
244  using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost;
245 
246  using OutputViewType = typename BasisBase::OutputViewType;
247  using PointViewType = typename BasisBase::PointViewType ;
248  using ScalarViewType = typename BasisBase::ScalarViewType;
249  protected:
250  int polyOrder_; // the maximum order of the polynomial
251  bool defineVertexFunctions_; // if true, first and second basis functions are x and 1-x. Otherwise, they are 1 and x.
252  EPointType pointType_;
253  public:
264  IntegratedLegendreBasis_HGRAD_LINE(int polyOrder, EPointType pointType=POINTTYPE_DEFAULT)
265  :
266  polyOrder_(polyOrder),
267  pointType_(pointType)
268  {
269  INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
270 
271  this->basisCardinality_ = polyOrder+1;
272  this->basisDegree_ = polyOrder;
273  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
274  this->basisType_ = BASIS_FEM_HIERARCHICAL;
275  this->basisCoordinates_ = COORDINATES_CARTESIAN;
276  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
277 
278  const int degreeLength = 1;
279  this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) line polynomial degree lookup", this->basisCardinality_, degreeLength);
280 
281  for (int i=0; i<this->basisCardinality_; i++)
282  {
283  // for H(grad) line, if defineVertexFunctions is false, first basis member is constant, second is first-degree, etc.
284  // if defineVertexFunctions is true, then the only difference is that the entry is also degree 1
285  this->fieldOrdinalPolynomialDegree_(i,0) = i;
286  }
287  if (defineVertexFunctions)
288  {
289  this->fieldOrdinalPolynomialDegree_(0,0) = 1;
290  }
291 
292  // initialize tags
293  {
294  const auto & cardinality = this->basisCardinality_;
295 
296  // Basis-dependent initializations
297  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
298  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
299  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
300  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
301 
302  OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
303 
304  if (defineVertexFunctions) {
305  {
306  const ordinal_type v0 = 0;
307  tagView(v0*tagSize+0) = 0; // vertex dof
308  tagView(v0*tagSize+1) = 0; // vertex id
309  tagView(v0*tagSize+2) = 0; // local dof id
310  tagView(v0*tagSize+3) = 1; // total number of dofs in this vertex
311 
312  const ordinal_type v1 = 1;
313  tagView(v1*tagSize+0) = 0; // vertex dof
314  tagView(v1*tagSize+1) = 1; // vertex id
315  tagView(v1*tagSize+2) = 0; // local dof id
316  tagView(v1*tagSize+3) = 1; // total number of dofs in this vertex
317 
318  const ordinal_type iend = cardinality - 2;
319  for (ordinal_type i=0;i<iend;++i) {
320  const auto e = i + 2;
321  tagView(e*tagSize+0) = 1; // edge dof
322  tagView(e*tagSize+1) = 0; // edge id
323  tagView(e*tagSize+2) = i; // local dof id
324  tagView(e*tagSize+3) = iend; // total number of dofs in this edge
325  }
326  }
327  } else {
328  for (ordinal_type i=0;i<cardinality;++i) {
329  tagView(i*tagSize+0) = 1; // edge dof
330  tagView(i*tagSize+1) = 0; // edge id
331  tagView(i*tagSize+2) = i; // local dof id
332  tagView(i*tagSize+3) = cardinality; // total number of dofs in this edge
333  }
334  }
335 
336  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
337  // tags are constructed on host
338  this->setOrdinalTagData(this->tagToOrdinal_,
339  this->ordinalToTag_,
340  tagView,
341  this->basisCardinality_,
342  tagSize,
343  posScDim,
344  posScOrd,
345  posDfOrd);
346  }
347  }
348 
353  const char* getName() const override {
354  return "Intrepid2_IntegratedLegendreBasis_HGRAD_LINE";
355  }
356 
359  virtual bool requireOrientation() const override {
360  return false;
361  }
362 
363  // since the getValues() below only overrides the FEM variant, we specify that
364  // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
365  // (It's an error to use the FVD variant on this basis.)
366  using BasisBase::getValues;
367 
386  virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
387  const EOperator operatorType = OPERATOR_VALUE ) const override
388  {
389  auto numPoints = inputPoints.extent_int(0);
390 
392 
393  FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
394 
395  const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
396  const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
397  const int vectorSize = std::max(outputVectorSize,pointVectorSize);
398  const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
399 
400  using ExecutionSpace = typename BasisBase::ExecutionSpace;
401 
402  auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
403  Kokkos::parallel_for( policy, functor, "Hierarchical_HGRAD_LINE_Functor");
404  }
405 
411  getHostBasis() const override {
412  using HostDeviceType = typename Kokkos::HostSpace::device_type;
414  return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
415  }
416  };
417 } // end namespace Intrepid2
418 
419 #endif /* Intrepid2_IntegratedLegendreBasis_HGRAD_LINE_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.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
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...
IntegratedLegendreBasis_HGRAD_LINE(int polyOrder, EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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...
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...
EPointType
Enumeration of types of point distributions in Intrepid.
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.
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
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.
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...
Functor for computing values for the IntegratedLegendreBasis_HGRAD_LINE class.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...