8 #ifndef Intrepid2_SerendipityBasis_h
9 #define Intrepid2_SerendipityBasis_h
19 template<
typename BasisBaseClass =
void>
25 using BasisBase = BasisBaseClass;
26 using BasisPtr = Teuchos::RCP<BasisBase>;
33 int numTensorialExtrusions_;
35 using DeviceType =
typename BasisBase::DeviceType;
36 using ExecutionSpace =
typename BasisBase::ExecutionSpace;
37 using OutputValueType =
typename BasisBase::OutputValueType;
38 using PointValueType =
typename BasisBase::PointValueType;
40 using OrdinalTypeArray1D =
typename BasisBase::OrdinalTypeArray1D;
41 using OrdinalTypeArray1DHost =
typename BasisBase::OrdinalTypeArray1DHost;
42 using OrdinalTypeArray2DHost =
typename BasisBase::OrdinalTypeArray2DHost;
43 using OutputViewType =
typename BasisBase::OutputViewType;
44 using PointViewType =
typename BasisBase::PointViewType;
46 OrdinalTypeArray1D ordinalMap_;
55 INTREPID2_TEST_FOR_EXCEPTION(fullBasis_->getBasisType() != BASIS_FEM_HIERARCHICAL, std::invalid_argument,
"SerendipityBasis only supports full bases whose type is BASIS_FEM_HIERARCHICAL");
56 this->basisType_ = fullBasis_->getBasisType();
58 this->functionSpace_ = fullBasis_->getFunctionSpace();
59 this->basisDegree_ = fullBasis_->getDegree();
62 std::ostringstream basisName;
63 basisName <<
"Serendipity(" << fullBasis_->getName() <<
")";
64 name_ = basisName.str();
68 vector<ordinal_type> fullBasisOrdinals;
69 vector<vector<ordinal_type>> polynomialDegreeOfField;
70 vector<vector<ordinal_type>> polynomialH1DegreeOfField;
71 ordinal_type fullCardinality = fullBasis_->getCardinality();
72 ordinal_type basisH1Degree = (this->functionSpace_ == FUNCTION_SPACE_HVOL) ? this->basisDegree_ + 1 : this->basisDegree_;
73 for (ordinal_type i=0; i<fullCardinality; i++)
75 vector<ordinal_type> fieldDegree = fullBasis_->getPolynomialDegreeOfFieldAsVector(i);
76 vector<ordinal_type> fieldH1Degree = fullBasis_->getH1PolynomialDegreeOfFieldAsVector(i);
77 ordinal_type superlinearDegree = 0;
78 for (
const ordinal_type & p : fieldH1Degree)
83 superlinearDegree += p;
86 if (superlinearDegree <= basisH1Degree)
89 fullBasisOrdinals.push_back(i);
90 polynomialDegreeOfField.push_back(fieldDegree);
91 polynomialH1DegreeOfField.push_back(fieldH1Degree);
94 this->basisCardinality_ = fullBasisOrdinals.size();
95 ordinalMap_ = OrdinalTypeArray1D(
"serendipity ordinal map",fullBasisOrdinals.size());
97 auto ordinalMapHost = Kokkos::create_mirror_view(ordinalMap_);
98 const ordinal_type fullBasisCardinality = fullBasisOrdinals.size();
99 for (ordinal_type i=0; i<fullBasisCardinality; i++)
101 ordinalMapHost(i) = fullBasisOrdinals[i];
103 Kokkos::deep_copy(ordinalMap_, ordinalMapHost);
106 this->basisCellTopology_ = fullBasis_->getBaseCellTopology();
107 this->numTensorialExtrusions_ = fullBasis_->getNumTensorialExtrusions();
108 const ordinal_type spaceDim = this->basisCellTopology_.getDimension() + numTensorialExtrusions_;
110 this->basisCoordinates_ = COORDINATES_CARTESIAN;
113 int degreeSize = fullBasis_->getPolynomialDegreeLength();
114 this->fieldOrdinalPolynomialDegree_ = OrdinalTypeArray2DHost(
"TensorBasis - field ordinal polynomial degree", this->basisCardinality_, degreeSize);
115 this->fieldOrdinalH1PolynomialDegree_ = OrdinalTypeArray2DHost(
"TensorBasis - field ordinal H^1 degree", this->basisCardinality_, degreeSize);
117 for (ordinal_type fieldOrdinal1 = 0; fieldOrdinal1 < this->basisCardinality_; fieldOrdinal1++)
119 for (
int d=0; d<degreeSize; d++)
121 this->fieldOrdinalPolynomialDegree_ (fieldOrdinal1,d) = polynomialDegreeOfField [fieldOrdinal1][d];
122 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinal1,d) = polynomialH1DegreeOfField[fieldOrdinal1][d];
127 const auto & cardinality = this->basisCardinality_;
129 const ordinal_type tagSize = 4;
130 const ordinal_type posScDim = 0;
131 const ordinal_type posScOrd = 1;
132 const ordinal_type posDfOrd = 2;
133 const ordinal_type posDfCnt = 3;
135 OrdinalTypeArray1DHost tagView(
"tag view", cardinality*tagSize);
137 const OrdinalTypeArray2DHost fullBasisOrdinalToTag = fullBasis->getAllDofTags();
140 vector<vector<ordinal_type>> subcellDofCount(spaceDim+1);
141 vector<vector<ordinal_type>> subcellDofOrdinal(spaceDim+1);
142 for (ordinal_type d=0; d<=spaceDim; d++)
144 const ordinal_type numSubcells = cellTopo->getSubcellCount(d);
145 subcellDofCount[d] = vector<ordinal_type>(numSubcells,0);
146 subcellDofOrdinal[d] = vector<ordinal_type>(numSubcells,0);
148 for (ordinal_type fieldOrdinal=0; fieldOrdinal<fullBasisCardinality; fieldOrdinal++)
150 const ordinal_type fullFieldOrdinal = fullBasisOrdinals[fieldOrdinal];
151 const ordinal_type subcellDim = fullBasisOrdinalToTag(fullFieldOrdinal,posScDim);
152 const ordinal_type subcellOrd = fullBasisOrdinalToTag(fullFieldOrdinal,posScOrd);
153 subcellDofCount[subcellDim][subcellOrd]++;
155 for (ordinal_type fieldOrdinal=0; fieldOrdinal<fullBasisCardinality; fieldOrdinal++)
157 const ordinal_type fullFieldOrdinal = fullBasisOrdinals[fieldOrdinal];
158 const ordinal_type subcellDim = fullBasisOrdinalToTag(fullFieldOrdinal,posScDim);
159 const ordinal_type subcellOrd = fullBasisOrdinalToTag(fullFieldOrdinal,posScOrd);
160 const ordinal_type subcellDfCnt = subcellDofCount[subcellDim][subcellOrd];
161 const ordinal_type subcellDfOrd = subcellDofOrdinal[subcellDim][subcellOrd]++;
163 tagView(tagSize*fieldOrdinal + posScDim) = subcellDim;
164 tagView(tagSize*fieldOrdinal + posScOrd) = subcellOrd;
165 tagView(tagSize*fieldOrdinal + posDfOrd) = subcellDfOrd;
166 tagView(tagSize*fieldOrdinal + posDfCnt) = subcellDfCnt;
170 this->setOrdinalTagData(this->tagToOrdinal_,
173 this->basisCardinality_,
180 virtual int getNumTensorialExtrusions()
const override
182 return numTensorialExtrusions_;
191 auto basisValues = fullBasis_->allocateBasisValues(points,operatorType);
192 basisValues.setOrdinalFilter(this->ordinalMap_);
199 using BasisBase::getValues;
208 return name_.c_str();
226 const EOperator operatorType = OPERATOR_VALUE )
const override
228 fullBasis_->getValues(outputValues, inputPoints, operatorType);
229 outputValues.setOrdinalFilter(this->ordinalMap_);
251 void getValues( OutputViewType outputValues,
const PointViewType inputPoints,
252 const EOperator operatorType = OPERATOR_VALUE )
const override
254 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
255 ">>> ERROR (Basis::getValues): this method is not supported by SerendipityBasis (use the getValues() method that accepts a BasisValues object instead).");
262 virtual HostBasisPtr<OutputValueType, PointValueType>
266 auto hostBasis = Teuchos::rcp(
new HostBasis(fullBasis_->getHostBasis()));
Implements arbitrary-dimensional extrusion of a base shards::CellTopology.
Serendipity Basis, defined as the sub-basis of a provided basis, consisting of basis elements for whi...
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
virtual const char * getName() const override
Returns basis name.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
virtual void getValues(BasisValues< OutputValueType, DeviceType > outputValues, const TensorPoints< PointValueType, DeviceType > inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell, using point and output value containers that allow pre...
The data containers in Intrepid2 that support sum factorization and other reduced-data optimizations ...
virtual BasisValues< OutputValueType, DeviceType > allocateBasisValues(TensorPoints< PointValueType, DeviceType > points, const EOperator operatorType=OPERATOR_VALUE) const override
Allocate BasisValues container suitable for passing to the getValues() variant that takes a TensorPoi...
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
BasisPtr getUnderlyingBasis() const
Returns a pointer to the underlying full basis.
OrdinalTypeArray1D ordinalMap() const
Returns the ordinal map from the Serendipity basis ordinal to the ordinal in the underlying full basi...
static CellTopoPtr cellTopology(const shards::CellTopology &shardsCellTopo, ordinal_type tensorialDegree=0)
static accessor that returns a CellTopoPtr; these are lazily constructed and cached.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
SerendipityBasis(BasisPtr fullBasis)
Constructor.