Intrepid2
Intrepid2_LegendreBasis_HVOL_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_LegendreBasis_HVOL_LINE_h
50 #define Intrepid2_LegendreBasis_HVOL_LINE_h
51 
52 #include <Kokkos_View.hpp>
53 #include <Kokkos_DynRankView.hpp>
54 
55 #include <Intrepid2_config.h>
56 
57 // Sacado header that defines some fad sizing…
58 #ifdef HAVE_INTREPID2_SACADO
59 #include <KokkosExp_View_Fad.hpp>
60 #endif
61 
64 
65 namespace Intrepid2
66 {
72  template<class DeviceType, class OutputScalar, class PointScalar,
73  class OutputFieldType, class InputPointsType>
75  {
76  using ExecutionSpace = typename DeviceType::execution_space;
77  using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
78  using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
79  using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
80 
81  using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
82  using TeamMember = typename TeamPolicy::member_type;
83 
84  OutputFieldType output_; // F,P
85  InputPointsType inputPoints_; // P,D
86 
87  int polyOrder_;
88  int numFields_, numPoints_;
89 
90  EOperator op_;
91 
92  size_t fad_size_output_;
93 
94  Hierarchical_HVOL_LINE_Functor(OutputFieldType output, InputPointsType inputPoints,
95  int polyOrder, EOperator op)
96  : output_(output), inputPoints_(inputPoints), polyOrder_(polyOrder), op_(op),
97  fad_size_output_(getScalarDimensionForView(output))
98  {
99  numFields_ = output.extent_int(0);
100  numPoints_ = output.extent_int(1);
101  INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
102  INTREPID2_TEST_FOR_EXCEPTION(numFields_ != polyOrder_+1, std::invalid_argument, "output field size does not match basis cardinality");
103  }
104 
105  KOKKOS_INLINE_FUNCTION
106  void operator()( const TeamMember & teamMember ) const
107  {
108  auto pointOrdinal = teamMember.league_rank();
109  OutputScratchView field_values_at_point;
110  if (fad_size_output_ > 0) {
111  field_values_at_point = OutputScratchView(teamMember.team_shmem(), numFields_, fad_size_output_);
112  }
113  else {
114  field_values_at_point = OutputScratchView(teamMember.team_shmem(), numFields_);
115  }
116 
117  const PointScalar x = inputPoints_(pointOrdinal,0);
118 
119  switch (op_)
120  {
121  case OPERATOR_VALUE:
122  Polynomials::legendreValues(field_values_at_point, polyOrder_, x);
123 
124  // note that because legendreValues determines field values recursively, there is not much
125  // opportunity at that level for further parallelism
126  break;
127  case OPERATOR_GRAD:
128  case OPERATOR_D1:
129  case OPERATOR_D2:
130  case OPERATOR_D3:
131  case OPERATOR_D4:
132  case OPERATOR_D5:
133  case OPERATOR_D6:
134  case OPERATOR_D7:
135  case OPERATOR_D8:
136  case OPERATOR_D9:
137  case OPERATOR_D10:
138  {
139  auto derivativeOrder = getOperatorOrder(op_);
140  Polynomials::legendreDerivativeValues(field_values_at_point, polyOrder_, x, derivativeOrder);
141  break;
142  }
143  default:
144  // unsupported operator type
145  device_assert(false);
146  }
147  // copy the values into the output container
148  for (int fieldOrdinal=0; fieldOrdinal<numFields_; fieldOrdinal++)
149  {
150  output_.access(fieldOrdinal,pointOrdinal,0) = field_values_at_point(fieldOrdinal);
151  }
152  }
153 
154  // Provide the shared memory capacity.
155  // This function takes the team_size as an argument,
156  // which allows team_size-dependent allocations.
157  size_t team_shmem_size (int team_size) const
158  {
159  // we want to use shared memory to create a fast buffer that we can use for basis computations
160  size_t shmem_size = 0;
161  if (fad_size_output_ > 0)
162  shmem_size += OutputScratchView::shmem_size(numFields_, fad_size_output_);
163  else
164  shmem_size += OutputScratchView::shmem_size(numFields_);
165 
166  return shmem_size;
167  }
168  };
169 
181  template<typename DeviceType,
182  typename OutputScalar = double,
183  typename PointScalar = double>
185  : public Basis<DeviceType,OutputScalar,PointScalar>
186  {
187  public:
190 
191  using OrdinalTypeArray1DHost = typename BasisBase::OrdinalTypeArray1DHost;
192  using OrdinalTypeArray2DHost = typename BasisBase::OrdinalTypeArray2DHost;
193 
194  using OutputViewType = typename BasisBase::OutputViewType;
195  using PointViewType = typename BasisBase::PointViewType ;
196  using ScalarViewType = typename BasisBase::ScalarViewType;
197  protected:
198  int polyOrder_; // the maximum order of the polynomial
199  EPointType pointType_;
200  public:
209  LegendreBasis_HVOL_LINE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
210  :
211  polyOrder_(polyOrder),
212  pointType_(pointType)
213  {
214  INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
215  this->basisCardinality_ = polyOrder+1;
216  this->basisDegree_ = polyOrder;
217  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
218  this->basisType_ = BASIS_FEM_HIERARCHICAL;
219  this->basisCoordinates_ = COORDINATES_CARTESIAN;
220  this->functionSpace_ = FUNCTION_SPACE_HVOL;
221 
222  const int degreeLength = 1;
223  this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) line polynomial degree lookup", this->basisCardinality_, degreeLength);
224 
225  for (int i=0; i<this->basisCardinality_; i++)
226  {
227  // for H(vol) line, first basis member is constant, second is first-degree, etc.
228  this->fieldOrdinalPolynomialDegree_(i,0) = i;
229  }
230 
231  // initialize tags
232  {
233  const auto & cardinality = this->basisCardinality_;
234 
235  // Basis-dependent initializations
236  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
237  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
238  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
239  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
240 
241  OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
242 
243  for (ordinal_type i=0;i<cardinality;++i) {
244  tagView(i*tagSize+0) = 1; // edge dof
245  tagView(i*tagSize+1) = 0; // edge id
246  tagView(i*tagSize+2) = i; // local dof id
247  tagView(i*tagSize+3) = cardinality; // total number of dofs in this edge
248  }
249 
250  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
251  // tags are constructed on host
252  this->setOrdinalTagData(this->tagToOrdinal_,
253  this->ordinalToTag_,
254  tagView,
255  this->basisCardinality_,
256  tagSize,
257  posScDim,
258  posScOrd,
259  posDfOrd);
260  }
261  }
262 
263  // since the getValues() below only overrides the FEM variant, we specify that
264  // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
265  // (It's an error to use the FVD variant on this basis.)
266  using BasisBase::getValues;
267 
272  virtual
273  const char*
274  getName() const override {
275  return "Intrepid2_LegendreBasis_HVOL_LINE";
276  }
277 
296  virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
297  const EOperator operatorType = OPERATOR_VALUE ) const override
298  {
299  auto numPoints = inputPoints.extent_int(0);
300 
302 
303  FunctorType functor(outputValues, inputPoints, polyOrder_, operatorType);
304 
305  const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
306  const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
307  const int vectorSize = std::max(outputVectorSize,pointVectorSize);
308  const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
309 
310  using ExecutionSpace = typename BasisBase::ExecutionSpace;
311 
312  auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
313  Kokkos::parallel_for( policy , functor, "Hierarchical_HVOL_LINE_Functor");
314  }
315 
321  getHostBasis() const override {
322  using HostDeviceType = typename Kokkos::HostSpace::device_type;
324  return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
325  }
326  };
327 } // end namespace Intrepid2
328 
329 #endif /* Intrepid2_LegendreBasis_HVOL_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
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 const char * getName() const override
Returns basis name.
Functor for computing values for the LegendreBasis_HVOL_LINE class.
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...
LegendreBasis_HVOL_LINE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Implementation of an assert that can safely be called from device code.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line...
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 outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
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...
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...