Intrepid2
Intrepid2_LegendreBasis_HVOL_TET.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_TET_h
50#define Intrepid2_LegendreBasis_HVOL_TET_h
51
52#include <Kokkos_DynRankView.hpp>
53
54#include <Intrepid2_config.h>
55
56#include "Intrepid2_Basis.hpp"
60#include "Intrepid2_Utils.hpp"
61
62namespace Intrepid2
63{
69 template<class DeviceType, class OutputScalar, class PointScalar,
70 class OutputFieldType, class InputPointsType>
72 {
73 using ExecutionSpace = typename DeviceType::execution_space;
74 using ScratchSpace = typename ExecutionSpace::scratch_memory_space;
75 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
76 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
77
78 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
79 using TeamMember = typename TeamPolicy::member_type;
80
81 EOperator opType_;
82
83 OutputFieldType output_; // F,P
84 InputPointsType inputPoints_; // P,D
85
86 int polyOrder_;
87 int numFields_, numPoints_;
88
89 size_t fad_size_output_;
90
91 Hierarchical_HVOL_TET_Functor(EOperator opType, OutputFieldType output, InputPointsType inputPoints, int polyOrder)
92 : opType_(opType), output_(output), inputPoints_(inputPoints),
93 polyOrder_(polyOrder),
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)*(polyOrder_+2)*(polyOrder_+3)/6, 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 // values are product of [P_i](lambda_0,lambda_1), [P^{2i+1}_j](lambda_0 + lambda_1, lambda_2), and [P^{2*(i+j+1)}_k](1-lambda_3,lambda_3),
106 // times ((grad lambda_1) x (grad lambda_2)) \cdot (grad lambda_3).
107 // For the canonical orientation (all we support), the last term evaluates to 1, and
108 // lambda_0 = 1 - x - y - z
109 // lambda_1 = x
110 // lambda_2 = y
111 // lambda_3 = z
112 // [P_i](lambda_0, lambda_1) = P_i(lambda_1; lambda_0 + lambda_1) = P_i(x; 1 - y - z) -- a shifted, scaled Legendre function
113 // [P^{2i+1}_j](lambda_0 + lambda_1, lambda_2) = P^{2i+1}_j(lambda_2; lambda_0 + lambda_1 + lambda_2) = P^{2i+1}_j(y; 1 - z) -- a shifted, scaled Jacobi function
114 // [P^{2*(i+j+1)}_k](1-lambda_3,lambda_3) = P^{2*(i+j+1)}_k(lambda_3; 1) = P^{2*(i+j+1)}_k(z; 1) -- another shifted, scaled Jacobi function
115 auto pointOrdinal = teamMember.league_rank();
116 OutputScratchView P, P_2p1, P_2ipjp1;
117 if (fad_size_output_ > 0) {
118 P = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
119 P_2p1 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
120 P_2ipjp1 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
121 }
122 else {
123 P = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
124 P_2p1 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
125 P_2ipjp1 = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
126 }
127
128 const auto & x = inputPoints_(pointOrdinal,0);
129 const auto & y = inputPoints_(pointOrdinal,1);
130 const auto & z = inputPoints_(pointOrdinal,2);
131
132 // write as barycentric coordinates:
133 const PointScalar lambda[4] = {1. - x - y - z, x, y, z};
134
135 switch (opType_)
136 {
137 case OPERATOR_VALUE:
138 {
139 // face functions
140 {
141 const PointScalar tLegendre = lambda[0] + lambda[1];
142 Polynomials::shiftedScaledLegendreValues(P, polyOrder_, lambda[1], tLegendre);
143
144 int fieldOrdinalOffset = 0;
145
146 const int min_i = 0;
147 const int min_j = 0;
148 const int min_k = 0;
149 const int min_ij = min_i + min_j;
150 const int min_ijk = min_ij + min_k;
151 for (int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= polyOrder_; totalPolyOrder_ijk++)
152 {
153 for (int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
154 {
155 for (int i=min_i; i <= totalPolyOrder_ij-min_j; i++)
156 {
157 const int j = totalPolyOrder_ij - i;
158 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
159
160 const double alpha1 = i * 2.0 + 1.;
161 const PointScalar tJacobi1 = lambda[0] + lambda[1] + lambda[2];
162 const PointScalar & xJacobi1 = lambda[2];
163 Polynomials::shiftedScaledJacobiValues(P_2p1, alpha1, polyOrder_, xJacobi1, tJacobi1);
164
165 const double alpha2 = 2. * (i + j + 1.);
166 const PointScalar tJacobi2 = 1.0; // 1 - lambda[3] + lambda[3]
167 const PointScalar & xJacobi2 = lambda[3];
168 Polynomials::shiftedScaledJacobiValues(P_2ipjp1, alpha2, polyOrder_, xJacobi2, tJacobi2);
169
170 const auto & P_i = P(i);
171 const auto & P_2p1_j = P_2p1(j);
172 const auto & P_2ipjp1_k = P_2ipjp1(k);
173
174 output_(fieldOrdinalOffset,pointOrdinal) = P_i * P_2p1_j * P_2ipjp1_k;
175 fieldOrdinalOffset++;
176 }
177 }
178 }
179 }
180 } // end OPERATOR_VALUE
181 break;
182 default:
183 // unsupported operator type
184 device_assert(false);
185 }
186 }
187
188 // Provide the shared memory capacity.
189 // This function takes the team_size as an argument,
190 // which allows team_size-dependent allocations.
191 size_t team_shmem_size (int team_size) const
192 {
193 // we will use shared memory to create a fast buffer for basis computations
194 size_t shmem_size = 0;
195 if (fad_size_output_ > 0)
196 shmem_size += 3 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
197 else
198 shmem_size += 3 * OutputScratchView::shmem_size(polyOrder_ + 1);
199
200 return shmem_size;
201 }
202 };
203
214 template<typename DeviceType,
215 typename OutputScalar = double,
216 typename PointScalar = double>
218 : public Basis<DeviceType,OutputScalar,PointScalar>
219 {
220 public:
222
225
226 using typename BasisBase::OutputViewType;
227 using typename BasisBase::PointViewType;
228 using typename BasisBase::ScalarViewType;
229
230 using typename BasisBase::ExecutionSpace;
231
232 protected:
233 int polyOrder_; // the maximum order of the polynomial
234 EPointType pointType_;
235 public:
242 LegendreBasis_HVOL_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
243 :
244 polyOrder_(polyOrder),
245 pointType_(pointType)
246 {
247 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,"PointType not supported");
248
249 this->basisCardinality_ = ((polyOrder+3) * (polyOrder+2) * (polyOrder+1)) / 6;
250 this->basisDegree_ = polyOrder;
251 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
252 this->basisType_ = BASIS_FEM_HIERARCHICAL;
253 this->basisCoordinates_ = COORDINATES_CARTESIAN;
254 this->functionSpace_ = FUNCTION_SPACE_HVOL;
255
256 const int degreeLength = 1;
257 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost("Integrated Legendre H(vol) triangle polynomial degree lookup", this->basisCardinality_, degreeLength);
258
259 int fieldOrdinalOffset = 0;
260 // **** volume/interior functions **** //
261 const int min_i = 0;
262 const int min_j = 0;
263 const int min_k = 0;
264 const int min_ij = min_i + min_j;
265 const int min_ijk = min_ij + min_k;
266 for (int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= polyOrder_; totalPolyOrder_ijk++)
267 {
268 for (int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
269 {
270 for (int i=min_i; i <= totalPolyOrder_ij-min_j; i++)
271 {
272 const int j = totalPolyOrder_ij - i;
273 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
274
275 this->fieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = i+j+k;
276 fieldOrdinalOffset++;
277 }
278 }
279 }
280 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
281
282 // initialize tags
283 {
284 const auto & cardinality = this->basisCardinality_;
285
286 // Basis-dependent initializations
287 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
288 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
289 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
290 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
291
292 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
293 const int volumeDim = 3;
294
295 for (ordinal_type i=0;i<cardinality;++i) {
296 tagView(i*tagSize+0) = volumeDim; // volume dimension
297 tagView(i*tagSize+1) = 0; // volume id
298 tagView(i*tagSize+2) = i; // local dof id
299 tagView(i*tagSize+3) = cardinality; // total number of dofs on this face
300 }
301
302 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
303 // tags are constructed on host
305 this->ordinalToTag_,
306 tagView,
307 this->basisCardinality_,
308 tagSize,
309 posScDim,
310 posScOrd,
311 posDfOrd);
312 }
313 }
314
319 const char* getName() const override {
320 return "Intrepid2_LegendreBasis_HVOL_TET";
321 }
322
325 virtual bool requireOrientation() const override {
326 return (this->getDegree() > 2);
327 }
328
329 // since the getValues() below only overrides the FEM variant, we specify that
330 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
331 // (It's an error to use the FVD variant on this basis.)
333
352 virtual void getValues( OutputViewType outputValues, const PointViewType inputPoints,
353 const EOperator operatorType = OPERATOR_VALUE ) const override
354 {
355 auto numPoints = inputPoints.extent_int(0);
356
358
359 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
360
361 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
362 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
363 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
364 const int teamSize = 1; // because of the way the basis functions are computed, we don't have a second level of parallelism...
365
366 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
367 Kokkos::parallel_for("Hierarchical_HVOL_TET_Functor", policy , functor);
368 }
369
375 getHostBasis() const override {
376 using HostDeviceType = typename Kokkos::HostSpace::device_type;
378 return Teuchos::rcp( new HostBasisType(polyOrder_, pointType_) );
379 }
380 };
381} // end namespace Intrepid2
382
383#endif /* Intrepid2_LegendreBasis_HVOL_TET_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(vol) basis on the line based on Legendre polynomials.
H(vol) basis on the triangle 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
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 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.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
const char * getName() const override
Returns basis name.
LegendreBasis_HVOL_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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.
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 LegendreBasis_HVOL_TET class.