Intrepid2
Intrepid2_SerendipityBasis.hpp
1//
2// Intrepid2_SerendipityBasis.hpp
3// intrepid2
4//
5// Created by Roberts, Nathan V on 1/4/22.
6//
7
8#ifndef Intrepid2_SerendipityBasis_h
9#define Intrepid2_SerendipityBasis_h
10
12
13namespace Intrepid2
14{
15
19template<typename BasisBaseClass = void>
21:
22public BasisBaseClass
23{
24public:
25 using BasisBase = BasisBaseClass;
26 using BasisPtr = Teuchos::RCP<BasisBase>;
27
28protected:
29 BasisPtr fullBasis_;
30
31 std::string name_; // name of the basis
32
33 int numTensorialExtrusions_; // relative to cell topo returned by getBaseCellTopology().
34public:
35 using DeviceType = typename BasisBase::DeviceType;
36 using ExecutionSpace = typename BasisBase::ExecutionSpace;
37 using OutputValueType = typename BasisBase::OutputValueType;
38 using PointValueType = typename BasisBase::PointValueType;
39
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;
45protected:
46 OrdinalTypeArray1D ordinalMap_; // our field ordinal --> full basis field ordinal
47public:
51 SerendipityBasis(BasisPtr fullBasis)
52 :
53 fullBasis_(fullBasis)
54 {
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(); // BASIS_FEM_HIERARCHICAL
57
58 this->functionSpace_ = fullBasis_->getFunctionSpace();
59 this->basisDegree_ = fullBasis_->getDegree();
60
61 {
62 std::ostringstream basisName;
63 basisName << "Serendipity(" << fullBasis_->getName() << ")";
64 name_ = basisName.str();
65 }
66
67 using std::vector;
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++)
74 {
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)
79 {
80 if (p > 1)
81 {
82 // superlinear; contributes to superlinearDegree
83 superlinearDegree += p;
84 }
85 }
86 if (superlinearDegree <= basisH1Degree)
87 {
88 // satisfies serendipity criterion
89 fullBasisOrdinals.push_back(i);
90 polynomialDegreeOfField.push_back(fieldDegree);
91 polynomialH1DegreeOfField.push_back(fieldH1Degree);
92 }
93 }
94 this->basisCardinality_ = fullBasisOrdinals.size();
95 ordinalMap_ = OrdinalTypeArray1D("serendipity ordinal map",fullBasisOrdinals.size());
96
97 auto ordinalMapHost = Kokkos::create_mirror_view(ordinalMap_);
98 const ordinal_type fullBasisCardinality = fullBasisOrdinals.size();
99 for (ordinal_type i=0; i<fullBasisCardinality; i++)
100 {
101 ordinalMapHost(i) = fullBasisOrdinals[i];
102 }
103 Kokkos::deep_copy(ordinalMap_, ordinalMapHost);
104
105 // set cell topology
106 this->basisCellTopology_ = fullBasis_->getBaseCellTopology();
107 this->numTensorialExtrusions_ = fullBasis_->getNumTensorialExtrusions();
108 const ordinal_type spaceDim = this->basisCellTopology_.getDimension() + numTensorialExtrusions_;
109
110 this->basisCoordinates_ = COORDINATES_CARTESIAN;
111
112 // fill in degree lookup:
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);
116
117 for (ordinal_type fieldOrdinal1 = 0; fieldOrdinal1 < this->basisCardinality_; fieldOrdinal1++)
118 {
119 for (int d=0; d<degreeSize; d++)
120 {
121 this->fieldOrdinalPolynomialDegree_ (fieldOrdinal1,d) = polynomialDegreeOfField [fieldOrdinal1][d];
122 this->fieldOrdinalH1PolynomialDegree_(fieldOrdinal1,d) = polynomialH1DegreeOfField[fieldOrdinal1][d];
123 }
124 }
125
126 // build tag view
127 const auto & cardinality = this->basisCardinality_;
128
129 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
130 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
131 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
132 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
133 const ordinal_type posDfCnt = 3; // position in the tag, counting from 0, of DoF count for the subcell
134
135 OrdinalTypeArray1DHost tagView("tag view", cardinality*tagSize);
136 auto cellTopo = CellTopology::cellTopology(this->basisCellTopology_, numTensorialExtrusions_);
137 const OrdinalTypeArray2DHost fullBasisOrdinalToTag = fullBasis->getAllDofTags();
138
139 using std::map;
140 vector<vector<ordinal_type>> subcellDofCount(spaceDim+1); // outer: dimension of subcell; inner: subcell ordinal. Entry: number of dofs.
141 vector<vector<ordinal_type>> subcellDofOrdinal(spaceDim+1); // outer: dimension of subcell; inner: subcell ordinal. Entry: ordinal relative to subcell.
142 for (ordinal_type d=0; d<=spaceDim; d++)
143 {
144 const ordinal_type numSubcells = cellTopo->getSubcellCount(d);
145 subcellDofCount[d] = vector<ordinal_type>(numSubcells,0);
146 subcellDofOrdinal[d] = vector<ordinal_type>(numSubcells,0);
147 }
148 for (ordinal_type fieldOrdinal=0; fieldOrdinal<fullBasisCardinality; fieldOrdinal++)
149 {
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]++;
154 }
155 for (ordinal_type fieldOrdinal=0; fieldOrdinal<fullBasisCardinality; fieldOrdinal++)
156 {
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]++;
162
163 tagView(tagSize*fieldOrdinal + posScDim) = subcellDim; // subcellDim
164 tagView(tagSize*fieldOrdinal + posScOrd) = subcellOrd; // subcell ordinal
165 tagView(tagSize*fieldOrdinal + posDfOrd) = subcellDfOrd; // ordinal of the specified DoF relative to the subcell
166 tagView(tagSize*fieldOrdinal + posDfCnt) = subcellDfCnt; // total number of DoFs associated with the subcell
167 }
168 // // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
169 // // tags are constructed on host
170 this->setOrdinalTagData(this->tagToOrdinal_,
171 this->ordinalToTag_,
172 tagView,
173 this->basisCardinality_,
174 tagSize,
175 posScDim,
176 posScOrd,
177 posDfOrd);
178 }
179
180 virtual int getNumTensorialExtrusions() const override
181 {
182 return numTensorialExtrusions_;
183 }
184
189 virtual BasisValues<OutputValueType,DeviceType> allocateBasisValues( TensorPoints<PointValueType,DeviceType> points, const EOperator operatorType = OPERATOR_VALUE) const override
190 {
191 auto basisValues = fullBasis_->allocateBasisValues(points,operatorType);
192 basisValues.setOrdinalFilter(this->ordinalMap_);
193 return basisValues;
194 }
195
196 // since the getValues() below only overrides the FEM variant, we specify that
197 // we use the base class's getValues(), which implements the FVD variant by throwing an exception.
198 // (It's an error to use the FVD variant on this basis.)
199 using BasisBase::getValues;
200
205 virtual
206 const char*
207 getName() const override {
208 return name_.c_str();
209 }
210
222 virtual
223 void
226 const EOperator operatorType = OPERATOR_VALUE ) const override
227 {
228 fullBasis_->getValues(outputValues, inputPoints, operatorType);
229 outputValues.setOrdinalFilter(this->ordinalMap_);
230 }
231
250 virtual
251 void getValues( OutputViewType outputValues, const PointViewType inputPoints,
252 const EOperator operatorType = OPERATOR_VALUE ) const override
253 {
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).");
256 }
257
263 getHostBasis() const override {
265 using HostBasis = SerendipityBasis<HostBasisBase>;
266 auto hostBasis = Teuchos::rcp(new HostBasis(fullBasis_->getHostBasis()));
267 return hostBasis;
268 }
269
274 BasisPtr getUnderlyingBasis() const
275 {
276 return fullBasis_;
277 }
278
283 OrdinalTypeArray1D ordinalMap() const
284 {
285 return ordinalMap_;
286 }
287}; // SerendipityBasis
288
289} // namespace Intrepid2
290#endif /* Intrepid2_SerendipityBasis_h */
BasisPtr< typename Kokkos::HostSpace::device_type, OutputType, PointType > HostBasisPtr
Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type),...
Implements arbitrary-dimensional extrusion of a base shards::CellTopology.
The data containers in Intrepid2 that support sum factorization and other reduced-data optimizations ...
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
static CellTopoPtr cellTopology(const shards::CellTopology &shardsCellTopo, ordinal_type tensorialDegree=0)
static accessor that returns a CellTopoPtr; these are lazily constructed and cached.
Serendipity Basis, defined as the sub-basis of a provided basis, consisting of basis elements for whi...
SerendipityBasis(BasisPtr fullBasis)
Constructor.
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...
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
BasisPtr getUnderlyingBasis() const
Returns a pointer to the underlying full basis.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
OrdinalTypeArray1D ordinalMap() const
Returns the ordinal map from the Serendipity basis ordinal to the ordinal in the underlying full basi...
virtual const char * getName() const override
Returns basis name.
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...
View-like interface to tensor points; point components are stored separately; the appropriate coordin...