49 #ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h 50 #define Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h 52 #include <Kokkos_View.hpp> 53 #include <Kokkos_DynRankView.hpp> 55 #include <Intrepid2_config.h> 67 template<
class DeviceType,
class OutputScalar,
class PointScalar,
68 class OutputFieldType,
class InputPointsType>
71 using ExecutionSpace =
typename DeviceType::execution_space;
72 using ScratchSpace =
typename ExecutionSpace::scratch_memory_space;
73 using OutputScratchView = Kokkos::View<OutputScalar*,ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
74 using PointScratchView = Kokkos::View<PointScalar*, ScratchSpace,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
76 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
77 using TeamMember =
typename TeamPolicy::member_type;
81 OutputFieldType output_;
82 InputPointsType inputPoints_;
85 bool defineVertexFunctions_;
86 int numFields_, numPoints_;
88 size_t fad_size_output_;
90 static const int numVertices = 3;
91 static const int numEdges = 3;
92 const int edge_start_[numEdges] = {0,1,0};
93 const int edge_end_[numEdges] = {1,2,2};
96 int polyOrder,
bool defineVertexFunctions)
97 : opType_(opType), output_(output), inputPoints_(inputPoints),
98 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
101 numFields_ = output.extent_int(0);
102 numPoints_ = output.extent_int(1);
103 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
104 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument,
"output field size does not match basis cardinality");
107 KOKKOS_INLINE_FUNCTION
108 void operator()(
const TeamMember & teamMember )
const 110 auto pointOrdinal = teamMember.league_rank();
111 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
112 if (fad_size_output_ > 0) {
113 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
114 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
115 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
116 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
119 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
120 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
121 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
122 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
125 const auto & x = inputPoints_(pointOrdinal,0);
126 const auto & y = inputPoints_(pointOrdinal,1);
129 const PointScalar lambda[3] = {1. - x - y, x, y};
130 const PointScalar lambda_dx[3] = {-1., 1., 0.};
131 const PointScalar lambda_dy[3] = {-1., 0., 1.};
133 const int num1DEdgeFunctions = polyOrder_ - 1;
140 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
142 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
144 if (!defineVertexFunctions_)
148 output_(0,pointOrdinal) = 1.0;
152 int fieldOrdinalOffset = 3;
153 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
155 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
156 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
158 Polynomials::shiftedScaledIntegratedLegendreValues(edge_field_values_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
159 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
162 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = edge_field_values_at_point(edgeFunctionOrdinal+2);
164 fieldOrdinalOffset += num1DEdgeFunctions;
170 const double jacobiScaling = 1.0;
172 for (
int i=2; i<polyOrder_; i++)
174 const int edgeBasisOrdinal = i+numVertices-2;
175 const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
176 const double alpha = i*2.0;
178 Polynomials::integratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-2, lambda[2], jacobiScaling);
179 for (
int j=1; i+j <= polyOrder_; j++)
181 const auto & jacobiValue = jacobi_values_at_point(j);
182 output_(fieldOrdinalOffset,pointOrdinal) = edgeValue * jacobiValue;
183 fieldOrdinalOffset++;
193 if (defineVertexFunctions_)
197 output_(0,pointOrdinal,0) = -1.0;
198 output_(0,pointOrdinal,1) = -1.0;
204 output_(0,pointOrdinal,0) = 0.0;
205 output_(0,pointOrdinal,1) = 0.0;
208 output_(1,pointOrdinal,0) = 1.0;
209 output_(1,pointOrdinal,1) = 0.0;
211 output_(2,pointOrdinal,0) = 0.0;
212 output_(2,pointOrdinal,1) = 1.0;
215 int fieldOrdinalOffset = 3;
227 auto & P_i_minus_1 = edge_field_values_at_point;
228 auto & L_i_dt = jacobi_values_at_point;
229 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
231 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
232 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
234 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
235 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
236 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
237 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
239 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
240 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
241 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
244 const int i = edgeFunctionOrdinal+2;
245 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
246 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
248 fieldOrdinalOffset += num1DEdgeFunctions;
269 auto & P_2i_j_minus_1 = edge_field_values_at_point;
270 auto & L_2i_j_dt = jacobi_values_at_point;
271 auto & L_i = other_values_at_point;
272 auto & L_2i_j = other_values2_at_point;
275 const double jacobiScaling = 1.0;
277 for (
int i=2; i<polyOrder_; i++)
280 const int edgeBasisOrdinal = i+numVertices-2;
281 const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
282 const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
284 const double alpha = i*2.0;
286 Polynomials::shiftedScaledIntegratedLegendreValues(L_i, polyOrder_, lambda[1], lambda[0]+lambda[1]);
287 Polynomials::integratedJacobiValues_dt( L_2i_j_dt, alpha, polyOrder_, lambda[2], jacobiScaling);
288 Polynomials::integratedJacobiValues ( L_2i_j, alpha, polyOrder_, lambda[2], jacobiScaling);
289 Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1, alpha, polyOrder_-1, lambda[2], jacobiScaling);
291 const auto & s0_dx = lambda_dx[0];
292 const auto & s0_dy = lambda_dy[0];
293 const auto & s1_dx = lambda_dx[1];
294 const auto & s1_dy = lambda_dy[1];
295 const auto & s2_dx = lambda_dx[2];
296 const auto & s2_dy = lambda_dy[2];
298 for (
int j=1; i+j <= polyOrder_; j++)
300 const OutputScalar basisValue_dx = L_2i_j(j) * grad_L_i_dx + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dx + L_2i_j_dt(j) * (s0_dx + s1_dx + s2_dx));
301 const OutputScalar basisValue_dy = L_2i_j(j) * grad_L_i_dy + L_i(i) * (P_2i_j_minus_1(j-1) * s2_dy + L_2i_j_dt(j) * (s0_dy + s1_dy + s2_dy));
303 output_(fieldOrdinalOffset,pointOrdinal,0) = basisValue_dx;
304 output_(fieldOrdinalOffset,pointOrdinal,1) = basisValue_dy;
305 fieldOrdinalOffset++;
320 INTREPID2_TEST_FOR_ABORT(
true,
321 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
324 device_assert(
false);
331 size_t team_shmem_size (
int team_size)
const 334 size_t shmem_size = 0;
335 if (fad_size_output_ > 0)
336 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
338 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
361 template<
typename DeviceType,
362 typename OutputScalar = double,
363 typename PointScalar = double,
364 bool defineVertexFunctions =
true>
366 :
public Basis<DeviceType,OutputScalar,PointScalar>
393 polyOrder_(polyOrder),
394 pointType_(pointType)
396 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,
"PointType not supported");
400 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
405 const int degreeLength = 1;
408 int fieldOrdinalOffset = 0;
411 const int numFunctionsPerVertex = 1;
412 const int numVertexFunctions = numVertices * numFunctionsPerVertex;
413 for (
int i=0; i<numVertexFunctions; i++)
419 if (!defineVertexFunctions)
423 fieldOrdinalOffset += numVertexFunctions;
426 const int numFunctionsPerEdge = polyOrder - 1;
428 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
430 for (
int i=0; i<numFunctionsPerEdge; i++)
434 fieldOrdinalOffset += numFunctionsPerEdge;
438 const int max_ij_sum = polyOrder;
439 for (
int i=2; i<max_ij_sum; i++)
441 for (
int j=1; i+j<=max_ij_sum; j++)
444 fieldOrdinalOffset++;
447 const int numFaces = 1;
448 const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
449 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
456 const ordinal_type tagSize = 4;
457 const ordinal_type posScDim = 0;
458 const ordinal_type posScOrd = 1;
459 const ordinal_type posDfOrd = 2;
461 OrdinalTypeArray1DHost tagView(
"tag view", cardinality*tagSize);
462 const int vertexDim = 0, edgeDim = 1, faceDim = 2;
464 if (defineVertexFunctions) {
467 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
469 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerVertex; functionOrdinal++)
471 tagView(tagNumber*tagSize+0) = vertexDim;
472 tagView(tagNumber*tagSize+1) = vertexOrdinal;
473 tagView(tagNumber*tagSize+2) = functionOrdinal;
474 tagView(tagNumber*tagSize+3) = numFunctionsPerVertex;
478 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
480 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
482 tagView(tagNumber*tagSize+0) = edgeDim;
483 tagView(tagNumber*tagSize+1) = edgeOrdinal;
484 tagView(tagNumber*tagSize+2) = functionOrdinal;
485 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge;
489 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
491 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
493 tagView(tagNumber*tagSize+0) = faceDim;
494 tagView(tagNumber*tagSize+1) = faceOrdinal;
495 tagView(tagNumber*tagSize+2) = functionOrdinal;
496 tagView(tagNumber*tagSize+3) = numFunctionsPerFace;
502 for (ordinal_type i=0;i<cardinality;++i) {
503 tagView(i*tagSize+0) = faceDim;
504 tagView(i*tagSize+1) = 0;
505 tagView(i*tagSize+2) = i;
506 tagView(i*tagSize+3) = cardinality;
528 return "Intrepid2_IntegratedLegendreBasis_HGRAD_TRI";
560 virtual void getValues( OutputViewType outputValues,
const PointViewType inputPoints,
561 const EOperator operatorType = OPERATOR_VALUE )
const override 563 auto numPoints = inputPoints.extent_int(0);
567 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
569 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
570 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
571 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
572 const int teamSize = 1;
576 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
577 Kokkos::parallel_for( policy , functor,
"Hierarchical_HGRAD_TRI_Functor");
590 if(subCellDim == 1) {
591 return Teuchos::rcp(
new 595 INTREPID2_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Input parameters out of bounds");
604 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
606 return Teuchos::rcp(
new HostBasisType(polyOrder_, pointType_) );
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual bool requireOrientation() const override
True if orientation is required.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
const char * getName() const override
Returns basis name.
IntegratedLegendreBasis_HGRAD_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
EFunctionSpace functionSpace_
The function space in which the basis is defined.
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::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
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. This method is useful for sizing scratch storage in hierarchically parallel kernels. Whereas get_dimension_scalar() returns 1 for POD types, this returns 0 for POD types.
EBasis basisType_
Type of the basis.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
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...
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
Functor for computing values for the IntegratedLegendreBasis_HGRAD_TRI class.
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line...
EPointType
Enumeration of types of point distributions in Intrepid.
ordinal_type getDegree() const
Returns the degree of the basis.
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.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
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.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...