Intrepid2
Intrepid2_LegendreBasis_HVOL_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_LegendreBasis_HVOL_TRI_h
50#define Intrepid2_LegendreBasis_HVOL_TRI_h
51
52#include <Kokkos_DynRankView.hpp>
53
54#include <Intrepid2_config.h>
55
56#include "Intrepid2_Basis.hpp"
59#include "Intrepid2_Utils.hpp"
60
61namespace Intrepid2
62{
68 template<class DeviceType, class OutputScalar, class PointScalar,
69 class OutputFieldType, class InputPointsType>
71 {
72 using ExecutionSpace = typename DeviceType::execution_space;
73 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
74 using OutputScratchView = 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 int numFields_, numPoints_;
87
88 size_t fad_size_output_;
89
90 Hierarchical_HVOL_TRI_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, int polyOrder)
91 : opType_(opType), output_(output), inputPoints_(inputPoints),
92 polyOrder_(polyOrder),
93 fad_size_output_(getScalarDimensionForView(output))
94 {
95 numFields_ = output.extent_int(0);
96 numPoints_ = output.extent_int(1);
97 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument, "point counts need to match!");
98 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument, "output field size does not match basis cardinality");
99 }
100
101 KOKKOS_INLINE_FUNCTION
102 void operator()( const TeamMember & teamMember ) const
103 {
104 auto pointOrdinal = teamMember.league_rank();
105 OutputScratchView legendre_field_values_at_point, jacobi_values_at_point;
106 if (fad_size_output_ > 0) {
107 legendre_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
108 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
109 }
110 else {
111 legendre_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
112 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
113 }
114
115 const auto & x = inputPoints_(pointOrdinal,0);
116 const auto & y = inputPoints_(pointOrdinal,1);
117
118 // write as barycentric coordinates:
119 const PointScalar lambda[3] = {1. - x - y, x, y};
120
121 switch (opType_)
122 {
123 case OPERATOR_VALUE:
124 {
125 // face functions
126 {
127 const PointScalar tLegendre = lambda[0] + lambda[1];
128 Polynomials::shiftedScaledLegendreValues(legendre_field_values_at_point, polyOrder_, lambda[1], tLegendre);
129
130 int fieldOrdinalOffset = 0;
131 const int max_ij_sum = polyOrder_;
132 for (int ij_sum=0; ij_sum<=max_ij_sum; ij_sum++)
133 {
134 for (int i=0; i<=ij_sum; i++)
135 {
136 const int j = ij_sum - i;
137 const auto & legendreValue = legendre_field_values_at_point(i);
138 const double alpha = i*2.0+1;
139
140 const PointScalar tJacobi = 1.0;// lambda[0] + lambda[1] + lambda[2];
141 Polynomials::shiftedScaledJacobiValues(jacobi_values_at_point, alpha, polyOrder_, lambda[2], tJacobi);
142
143 const auto & jacobiValue = jacobi_values_at_point(j);
144 output_(fieldOrdinalOffset,pointOrdinal) = legendreValue * jacobiValue;
145 fieldOrdinalOffset++;
146 }
147 }
148 }
149 } // end OPERATOR_VALUE
150 break;
151 default:
152 // unsupported operator type
153 device_assert(false);
154 }
155 }
156
157 // Provide the shared memory capacity.
158 // This function takes the team_size as an argument,
159 // which allows team_size-dependent allocations.
160 size_t team_shmem_size (int team_size) const
161 {
162 // TODO: edit this to match scratch that we actually need. (What's here is copied from H^1 basis on triangles…)
163 // we will use shared memory to create a fast buffer for basis computations
164 size_t shmem_size = 0;
165 if (fad_size_output_ > 0)
166 shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
167 else
168 shmem_size += 2 * OutputScratchView::shmem_size(polyOrder_ + 1);
169
170 return shmem_size;
171 }
172 };
173
184 template<typename DeviceType,
185 typename OutputScalar = double,
186 typename PointScalar = double>
188 : public Basis<DeviceType,OutputScalar,PointScalar>
189 {
190 public:
193
196
197 using typename BasisBase::OutputViewType;
198 using typename BasisBase::PointViewType;
199 using typename BasisBase::ScalarViewType;
200
201 using typename BasisBase::ExecutionSpace;
202
203 protected:
204 int polyOrder_; // the maximum order of the polynomial
205 EPointType pointType_;
206 public:
213 LegendreBasis_HVOL_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
214 :
215 polyOrder_(polyOrder),
216 pointType_(pointType)
217 {
218 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
219
220 this->basisCardinality_ = ((polyOrder+2) * (polyOrder+1)) / 2;
221 this->basisDegree_ = polyOrder;
222 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
223 this->basisType_ = BASIS_FEM_HIERARCHICAL;
224 this->basisCoordinates_ = COORDINATES_CARTESIAN;
225 this->functionSpace_ = FUNCTION_SPACE_HVOL;
226
227 const int degreeLength = 1;
228 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Legendre H(vol) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
229 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost("Legendre H(grad) line polynomial degree lookup", this->basisCardinality_, degreeLength);
230
231 int fieldOrdinalOffset = 0;
232 // **** face functions **** //
233 const int max_ij_sum = polyOrder;
234 for (int ij_sum=0; ij_sum<=max_ij_sum; ij_sum++)
235 {
236 for (int i=0; i<=ij_sum; i++)
237 {
238 const int j = ij_sum - i;
239 this->fieldOrdinalPolynomialDegree_ (fieldOrdinalOffset,0) = i+j;
240 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinalOffset,0) = i+j+1; // H^1 degree is one greater
241 fieldOrdinalOffset++;
242 }
243 }
244 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
245
246 // initialize tags
247 {
248 const auto & cardinality = this->basisCardinality_;
249
250 // Basis-dependent initializations
251 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
252 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
253 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
254 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
255
256 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
257 const int faceDim = 2;
258
259 for (ordinal_type i=0;i<cardinality;++i) {
260 tagView(i*tagSize+0) = faceDim; // face dimension
261 tagView(i*tagSize+1) = 0; // face id
262 tagView(i*tagSize+2) = i; // local dof id
263 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
264 }
265
266 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
267 // tags are constructed on host
269 this->ordinalToTag_,
270 tagView,
271 this->basisCardinality_,
272 tagSize,
273 posScDim,
274 posScOrd,
275 posDfOrd);
276 }
277 }
278
283 const char* getName() const override {
284 return "Intrepid2_LegendreBasis_HVOL_TRI";
285 }
286
289 virtual bool requireOrientation() const override {
290 return (this->getDegree() > 2);
291 }
292
293 // since the getValues() below only overrides the FEM variant, we specify that
294 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
295 // (It's an error to use the FVD variant on this basis.)
297
316 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
317 const EOperator operatorType = OPERATOR_VALUE ) const override
318 {
319 auto numPoints = inputPoints.extent_int(0);
320
322
323 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
324
325 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
326 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
327 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
328 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
329
330 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
331 Kokkos::parallel_for("Hierarchical_HVOL_TRI_Functor", policy, functor);
332 }
333
343 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
344 if(subCellDim == 1) {
345 return Teuchos::rcp(new
347 (this->basisDegree_));
348 }
349 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
350 }
351
357 getHostBasis() const override {
358 using HostDeviceType = typename Kokkos::HostSpace::device_type;
360 return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
361 }
362 };
363} // end namespace Intrepid2
364
365#endif /* Intrepid2_LegendreBasis_HVOL_TRI_h */
Header file for the abstract base class Intrepid2::Basis.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
H(grad) basis on the line based on integrated Legendre polynomials.
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_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.
ordinal_type getDegree() const
Returns the degree of 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 integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line.
Basis defining Legendre basis on the line, a polynomial subspace of H(vol) on the line: extension to ...
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.
const char * getName() const override
Returns basis name.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
LegendreBasis_HVOL_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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...
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
virtual bool requireOrientation() const override
True if orientation is required.
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_TRI class.