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_DynRankView.hpp>
53
54#include <Intrepid2_config.h>
55
56// Sacado header that defines some fad sizing…
57#ifdef HAVE_INTREPID2_SACADO
58#include <KokkosExp_View_Fad.hpp>
59#endif
60
61#include "Intrepid2_Basis.hpp"
64
65namespace 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
193
194 using typename BasisBase::OutputViewType;
195 using typename BasisBase::PointViewType;
196 using typename BasisBase::ScalarViewType;
197
198 using typename BasisBase::ExecutionSpace;
199
200 protected:
201 int polyOrder_; // the maximum order of the polynomial
202 EPointType pointType_;
203 public:
212 LegendreBasis_HVOL_LINE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
213 :
214 polyOrder_(polyOrder),
215 pointType_(pointType)
216 {
217 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
218 this->basisCardinality_ = polyOrder+1;
219 this->basisDegree_ = polyOrder;
220 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
221 this->basisType_ = BASIS_FEM_HIERARCHICAL;
222 this->basisCoordinates_ = COORDINATES_CARTESIAN;
223 this->functionSpace_ = FUNCTION_SPACE_HVOL;
224
225 const int degreeLength = 1;
226 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) line polynomial degree lookup", this->basisCardinality_, degreeLength);
227 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(grad) line polynomial degree lookup", this->basisCardinality_, degreeLength);
228
229 for (int i=0; i<this->basisCardinality_; i++)
230 {
231 // for H(vol) line, first basis member is constant, second is first-degree, etc.
232 this->fieldOrdinalPolynomialDegree_ (i,0) = i;
233 this->fieldOrdinalH1PolynomialDegree_(i,0) = i+1; // H^1 degree is one greater
234 }
235
236 // initialize tags
237 {
238 const auto & cardinality = this->basisCardinality_;
239
240 // Basis-dependent initializations
241 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
242 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
243 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
244 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
245
246 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
247
248 for (ordinal_type i=0;i<cardinality;++i) {
249 tagView(i*tagSize+0) = 1; // edge dof
250 tagView(i*tagSize+1) = 0; // edge id
251 tagView(i*tagSize+2) = i; // local dof id
252 tagView(i*tagSize+3) = cardinality; // total number of dofs in this edge
253 }
254
255 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
256 // tags are constructed on host
258 this->ordinalToTag_,
259 tagView,
260 this->basisCardinality_,
261 tagSize,
262 posScDim,
263 posScOrd,
264 posDfOrd);
265 }
266 }
267
268 // since the getValues() below only overrides the FEM variant, we specify that
269 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
270 // (It's an error to use the FVD variant on this basis.)
272
277 virtual
278 const char*
279 getName() const override {
280 return "Intrepid2_LegendreBasis_HVOL_LINE";
281 }
282
301 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
302 const EOperator operatorType = OPERATOR_VALUE ) const override
303 {
304 auto numPoints = inputPoints.extent_int(0);
305
307
308 FunctorType functor(outputValues, inputPoints, polyOrder_, operatorType);
309
310 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
311 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
312 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
313 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
314
315 auto policy =
316 Kokkos::TeamPolicy<ExecutionSpace>(numPoints, teamSize, vectorSize);
317 Kokkos::parallel_for("Hierarchical_HVOL_LINE_Functor", policy , functor);
318 }
319
325 getHostBasis() const override {
326 using HostDeviceType = typename Kokkos::HostSpace::device_type;
328 return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
329 }
330 };
331} // end namespace Intrepid2
332
333#endif /* Intrepid2_LegendreBasis_HVOL_LINE_h */
Header file for the abstract base class Intrepid2::Basis.
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
Implementation of an assert that can safely be called from device code.
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
Free functions, callable from device code, that implement various polynomials useful in basis definit...
KOKKOS_INLINE_FUNCTION constexpr unsigned getScalarDimensionForView(const ViewType &view)
Returns the size of the Scalar dimension for the View. This is 0 for non-AD types....
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
EBasis basisType_
Type of the basis.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
void setOrdinalTagData(OrdinalTypeView3D &tagToOrdinal, OrdinalTypeView2D &ordinalToTag, const OrdinalTypeView1D tags, const ordinal_type basisCard, const ordinal_type tagSize, const ordinal_type posScDim, const ordinal_type posScOrd, const ordinal_type posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
OrdinalTypeArray2DHost fieldOrdinalH1PolynomialDegree_
H^1 polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
EFunctionSpace functionSpace_
The function space in which the basis is defined.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
virtual const char * getName() const override
Returns basis name.
virtual BasisPtr< typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
LegendreBasis_HVOL_LINE(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Functor for computing values for the LegendreBasis_HVOL_LINE class.