49 #ifndef Intrepid2_LegendreBasis_HVOL_TET_h
50 #define Intrepid2_LegendreBasis_HVOL_TET_h
52 #include <Kokkos_DynRankView.hpp>
54 #include <Intrepid2_config.h>
69 template<
class DeviceType,
class OutputScalar,
class PointScalar,
70 class OutputFieldType,
class InputPointsType>
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>>;
78 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
79 using TeamMember =
typename TeamPolicy::member_type;
83 OutputFieldType output_;
84 InputPointsType inputPoints_;
87 int numFields_, numPoints_;
89 size_t fad_size_output_;
92 : opType_(opType), output_(output), inputPoints_(inputPoints),
93 polyOrder_(polyOrder),
94 fad_size_output_(getScalarDimensionForView(output))
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");
102 KOKKOS_INLINE_FUNCTION
103 void operator()(
const TeamMember & teamMember )
const
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_);
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);
128 const auto & x = inputPoints_(pointOrdinal,0);
129 const auto & y = inputPoints_(pointOrdinal,1);
130 const auto & z = inputPoints_(pointOrdinal,2);
133 const PointScalar lambda[4] = {1. - x - y - z, x, y, z};
141 const PointScalar tLegendre = lambda[0] + lambda[1];
142 Polynomials::shiftedScaledLegendreValues(P, polyOrder_, lambda[1], tLegendre);
144 int fieldOrdinalOffset = 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++)
153 for (
int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
155 for (
int i=min_i; i <= totalPolyOrder_ij-min_j; i++)
157 const int j = totalPolyOrder_ij - i;
158 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
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);
165 const double alpha2 = 2. * (i + j + 1.);
166 const PointScalar tJacobi2 = 1.0;
167 const PointScalar & xJacobi2 = lambda[3];
168 Polynomials::shiftedScaledJacobiValues(P_2ipjp1, alpha2, polyOrder_, xJacobi2, tJacobi2);
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);
174 output_(fieldOrdinalOffset,pointOrdinal) = P_i * P_2p1_j * P_2ipjp1_k;
175 fieldOrdinalOffset++;
184 device_assert(
false);
191 size_t team_shmem_size (
int team_size)
const
194 size_t shmem_size = 0;
195 if (fad_size_output_ > 0)
196 shmem_size += 3 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
198 shmem_size += 3 * OutputScratchView::shmem_size(polyOrder_ + 1);
214 template<
typename DeviceType,
215 typename OutputScalar = double,
216 typename PointScalar =
double>
218 :
public Basis<DeviceType,OutputScalar,PointScalar>
234 EPointType pointType_;
244 polyOrder_(polyOrder),
245 pointType_(pointType)
247 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,
"PointType not supported");
251 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<> >() );
256 const int degreeLength = 1;
259 int fieldOrdinalOffset = 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++)
268 for (
int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_j; totalPolyOrder_ij++)
270 for (
int i=min_i; i <= totalPolyOrder_ij-min_j; i++)
272 const int j = totalPolyOrder_ij - i;
273 const int k = totalPolyOrder_ijk - totalPolyOrder_ij;
276 fieldOrdinalOffset++;
280 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
287 const ordinal_type tagSize = 4;
288 const ordinal_type posScDim = 0;
289 const ordinal_type posScOrd = 1;
290 const ordinal_type posDfOrd = 2;
293 const int volumeDim = 3;
295 for (ordinal_type i=0;i<cardinality;++i) {
296 tagView(i*tagSize+0) = volumeDim;
297 tagView(i*tagSize+1) = 0;
298 tagView(i*tagSize+2) = i;
299 tagView(i*tagSize+3) = cardinality;
320 return "Intrepid2_LegendreBasis_HVOL_TET";
353 const EOperator operatorType = OPERATOR_VALUE )
const override
355 auto numPoints = inputPoints.extent_int(0);
359 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
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;
366 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
367 Kokkos::parallel_for(
"Hierarchical_HVOL_TET_Functor", policy , functor);
374 virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>
376 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
378 return Teuchos::rcp(
new HostBasisType(polyOrder_, pointType_) );
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
const char * getName() const override
Returns basis name.
LegendreBasis_HVOL_TET(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual void getValues(const ExecutionSpace &, OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
EBasis basisType_
Type of the basis.
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...
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
virtual bool requireOrientation() const override
True if orientation is required.
Functor for computing values for the LegendreBasis_HVOL_TET class.
Free functions, callable from device code, that implement various polynomials useful in basis definit...
EFunctionSpace functionSpace_
The function space in which the basis is defined.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
Header function for Intrepid2::Util class and other utility functions.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
H(vol) basis on the line based on Legendre polynomials.
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.
Basis defining Legendre basis on the line, a polynomial subspace of H(vol) on the line: extension to ...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
ordinal_type getDegree() const
Returns the degree of the basis.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
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.
H(vol) basis on the triangle based on integrated Legendre polynomials.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Header file for the abstract base class Intrepid2::Basis.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.