1 #ifndef INTREPID_HGRAD_QUAD_C2_FEMDEF_HPP 2 #define INTREPID_HGRAD_QUAD_C2_FEMDEF_HPP 53 template<
class Scalar,
class ArrayScalar>
56 this -> basisCardinality_ = 9;
57 this -> basisDegree_ = 2;
58 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
59 this -> basisType_ = BASIS_FEM_DEFAULT;
60 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
61 this -> basisTagsAreSet_ =
false;
65 template<
class Scalar,
class ArrayScalar>
75 int tags[] = { 0, 0, 0, 1,
89 this -> ordinalToTag_,
91 this -> basisCardinality_,
100 template<
class Scalar,
class ArrayScalar>
102 const ArrayScalar & inputPoints,
106 #ifdef HAVE_INTREPID_DEBUG 107 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
110 this -> getBaseCellTopology(),
111 this -> getCardinality() );
115 int dim0 = inputPoints.dimension(0);
121 switch (operatorType) {
124 for (
int i0 = 0; i0 < dim0; i0++) {
125 x = inputPoints(i0, 0);
126 y = inputPoints(i0, 1);
129 outputValues(0, i0) = x*(x - 1.0)*y*(y - 1.0)/4.0;
130 outputValues(1, i0) = x*(x + 1.0)*y*(y - 1.0)/4.0;
131 outputValues(2, i0) = x*(x + 1.0)*y*(y + 1.0)/4.0;
132 outputValues(3, i0) = x*(x - 1.0)*y*(y + 1.0)/4.0;
134 outputValues(4, i0) = (1.0 - x)*(1.0 + x)*y*(y - 1.0)/2.0;
135 outputValues(5, i0) = x*(x + 1.0)*(1.0 - y)*(1.0 + y)/2.0;
136 outputValues(6, i0) = (1.0 - x)*(1.0 + x)*y*(y + 1.0)/2.0;
137 outputValues(7, i0) = x*(x - 1.0)*(1.0 - y)*(1.0 + y)/2.0;
139 outputValues(8, i0) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y);
145 for (
int i0 = 0; i0 < dim0; i0++) {
146 x = inputPoints(i0,0);
147 y = inputPoints(i0,1);
150 outputValues(0, i0, 0) = (-0.25 + 0.5*x)*(-1. + y)*y;
151 outputValues(0, i0, 1) = (-1.0 + x)*x*(-0.25 + 0.5*y);
153 outputValues(1, i0, 0) = (0.25 + 0.5*x)*(-1. + y)*y;
154 outputValues(1, i0, 1) = x*(1. + x)*(-0.25 + 0.5*y);
156 outputValues(2, i0, 0) = (0.25 + 0.5*x)*y*(1. + y);
157 outputValues(2, i0, 1) = x*(1. + x)*(0.25 + 0.5*y);
159 outputValues(3, i0, 0) = (-0.25 + 0.5*x)*y*(1. + y);
160 outputValues(3, i0, 1) = (-1. + x)*x*(0.25 + 0.5*y);
162 outputValues(4, i0, 0) = x*(1.0 - y)*y;
163 outputValues(4, i0, 1) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
165 outputValues(5, i0, 0) = 0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
166 outputValues(5, i0, 1) =-x*(1.0 + x)*y;
168 outputValues(6, i0, 0) =-y*(1.0 + y)*x;
169 outputValues(6, i0, 1) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
171 outputValues(7, i0, 0) = 0.5*(1.0 - y)*(1.0+ y)*(-1.0 + 2.0*x);
172 outputValues(7, i0, 1) = (1.0 - x)*x*y;
174 outputValues(8, i0, 0) =-2.0*(1.0 - y)*(1.0 + y)*x;
175 outputValues(8, i0, 1) =-2.0*(1.0 - x)*(1.0 + x)*y;
180 for (
int i0 = 0; i0 < dim0; i0++) {
181 x = inputPoints(i0,0);
182 y = inputPoints(i0,1);
186 outputValues(0, i0, 1) =-(-0.25 + 0.5*x)*(-1. + y)*y;
187 outputValues(0, i0, 0) = (-1.0 + x)*x*(-0.25 + 0.5*y);
189 outputValues(1, i0, 1) =-(0.25 + 0.5*x)*(-1. + y)*y;
190 outputValues(1, i0, 0) = x*(1. + x)*(-0.25 + 0.5*y);
192 outputValues(2, i0, 1) =-(0.25 + 0.5*x)*y*(1. + y);
193 outputValues(2, i0, 0) = x*(1. + x)*(0.25 + 0.5*y);
195 outputValues(3, i0, 1) =-(-0.25 + 0.5*x)*y*(1. + y);
196 outputValues(3, i0, 0) = (-1. + x)*x*(0.25 + 0.5*y);
198 outputValues(4, i0, 1) =-x*(1.0 - y)*y;
199 outputValues(4, i0, 0) = 0.5*(1.0 - x)*(1.0 + x)*(-1.0 + 2.0*y);
201 outputValues(5, i0, 1) =-0.5*(1.0 - y)*(1.0 + y)*(1.0 + 2.0*x);
202 outputValues(5, i0, 0) =-x*(1.0 + x)*y;
204 outputValues(6, i0, 1) = y*(1.0 + y)*x;
205 outputValues(6, i0, 0) = 0.5*(1.0 - x)*(1.0 + x)*(1.0 + 2.0*y);
207 outputValues(7, i0, 1) =-0.5*(1.0 - y)*(1.0 + y)*(-1.0 + 2.0*x);
208 outputValues(7, i0, 0) = (1.0 - x)*x*y;
210 outputValues(8, i0, 1) = 2.0*(1.0 - y)*(1.0 + y)*x;
211 outputValues(8, i0, 0) =-2.0*(1.0 - x)*(1.0 + x)*y;
216 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
217 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
221 for (
int i0 = 0; i0 < dim0; i0++) {
222 x = inputPoints(i0,0);
223 y = inputPoints(i0,1);
226 outputValues(0, i0, 0) = 0.5*(-1.0 + y)*y;
227 outputValues(0, i0, 1) = 0.25 - 0.5*y + x*(-0.5 + 1.*y);
228 outputValues(0, i0, 2) = 0.5*(-1.0 + x)*x;
230 outputValues(1, i0, 0) = 0.5*(-1.0 + y)*y;
231 outputValues(1, i0, 1) =-0.25 + 0.5*y + x*(-0.5 + 1.*y);
232 outputValues(1, i0, 2) = 0.5*x*(1.0 + x);
234 outputValues(2, i0, 0) = 0.5*y*(1.0 + y);
235 outputValues(2, i0, 1) = 0.25 + 0.5*y + x*(0.5 + 1.*y);
236 outputValues(2, i0, 2) = 0.5*x*(1.0 + x);
238 outputValues(3, i0, 0) = 0.5*y*(1.0 + y);
239 outputValues(3, i0, 1) =-0.25 - 0.5*y + x*(0.5 + 1.*y);
240 outputValues(3, i0, 2) = 0.5*(-1.0 + x)*x;
242 outputValues(4, i0, 0) = (1.0 - y)*y;
243 outputValues(4, i0, 1) = x*(1. - 2.*y);
244 outputValues(4, i0, 2) = (1.0 - x)*(1.0 + x);
246 outputValues(5, i0, 0) = (1.0 - y)*(1.0 + y);
247 outputValues(5, i0, 1) = x*(0. - 2.*y) - 1.*y;
248 outputValues(5, i0, 2) =-x*(1.0 + x);
250 outputValues(6, i0, 0) =-y*(1.0 + y);
251 outputValues(6, i0, 1) = x*(-1. - 2.*y);
252 outputValues(6, i0, 2) = (1.0 - x)*(1.0 + x);
254 outputValues(7, i0, 0) = (1.0 - y)*(1.0 + y);
255 outputValues(7, i0, 1) = x*(0. - 2.*y) + 1.*y;
256 outputValues(7, i0, 2) = (1.0 - x)*x;
258 outputValues(8, i0, 0) =-2.0 + 2.0*y*y;
259 outputValues(8, i0, 1) = 4*x*y;
260 outputValues(8, i0, 2) =-2.0 + 2.0*x*x;
267 for (
int i0 = 0; i0 < dim0; i0++) {
268 x = inputPoints(i0,0);
269 y = inputPoints(i0,1);
271 outputValues(0, i0, 0) = 0.0;
272 outputValues(0, i0, 1) =-0.5 + y;
273 outputValues(0, i0, 2) =-0.5 + x;
274 outputValues(0, i0, 3) = 0.0;
276 outputValues(1, i0, 0) = 0.0;
277 outputValues(1, i0, 1) =-0.5 + y;
278 outputValues(1, i0, 2) = 0.5 + x;
279 outputValues(1, i0, 3) = 0.0;
281 outputValues(2, i0, 0) = 0.0;
282 outputValues(2, i0, 1) = 0.5 + y;
283 outputValues(2, i0, 2) = 0.5 + x;
284 outputValues(2, i0, 3) = 0.0;
286 outputValues(3, i0, 0) = 0.0;
287 outputValues(3, i0, 1) = 0.5 + y;
288 outputValues(3, i0, 2) =-0.5 + x;
289 outputValues(3, i0, 3) = 0.0;
291 outputValues(4, i0, 0) = 0.0;
292 outputValues(4, i0, 1) = 1.0 - 2.0*y;
293 outputValues(4, i0, 2) =-2.0*x;
294 outputValues(4, i0, 3) = 0.0;
296 outputValues(5, i0, 0) = 0.0;
297 outputValues(5, i0, 1) =-2.0*y;
298 outputValues(5, i0, 2) =-1.0 - 2.0*x;
299 outputValues(5, i0, 3) = 0.0;
301 outputValues(6, i0, 0) = 0.0;
302 outputValues(6, i0, 1) =-1.0 - 2.0*y;
303 outputValues(6, i0, 2) =-2.0*x;
304 outputValues(6, i0, 3) = 0.0;
306 outputValues(7, i0, 0) = 0.0;
307 outputValues(7, i0, 1) =-2.0*y;
308 outputValues(7, i0, 2) = 1.0 - 2.0*x;
309 outputValues(7, i0, 3) = 0.0;
311 outputValues(8, i0, 0) = 0.0;
312 outputValues(8, i0, 1) = 4.0*y;
313 outputValues(8, i0, 2) = 4.0*x;
314 outputValues(8, i0, 3) = 0.0;
320 for (
int i0 = 0; i0 < dim0; i0++) {
322 outputValues(0, i0, 0) = 0.0;
323 outputValues(0, i0, 1) = 0.0;
324 outputValues(0, i0, 2) = 1.0;
325 outputValues(0, i0, 3) = 0.0;
326 outputValues(0, i0, 4) = 0.0;
328 outputValues(1, i0, 0) = 0.0;
329 outputValues(1, i0, 1) = 0.0;
330 outputValues(1, i0, 2) = 1.0;
331 outputValues(1, i0, 3) = 0.0;
332 outputValues(1, i0, 4) = 0.0;
334 outputValues(2, i0, 0) = 0.0;
335 outputValues(2, i0, 1) = 0.0;
336 outputValues(2, i0, 2) = 1.0;
337 outputValues(2, i0, 3) = 0.0;
338 outputValues(2, i0, 4) = 0.0;
340 outputValues(3, i0, 0) = 0.0;
341 outputValues(3, i0, 1) = 0.0;
342 outputValues(3, i0, 2) = 1.0;
343 outputValues(3, i0, 3) = 0.0;
344 outputValues(3, i0, 4) = 0.0;
346 outputValues(4, i0, 0) = 0.0;
347 outputValues(4, i0, 1) = 0.0;
348 outputValues(4, i0, 2) =-2.0;
349 outputValues(4, i0, 3) = 0.0;
350 outputValues(4, i0, 4) = 0.0;
352 outputValues(5, i0, 0) = 0.0;
353 outputValues(5, i0, 1) = 0.0;
354 outputValues(5, i0, 2) =-2.0;
355 outputValues(5, i0, 3) = 0.0;
356 outputValues(5, i0, 4) = 0.0;
358 outputValues(6, i0, 0) = 0.0;
359 outputValues(6, i0, 1) = 0.0;
360 outputValues(6, i0, 2) =-2.0;
361 outputValues(6, i0, 3) = 0.0;
362 outputValues(6, i0, 4) = 0.0;
364 outputValues(7, i0, 0) = 0.0;
365 outputValues(7, i0, 1) = 0.0;
366 outputValues(7, i0, 2) =-2.0;
367 outputValues(7, i0, 3) = 0.0;
368 outputValues(7, i0, 4) = 0.0;
370 outputValues(8, i0, 0) = 0.0;
371 outputValues(8, i0, 1) = 0.0;
372 outputValues(8, i0, 2) = 4.0;
373 outputValues(8, i0, 3) = 0.0;
374 outputValues(8, i0, 4) = 0.0;
387 this -> basisCellTopology_.getDimension() );
388 for(
int dofOrd = 0; dofOrd <
this -> basisCardinality_; dofOrd++) {
389 for (
int i0 = 0; i0 < dim0; i0++) {
390 for(
int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
391 outputValues(dofOrd, i0, dkOrd) = 0.0;
400 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type");
406 template<
class Scalar,
class ArrayScalar>
408 const ArrayScalar & inputPoints,
409 const ArrayScalar & cellVertices,
411 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
412 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): FEM Basis calling an FVD member function");
417 template<
class Scalar,
class ArrayScalar>
419 #ifdef HAVE_INTREPID_DEBUG 421 TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
422 ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) rank = 2 required for DofCoords array");
424 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) ==
this -> basisCardinality_ ), std::invalid_argument,
425 ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
427 TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(
this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
428 ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
431 DofCoords(0,0) = -1.0; DofCoords(0,1) = -1.0;
432 DofCoords(1,0) = 1.0; DofCoords(1,1) = -1.0;
433 DofCoords(2,0) = 1.0; DofCoords(2,1) = 1.0;
434 DofCoords(3,0) = -1.0; DofCoords(3,1) = 1.0;
436 DofCoords(4,0) = 0.0; DofCoords(4,1) = -1.0;
437 DofCoords(5,0) = 1.0; DofCoords(5,1) = 0.0;
438 DofCoords(6,0) = 0.0; DofCoords(6,1) = 1.0;
439 DofCoords(7,0) = -1.0; DofCoords(7,1) = 0.0;
441 DofCoords(8,0) = 0.0; DofCoords(8,1) = 0.0;
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...
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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.
Basis_HGRAD_QUAD_C2_FEM()
Constructor.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.