49 #ifndef Intrepid2_HierarchicalBasis_HCURL_TRI_h
50 #define Intrepid2_HierarchicalBasis_HCURL_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 int numFields_, numPoints_;
88 size_t fad_size_output_;
90 static const int numVertices = 3;
91 static const int numEdges = 3;
92 static const int numFaceFamilies = 2;
93 const int edge_start_[numEdges] = {0,1,0};
94 const int edge_end_[numEdges] = {1,2,2};
95 const int face_family_start_[numFaceFamilies] = {0,1};
96 const int face_family_middle_[numFaceFamilies] = {1,2};
97 const int face_family_end_ [numFaceFamilies] = {2,0};
100 : opType_(opType), output_(output), inputPoints_(inputPoints),
101 polyOrder_(polyOrder),
102 fad_size_output_(getScalarDimensionForView(output))
104 numFields_ = output.extent_int(0);
105 numPoints_ = output.extent_int(1);
106 const int expectedCardinality = 3 * polyOrder_ + polyOrder_ * (polyOrder-1);
108 INTREPID2_TEST_FOR_EXCEPTION(numPoints_ != inputPoints.extent_int(0), std::invalid_argument,
"point counts need to match!");
109 INTREPID2_TEST_FOR_EXCEPTION(numFields_ != expectedCardinality, std::invalid_argument,
"output field size does not match basis cardinality");
112 KOKKOS_INLINE_FUNCTION
113 void operator()(
const TeamMember & teamMember )
const
115 auto pointOrdinal = teamMember.league_rank();
116 OutputScratchView edge_field_values_at_point, jacobi_values_at_point, other_values_at_point, other_values2_at_point;
117 if (fad_size_output_ > 0) {
118 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
119 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
120 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
121 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_, fad_size_output_);
124 edge_field_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
125 jacobi_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
126 other_values_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
127 other_values2_at_point = OutputScratchView(teamMember.team_shmem(), polyOrder_);
130 const auto & x = inputPoints_(pointOrdinal,0);
131 const auto & y = inputPoints_(pointOrdinal,1);
134 const PointScalar lambda[3] = {1. - x - y, x, y};
135 const PointScalar lambda_dx[3] = {-1., 1., 0.};
136 const PointScalar lambda_dy[3] = {-1., 0., 1.};
138 const int num1DEdgeFunctions = polyOrder_;
145 int fieldOrdinalOffset = 0;
146 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
148 const auto & s0 = lambda [edge_start_[edgeOrdinal]];
149 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
150 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
152 const auto & s1 = lambda [ edge_end_[edgeOrdinal]];
153 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
154 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
156 Polynomials::shiftedScaledLegendreValues(edge_field_values_at_point, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
157 for (
int edgeFunctionOrdinal=0; edgeFunctionOrdinal<num1DEdgeFunctions; edgeFunctionOrdinal++)
159 const auto & legendreValue = edge_field_values_at_point(edgeFunctionOrdinal);
160 const PointScalar xWeight = s0 * s1_dx - s1 * s0_dx;
161 const PointScalar yWeight = s0 * s1_dy - s1 * s0_dy;
162 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,0) = legendreValue * xWeight;
163 output_(edgeFunctionOrdinal+fieldOrdinalOffset,pointOrdinal,1) = legendreValue * yWeight;
165 fieldOrdinalOffset += num1DEdgeFunctions;
171 const double jacobiScaling = 1.0;
173 const int max_ij_sum = polyOrder_ - 1;
176 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
177 for (
int familyOrdinal=1; familyOrdinal<=2; familyOrdinal++)
179 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
180 const auto &s2 = lambda[ face_family_end_[familyOrdinal-1]];
181 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
183 for (
int i=0; i<ij_sum; i++)
185 const int j = ij_sum - i;
187 const int edgeBasisOrdinal = i + (familyOrdinal-1)*num1DEdgeFunctions;
188 const auto & edgeValue_x = output_(edgeBasisOrdinal,pointOrdinal,0);
189 const auto & edgeValue_y = output_(edgeBasisOrdinal,pointOrdinal,1);
190 const double alpha = i*2.0 + 1;
192 Polynomials::shiftedScaledIntegratedJacobiValues(jacobi_values_at_point, alpha, polyOrder_-1, s2, jacobiScaling);
193 const auto & jacobiValue = jacobi_values_at_point(j);
194 output_(fieldOrdinal,pointOrdinal,0) = edgeValue_x * jacobiValue;
195 output_(fieldOrdinal,pointOrdinal,1) = edgeValue_y * jacobiValue;
207 int fieldOrdinalOffset = 0;
214 auto & P_i = edge_field_values_at_point;
215 auto & L_2ip1_j = jacobi_values_at_point;
216 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
218 const auto & s0 = lambda[edge_start_[edgeOrdinal]];
219 const auto & s1 = lambda[ edge_end_[edgeOrdinal]];
221 const auto & s0_dx = lambda_dx[edge_start_[edgeOrdinal]];
222 const auto & s0_dy = lambda_dy[edge_start_[edgeOrdinal]];
223 const auto & s1_dx = lambda_dx[ edge_end_[edgeOrdinal]];
224 const auto & s1_dy = lambda_dy[ edge_end_[edgeOrdinal]];
226 const OutputScalar grad_s0_cross_grad_s1 = s0_dx * s1_dy - s1_dx * s0_dy;
228 Polynomials::shiftedScaledLegendreValues(P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
229 for (
int i=0; i<num1DEdgeFunctions; i++)
231 output_(i+fieldOrdinalOffset,pointOrdinal) = (i+2) * P_i(i) * grad_s0_cross_grad_s1;
233 fieldOrdinalOffset += num1DEdgeFunctions;
254 auto & P_2ip1_j = other_values_at_point;
257 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
258 for (
int familyOrdinal=1; familyOrdinal<=2; familyOrdinal++)
260 int fieldOrdinal = faceFieldOrdinalOffset + familyOrdinal - 1;
262 const auto &s0_index = face_family_start_ [familyOrdinal-1];
263 const auto &s1_index = face_family_middle_[familyOrdinal-1];
264 const auto &s2_index = face_family_end_ [familyOrdinal-1];
265 const auto &s0 = lambda[s0_index];
266 const auto &s1 = lambda[s1_index];
267 const auto &s2 = lambda[s2_index];
268 const double jacobiScaling = 1.0;
270 const auto & s0_dx = lambda_dx[s0_index];
271 const auto & s0_dy = lambda_dy[s0_index];
272 const auto & s1_dx = lambda_dx[s1_index];
273 const auto & s1_dy = lambda_dy[s1_index];
274 const auto & s2_dx = lambda_dx[s2_index];
275 const auto & s2_dy = lambda_dy[s2_index];
277 const OutputScalar grad_s0_cross_grad_s1 = s0_dx * s1_dy - s1_dx * s0_dy;
279 Polynomials::shiftedScaledLegendreValues (P_i, polyOrder_-1, PointScalar(s1), PointScalar(s0+s1));
283 const PointScalar xEdgeWeight = s0 * s1_dx - s1 * s0_dx;
284 const PointScalar yEdgeWeight = s0 * s1_dy - s1 * s0_dy;
285 OutputScalar grad_s2_cross_xy_edgeWeight = s2_dx * yEdgeWeight - xEdgeWeight * s2_dy;
287 const int max_ij_sum = polyOrder_ - 1;
288 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
290 for (
int i=0; i<ij_sum; i++)
292 const int j = ij_sum - i;
293 const OutputScalar edgeCurl = (i+2.) * P_i(i) * grad_s0_cross_grad_s1;
295 const double alpha = i*2.0 + 1;
297 Polynomials::shiftedScaledJacobiValues(P_2ip1_j, alpha, polyOrder_-1, PointScalar(s2), jacobiScaling);
298 Polynomials::shiftedScaledIntegratedJacobiValues(L_2ip1_j, alpha, polyOrder_-1, s2, jacobiScaling);
300 const PointScalar & edgeValue = P_i(i);
301 output_(fieldOrdinal,pointOrdinal) = L_2ip1_j(j) * edgeCurl + P_2ip1_j(j-1) * edgeValue * grad_s2_cross_xy_edgeWeight;
320 INTREPID2_TEST_FOR_ABORT(
true,
321 ">>> ERROR: (Intrepid2::Hierarchical_HCURL_TRI_Functor) Unsupported differential operator");
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);
358 template<
typename DeviceType,
359 typename OutputScalar = double,
360 typename PointScalar = double,
361 bool useCGBasis =
true>
363 :
public Basis<DeviceType,OutputScalar,PointScalar>
388 polyOrder_(polyOrder)
390 const int numEdgeFunctions = polyOrder * 3;
391 const int numFaceFunctions = polyOrder * (polyOrder-1);
394 this->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
399 const int degreeLength = 1;
403 int fieldOrdinalOffset = 0;
408 const int numFunctionsPerEdge = polyOrder;
410 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
412 for (
int i=0; i<numFunctionsPerEdge; i++)
417 fieldOrdinalOffset += numFunctionsPerEdge;
419 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != numEdgeFunctions, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
422 const int max_ij_sum = polyOrder-1;
423 const int faceFieldOrdinalOffset = fieldOrdinalOffset;
424 for (
int faceFamilyOrdinal=1; faceFamilyOrdinal<=2; faceFamilyOrdinal++)
427 int fieldOrdinal = faceFieldOrdinalOffset + faceFamilyOrdinal - 1;
428 for (
int ij_sum=1; ij_sum <= max_ij_sum; ij_sum++)
430 for (
int i=0; i<ij_sum; i++)
437 fieldOrdinalOffset = fieldOrdinal - 1;
440 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument,
"Internal error: basis enumeration is incorrect");
447 const ordinal_type tagSize = 4;
448 const ordinal_type posScDim = 0;
449 const ordinal_type posScOrd = 1;
450 const ordinal_type posDfOrd = 2;
453 const ordinal_type edgeDim = 1, faceDim = 2;
458 for (
int edgeOrdinal=0; edgeOrdinal<numEdges; edgeOrdinal++)
460 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerEdge; functionOrdinal++)
462 tagView(tagNumber*tagSize+0) = edgeDim;
463 tagView(tagNumber*tagSize+1) = edgeOrdinal;
464 tagView(tagNumber*tagSize+2) = functionOrdinal;
465 tagView(tagNumber*tagSize+3) = numFunctionsPerEdge;
469 const int numFunctionsPerFace = numFaceFunctions;
470 const int numFaces = 1;
471 for (
int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
473 for (
int functionOrdinal=0; functionOrdinal<numFunctionsPerFace; functionOrdinal++)
475 tagView(tagNumber*tagSize+0) = faceDim;
476 tagView(tagNumber*tagSize+1) = faceOrdinal;
477 tagView(tagNumber*tagSize+2) = functionOrdinal;
478 tagView(tagNumber*tagSize+3) = numFunctionsPerFace;
487 for (ordinal_type i=0;i<cardinality;++i) {
488 tagView(i*tagSize+0) = faceDim;
489 tagView(i*tagSize+1) = 0;
490 tagView(i*tagSize+2) = i;
491 tagView(i*tagSize+3) = cardinality;
500 this->basisCardinality_,
513 return "Intrepid2_HierarchicalBasis_HCURL_TRI";
546 const EOperator operatorType = OPERATOR_VALUE )
const override
548 auto numPoints = inputPoints.extent_int(0);
552 FunctorType functor(operatorType, outputValues, inputPoints, polyOrder_);
554 const int outputVectorSize = getVectorSizeForHierarchicalParallelism<OutputScalar>();
555 const int pointVectorSize = getVectorSizeForHierarchicalParallelism<PointScalar>();
556 const int vectorSize = std::max(outputVectorSize,pointVectorSize);
557 const int teamSize = 1;
559 auto policy = Kokkos::TeamPolicy<ExecutionSpace>(numPoints,teamSize,vectorSize);
560 Kokkos::parallel_for(
"Hierarchical_HCURL_TRI_Functor", policy, functor);
571 BasisPtr<DeviceType,OutputScalar,PointScalar>
573 if(subCellDim == 1) {
574 return Teuchos::rcp(
new
578 INTREPID2_TEST_FOR_EXCEPTION(
true,std::invalid_argument,
"Input parameters out of bounds");
585 virtual BasisPtr<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar>
587 using HostDeviceType =
typename Kokkos::HostSpace::device_type;
589 return Teuchos::rcp(
new HostBasisType(polyOrder_) );
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 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.
Functor for computing values for the HierarchicalBasis_HCURL_TRI class.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
HierarchicalBasis_HCURL_TRI(int polyOrder, const EPointType pointType=POINTTYPE_DEFAULT)
Constructor.
virtual bool requireOrientation() const override
True if orientation is required.
Free functions, callable from device code, that implement various polynomials useful in basis definit...
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
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.
For mathematical details of the construction, see:
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line...
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
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< 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. ...
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
BasisPtr< DeviceType, OutputScalar, PointScalar > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
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...
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.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
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.
Header file for the abstract base class Intrepid2::Basis.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.