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_,