Intrepid
Intrepid_HGRAD_QUAD_C2_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_QUAD_C2_FEMDEF_HPP
2 #define INTREPID_HGRAD_QUAD_C2_FEMDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53  template<class Scalar, class ArrayScalar>
55  {
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;
62  }
63 
64 
65 template<class Scalar, class ArrayScalar>
67 
68  // Basis-dependent intializations
69  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
70  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
71  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
72  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
73 
74  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
75  int tags[] = { 0, 0, 0, 1,
76  0, 1, 0, 1,
77  0, 2, 0, 1,
78  0, 3, 0, 1,
79  // edge midpoints
80  1, 0, 0, 1,
81  1, 1, 0, 1,
82  1, 2, 0, 1,
83  1, 3, 0, 1,
84  // quad center
85  2, 0, 0, 1};
86 
87  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
88  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
89  this -> ordinalToTag_,
90  tags,
91  this -> basisCardinality_,
92  tagSize,
93  posScDim,
94  posScOrd,
95  posDfOrd);
96 }
97 
98 
99 
100 template<class Scalar, class ArrayScalar>
102  const ArrayScalar & inputPoints,
103  const EOperator operatorType) const {
104 
105  // Verify arguments
106 #ifdef HAVE_INTREPID_DEBUG
107  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
108  inputPoints,
109  operatorType,
110  this -> getBaseCellTopology(),
111  this -> getCardinality() );
112 #endif
113 
114  // Number of evaluation points = dim 0 of inputPoints
115  int dim0 = inputPoints.dimension(0);
116 
117  // Temporaries: (x,y) coordinates of the evaluation point
118  Scalar x = 0.0;
119  Scalar y = 0.0;
120 
121  switch (operatorType) {
122 
123  case OPERATOR_VALUE:
124  for (int i0 = 0; i0 < dim0; i0++) {
125  x = inputPoints(i0, 0);
126  y = inputPoints(i0, 1);
127 
128  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
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;
133  // edge midpoints basis functions
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;
138  // quad bubble basis function
139  outputValues(8, i0) = (1.0 - x)*(1.0 + x)*(1.0 - y)*(1.0 + y);
140  }
141  break;
142 
143  case OPERATOR_GRAD:
144  case OPERATOR_D1:
145  for (int i0 = 0; i0 < dim0; i0++) {
146  x = inputPoints(i0,0);
147  y = inputPoints(i0,1);
148 
149  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
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);
152 
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);
155 
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);
158 
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);
161 
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);
164 
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;
167 
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);
170 
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;
173 
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;
176  }
177  break;
178 
179  case OPERATOR_CURL:
180  for (int i0 = 0; i0 < dim0; i0++) {
181  x = inputPoints(i0,0);
182  y = inputPoints(i0,1);
183 
184  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
185  // CURL(u) = (u_y, -u_x), is rotated GRAD
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);
188 
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);
191 
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);
194 
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);
197 
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);
200 
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;
203 
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);
206 
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;
209 
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;
212  }
213  break;
214 
215  case OPERATOR_DIV:
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");
218  break;
219 
220  case OPERATOR_D2:
221  for (int i0 = 0; i0 < dim0; i0++) {
222  x = inputPoints(i0,0);
223  y = inputPoints(i0,1);
224 
225  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality=3)
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;
229 
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);
233 
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);
237 
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;
241 
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);
245 
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);
249 
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);
253 
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;
257 
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;
261 
262  }
263  break;
264 
265  case OPERATOR_D3:
266  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D3Cardinality=4)
267  for (int i0 = 0; i0 < dim0; i0++) {
268  x = inputPoints(i0,0);
269  y = inputPoints(i0,1);
270 
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;
275 
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;
280 
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;
285 
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;
290 
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;
295 
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;
300 
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;
305 
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;
310 
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;
315  }
316  break;
317 
318  case OPERATOR_D4:
319  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D4Cardinality=5)
320  for (int i0 = 0; i0 < dim0; i0++) {
321 
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;
327 
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;
333 
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;
339 
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;
345 
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;
351 
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;
357 
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;
363 
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;
369 
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;
375  }
376  break;
377 
378  case OPERATOR_D5:
379  case OPERATOR_D6:
380  case OPERATOR_D7:
381  case OPERATOR_D8:
382  case OPERATOR_D9:
383  case OPERATOR_D10:
384  {
385  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
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;
392  }
393  }
394  }
395  }
396  break;
397 
398  default:
399  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
400  ">>> ERROR (Basis_HGRAD_QUAD_C2_FEM): Invalid operator type");
401  }
402 }
403 
404 
405 
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");
413 }
414 
415 
416 
417 template<class Scalar, class ArrayScalar>
419 #ifdef HAVE_INTREPID_DEBUG
420  // Verify rank of output array.
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");
423  // Verify 0th dimension of output 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");
426  // Verify 1st dimension of output 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");
429 #endif
430 
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;
435 
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;
440 
441  DofCoords(8,0) = 0.0; DofCoords(8,1) = 0.0;
442 
443 }
444 
445 }// namespace Intrepid
446 #endif
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.