49 #ifndef Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h
50 #define Intrepid2_IntegratedLegendreBasis_HGRAD_TRI_h
52 #include <Kokkos_DynRankView.hpp>
54 #include <Intrepid2_config.h>
68 template<
class DeviceType,
class OutputScalar,
class PointScalar,
69 class OutputFieldType,
class InputPointsType>
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>>;
77 using TeamPolicy = Kokkos::TeamPolicy<ExecutionSpace>;
78 using TeamMember =
typename TeamPolicy::member_type;
82 OutputFieldType output_;
83 InputPointsType inputPoints_;
86 bool defineVertexFunctions_;
87 int numFields_, numPoints_;
89 size_t fad_size_output_;
91 static const int numVertices = 3;
92 static const int numEdges = 3;
93 const int edge_start_[numEdges] = {0,1,0};
94 const int edge_end_[numEdges] = {1,2,2};
97 int polyOrder,
bool defineVertexFunctions)
98 : opType_(opType), output_(output), inputPoints_(inputPoints),
99 polyOrder_(polyOrder), defineVertexFunctions_(defineVertexFunctions),
100 fad_size_output_(getScalarDimensionForView(output))
102 numFields_ = output.extent_int(0);
103 numPoints_ = output.extent_int(1);
104 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
105 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != (polyOrder_+1)*(polyOrder_+2)/2, std::invalid_argument,
"output field size does not match basis cardinality");
108 KOKKOS_INLINE_FUNCTION
109 void operator()(
const TeamMember & teamMember )
const
111 auto pointOrdinal = teamMember.league_rank();
112 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
113 if (fad_size_output_ > 0) {
114 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
115 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
116 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
117 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1, fad_size_output_);
120 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
121 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
122 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
123 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_ + 1);
126 const auto & x = inputPoints_(pointOrdinal,0);
127 const auto & y = inputPoints_(pointOrdinal,1);
130 const PointScalar lambda[3] = {1. - x - y, x, y};
131 const PointScalar lambda_dx[3] = {-1., 1., 0.};
132 const PointScalar lambda_dy[3] = {-1., 0., 1.};
134 const int num1DEdgeFunctions = polyOrder_ - 1;
141 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
143 output_(vertexOrdinal,pointOrdinal) = lambda[vertexOrdinal];
145 if (!defineVertexFunctions_)
149 output_(0,pointOrdinal) = 1.0;
153 int fieldOrdinalOffset = 3;
154 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
156 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
157 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
159 Polynomials::shiftedScaledIntegratedLegendreValues(edge_field_values_at_point, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
160 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
163 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal) = edge_field_values_at_point(edgeFunctionOrdinal+2);
165 fieldOrdinalOffset += num1DEdgeFunctions;
171 const double jacobiScaling = 1.0;
173 const int max_ij_sum = polyOrder_;
176 const int min_ij_sum = min_i + min_j;
177 for (
int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
179 for (
int i=min_i; i<=ij_sum-min_j; i++)
181 const int j = ij_sum - i;
182 const int edgeBasisOrdinal = i+numVertices-2;
183 const auto & edgeValue = output_(edgeBasisOrdinal,pointOrdinal);
184 const double alpha = i*2.0;
186 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-2, lambda[2], jacobiScaling);
187 const auto & jacobiValue = jacobi_values_at_point(j);
188 output_(fieldOrdinalOffset,pointOrdinal) = edgeValue * jacobiValue;
189 fieldOrdinalOffset++;
199 if (defineVertexFunctions_)
203 output_(0,pointOrdinal,0) = -1.0;
204 output_(0,pointOrdinal,1) = -1.0;
210 output_(0,pointOrdinal,0) = 0.0;
211 output_(0,pointOrdinal,1) = 0.0;
214 output_(1,pointOrdinal,0) = 1.0;
215 output_(1,pointOrdinal,1) = 0.0;
217 output_(2,pointOrdinal,0) = 0.0;
218 output_(2,pointOrdinal,1) = 1.0;
221 int fieldOrdinalOffset = 3;
233 auto & P_i_minus_1 = edge_field_values_at_point;
234 auto & L_i_dt = jacobi_values_at_point;
235 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
237 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
238 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
240 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
241 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
242 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
243 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
245 Polynomials::shiftedScaledLegendreValues (P_i_minus_1, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
246 Polynomials::shiftedScaledIntegratedLegendreValues_dt(L_i_dt, polyOrder_, PointScalar(s1), PointScalar(s0+s1));
247 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
250 const int i = edgeFunctionOrdinal+2;
251 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = P_i_minus_1(i-1) * s1_dx + L_i_dt(i) * (s1_dx + s0_dx);
252 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = P_i_minus_1(i-1) * s1_dy + L_i_dt(i) * (s1_dy + s0_dy);
254 fieldOrdinalOffset += num1DEdgeFunctions;
275 auto & P_2i_j_minus_1 = edge_field_values_at_point;
276 auto & L_2i_j_dt = jacobi_values_at_point;
277 auto & L_i = other_values_at_point;
278 auto & L_2i_j = other_values2_at_point;
281 const double jacobiScaling = 1.0;
283 const int max_ij_sum = polyOrder_;
286 const int min_ij_sum = min_i + min_j;
287 for (
int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
289 for (
int i=min_i; i<=ij_sum-min_j; i++)
291 const int j = ij_sum - i;
293 const int edgeBasisOrdinal = i+numVertices-2;
294 const auto & grad_L_i_dx = output_(edgeBasisOrdinal,pointOrdinal,0);
295 const auto & grad_L_i_dy = output_(edgeBasisOrdinal,pointOrdinal,1);
297 const double alpha = i*2.0;
299 Polynomials::shiftedScaledIntegratedLegendreValues (L_i, polyOrder_, lambda[1], lambda[0]+lambda[1]);
300 Polynomials::shiftedScaledIntegratedJacobiValues_dt(L_2i_j_dt, alpha, polyOrder_, lambda[2], jacobiScaling);
301 Polynomials::shiftedScaledIntegratedJacobiValues ( L_2i_j, alpha, polyOrder_, lambda[2], jacobiScaling);
302 Polynomials::shiftedScaledJacobiValues(P_2i_j_minus_1, alpha, polyOrder_-1, lambda[2], jacobiScaling);
304 const auto & s0_dx = lambda_dx[0];
305 const auto & s0_dy = lambda_dy[0];
306 const auto & s1_dx = lambda_dx[1];
307 const auto & s1_dy = lambda_dy[1];
308 const auto & s2_dx = lambda_dx[2];
309 const auto & s2_dy = lambda_dy[2];
311 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));
312 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));
314 output_(fieldOrdinalOffset,pointOrdinal,0) = basisValue_dx;
315 output_(fieldOrdinalOffset,pointOrdinal,1) = basisValue_dy;
316 fieldOrdinalOffset++;
331 INTREPID2_TEST_FOR_ABORT(
true,
332 ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
335 device_assert(
false);
342 size_t team_shmem_size (
int team_size)
const
345 size_t shmem_size = 0;
346 if (fad_size_output_ > 0)
347 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1, fad_size_output_);
349 shmem_size += 4 * OutputScratchView::shmem_size(polyOrder_ + 1);
371 template<
typename DeviceType,
372 typename OutputScalar = double,
373 typename PointScalar = double,
374 bool defineVertexFunctions =
true>
376 :
public Basis<DeviceType,OutputScalar,PointScalar>
393 EPointType pointType_;
407 polyOrder_(polyOrder),
408 pointType_(pointType)
410 INTREPID2_TEST_FOR_EXCEPTION(pointType!=POINTTYPE_DEFAULT,std::invalid_argument,
"PointType not supported");
414 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
419 const int degreeLength = 1;
423 int fieldOrdinalOffset = 0;
426 const int numFunctionsPerVertex = 1;
427 const int numVertexFunctions = numVertices * numFunctionsPerVertex;
428 for (
int i=0; i<numVertexFunctions; i++)
435 if (!defineVertexFunctions)
440 fieldOrdinalOffset += numVertexFunctions;
443 const int numFunctionsPerEdge = polyOrder - 1;
445 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
447 for (
int i=0; i<numFunctionsPerEdge; i++)
452 fieldOrdinalOffset += numFunctionsPerEdge;
456 const int max_ij_sum = polyOrder;
459 const int min_ij_sum = min_i + min_j;
460 for (
int ij_sum = min_ij_sum; ij_sum <= max_ij_sum; ij_sum++)
462 for (
int i=min_i; i<=ij_sum-min_j; i++)
464 const int j = ij_sum - i;
467 fieldOrdinalOffset++;
470 const int numFaces = 1;
471 const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
472 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->
basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
479 const ordinal_type tagSize = 4;
480 const ordinal_type posScDim = 0;
481 const ordinal_type posScOrd = 1;
482 const ordinal_type posDfOrd = 2;
485 const int vertexDim = 0, edgeDim = 1, faceDim = 2;
487 if (defineVertexFunctions) {
490 for (
int vertexOrdinal=0; vertexOrdinal<numVertices; vertexOrdinal++)
492 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerVertex; functionOrdinal++)
494 tagView(tagNumber*tagSize+0) = vertexDim;
495 tagView(tagNumber*tagSize+1) = vertexOrdinal;
496 tagView(tagNumber*tagSize+2) = functionOrdinal;
497 tagView(tagNumber*tagSize+3) = numFunctionsPerVertex;
501 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
503 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
505 tagView(tagNumber*tagSize+0) = edgeDim;
506 tagView(tagNumber*tagSize+1) = edgeOrdinal;
507 tagView(tagNumber*tagSize+2) = functionOrdinal;
508 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge;
512 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
514 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
516 tagView(tagNumber*tagSize+0) = faceDim;
517 tagView(tagNumber*tagSize+1) = faceOrdinal;
518 tagView(tagNumber*tagSize+2) = functionOrdinal;
519 tagView(tagNumber*tagSize+3) = numFunctionsPerFace;
525 for (ordinal_type i=0;i<cardinality;++i) {
526 tagView(i*tagSize+0) = faceDim;
527 tagView(i*tagSize+1) = 0;
528 tagView(i*tagSize+2) = i;
529 tagView(i*tagSize+3) = cardinality;
551 return "Intrepid2_IntegratedLegendreBasis_HGRAD_TRI";
584 const EOperator operatorType = OPERATOR_VALUE )
const override
586 auto numPoints = inputPoints.extent_int(0);
590 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_, defineVertexFunctions);
592 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
593 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
594 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
595 const int teamSize = 1;
597 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
598 Kokkos::parallel_for(
"Hierarchical_HGRAD_TRI_Functor", policy, functor);
609 BasisPtr<DeviceType,OutputScalar,PointScalar>
611 if(subCellDim == 1) {
612 return Teuchos::rcp(
new
616 INTREPID2_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Input parameters out of bounds");
623 virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>
625 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
627 return Teuchos::rcp(
new HostBasisType(polyOrder_, pointType_) );
H(grad) basis on the line based on integrated Legendre polynomials.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
virtual bool requireOrientation() const override
True if orientation is required.
virtual void getValues(const ExecutionSpace &, OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
const char * getName() const override
Returns basis name.
EBasis basisType_
Type of the basis.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
IntegratedLegendreBasis_HGRAD_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
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.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
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.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
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.
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: e...
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...
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.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
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 fieldOrdinalH1PolynomialDegree_
H^1 polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
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< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
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.