1 #ifndef INTREPID_HGRAD_PYR_C1_FEMDEF_HPP 2 #define INTREPID_HGRAD_PYR_C1_FEMDEF_HPP 56 template<
class Scalar,
class ArrayScalar>
59 this -> basisCardinality_ = 5;
60 this -> basisDegree_ = 1;
61 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() );
62 this -> basisType_ = BASIS_FEM_DEFAULT;
63 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
64 this -> basisTagsAreSet_ =
false;
68 template<
class Scalar,
class ArrayScalar>
78 int tags[] = { 0, 0, 0, 1,
86 this -> ordinalToTag_,
88 this -> basisCardinality_,
97 template<
class Scalar,
class ArrayScalar>
99 const ArrayScalar & inputPoints,
103 #ifdef HAVE_INTREPID_DEBUG 104 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
107 this -> getBaseCellTopology(),
108 this -> getCardinality() );
112 int dim0 = inputPoints.dimension(0);
118 const Scalar eps = std::numeric_limits<Scalar>::epsilon( );
120 switch (operatorType) {
123 for (
int i0 = 0; i0 < dim0; i0++) {
124 x = inputPoints(i0, 0);
125 y = inputPoints(i0, 1);
126 z = inputPoints(i0, 2);
129 if(fabs(z-1.0) < eps) {
130 if(z <= 1.0) z = 1.0-eps;
135 Scalar zTerm = 0.25/(1.0 - z);
138 outputValues(0, i0) = (1.0 - x - z) * (1.0 - y - z) * zTerm;
139 outputValues(1, i0) = (1.0 + x - z) * (1.0 - y - z) * zTerm;
140 outputValues(2, i0) = (1.0 + x - z) * (1.0 + y - z) * zTerm;
141 outputValues(3, i0) = (1.0 - x - z) * (1.0 + y - z) * zTerm;
142 outputValues(4, i0) = z;
148 for (
int i0 = 0; i0 < dim0; i0++) {
150 x = inputPoints(i0, 0);
151 y = inputPoints(i0, 1);
152 z = inputPoints(i0, 2);
157 if(fabs(z-1.0) < eps) {
158 if(z <= 1.0) z = 1.0-eps;
163 Scalar zTerm = 0.25/(1.0 - z);
164 Scalar zTerm2 = 4.0 * zTerm * zTerm;
167 outputValues(0, i0, 0) = (y + z - 1.0) * zTerm;
168 outputValues(0, i0, 1) = (x + z - 1.0) * zTerm;
169 outputValues(0, i0, 2) = x * y * zTerm2 - 0.25;
171 outputValues(1, i0, 0) = (1.0 - y - z) * zTerm;
172 outputValues(1, i0, 1) = (z - x - 1.0) * zTerm;
173 outputValues(1, i0, 2) = - x*y * zTerm2 - 0.25;
175 outputValues(2, i0, 0) = (1.0 + y - z) * zTerm;
176 outputValues(2, i0, 1) = (1.0 + x - z) * zTerm;
177 outputValues(2, i0, 2) = x * y * zTerm2 - 0.25;
179 outputValues(3, i0, 0) = (z - y - 1.0) * zTerm;
180 outputValues(3, i0, 1) = (1.0 - x - z) * zTerm;
181 outputValues(3, i0, 2) = - x*y * zTerm2 - 0.25;
183 outputValues(4, i0, 0) = 0.0;
184 outputValues(4, i0, 1) = 0.0;
185 outputValues(4, i0, 2) = 1;
190 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
191 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
195 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
196 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
200 for (
int i0 = 0; i0 < dim0; i0++) {
201 x = inputPoints(i0,0);
202 y = inputPoints(i0,1);
203 z = inputPoints(i0,2);
207 if(fabs(z-1.0) < eps) {
208 if(z <= 1.0) z = 1.0-eps;
213 Scalar zTerm = 0.25/(1.0 - z);
214 Scalar zTerm2 = 4.0 * zTerm * zTerm;
215 Scalar zTerm3 = 8.0 * zTerm * zTerm2;
218 outputValues(0, i0, 0) = 0.0;
219 outputValues(0, i0, 1) = zTerm;
220 outputValues(0, i0, 2) = y*zTerm2;
221 outputValues(0, i0, 3) = 0.0;
222 outputValues(0, i0, 4) = x*zTerm2;
223 outputValues(0, i0, 5) = x*y*zTerm3;
225 outputValues(1, i0, 0) = 0.0;
226 outputValues(1, i0, 1) = -zTerm;
227 outputValues(1, i0, 2) = -y*zTerm2;
228 outputValues(1, i0, 3) = 0.0;
229 outputValues(1, i0, 4) = -x*zTerm2;
230 outputValues(1, i0, 5) = -x*y*zTerm3;
232 outputValues(2, i0, 0) = 0.0;
233 outputValues(2, i0, 1) = zTerm;
234 outputValues(2, i0, 2) = y*zTerm2;
235 outputValues(2, i0, 3) = 0.0;
236 outputValues(2, i0, 4) = x*zTerm2;
237 outputValues(2, i0, 5) = x*y*zTerm3;
239 outputValues(3, i0, 0) = 0.0;
240 outputValues(3, i0, 1) = -zTerm;
241 outputValues(3, i0, 2) = -y*zTerm2;
242 outputValues(3, i0, 3) = 0.0;
243 outputValues(3, i0, 4) = -x*zTerm2;
244 outputValues(3, i0, 5) = -x*y*zTerm3;
246 outputValues(4, i0, 0) = 0.0;
247 outputValues(4, i0, 1) = 0.0;
248 outputValues(4, i0, 2) = 0.0;
249 outputValues(4, i0, 3) = 0.0;
250 outputValues(4, i0, 4) = 0.0;
251 outputValues(4, i0, 5) = 0.0;
266 this -> basisCellTopology_.getDimension() );
267 for(
int dofOrd = 0; dofOrd <
this -> basisCardinality_; dofOrd++) {
268 for (
int i0 = 0; i0 < dim0; i0++) {
269 for(
int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
270 outputValues(dofOrd, i0, dkOrd) = 0.0;
278 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): Invalid operator type");
284 template<
class Scalar,
class ArrayScalar>
286 const ArrayScalar & inputPoints,
287 const ArrayScalar & cellVertices,
289 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
290 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): FEM Basis calling an FVD member function");
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Pyramid cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Basis_HGRAD_PYR_C1_FEM()
Constructor.