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,
88 Intrepid::setOrdinalTagData(
this -> tagToOrdinal_,
89 this -> ordinalToTag_,
91 this -> basisCardinality_,
100 template<
class Scalar,
class ArrayScalar>
102 const ArrayScalar & inputPoints,
103 const EOperator operatorType)
const {
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;
386 int DkCardinality = Intrepid::getDkCardinality(operatorType,
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;
399 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
400 ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type");
406 template<
class Scalar,
class ArrayScalar>
408 const ArrayScalar & inputPoints,
409 const ArrayScalar & cellVertices,
410 const EOperator operatorType)
const {
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;
448 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
450 #warning "The Intrepid package is deprecated"
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.
Basis_HGRAD_QUAD_C2_FEM()
Constructor.
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.