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),
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++;
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);
376 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
378 return Teuchos::rcp(
new HostBasisType(polyOrder_, pointType_) );
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.