Intrepid
Intrepid_HGRAD_QUAD_C1_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_QUAD_C1_FEMDEF_HPP
2 #define INTREPID_HGRAD_QUAD_C1_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_ = 4;
57  this -> basisDegree_ = 1;
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 
80  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
81  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
82  this -> ordinalToTag_,
83  tags,
84  this -> basisCardinality_,
85  tagSize,
86  posScDim,
87  posScOrd,
88  posDfOrd);
89 }
90 
91 
92 
93 template<class Scalar, class ArrayScalar>
95  const ArrayScalar & inputPoints,
96  const EOperator operatorType) const {
97 
98  // Verify arguments
99 #ifdef HAVE_INTREPID_DEBUG
100  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
101  inputPoints,
102  operatorType,
103  this -> getBaseCellTopology(),
104  this -> getCardinality() );
105 #endif
106 
107  // Number of evaluation points = dim 0 of inputPoints
108  int dim0 = inputPoints.dimension(0);
109 
110  // Temporaries: (x,y) coordinates of the evaluation point
111  Scalar x = 0.0;
112  Scalar y = 0.0;
113 
114  switch (operatorType) {
115 
116  case OPERATOR_VALUE:
117  for (int i0 = 0; i0 < dim0; i0++) {
118  x = inputPoints(i0, 0);
119  y = inputPoints(i0, 1);
120 
121  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
122  outputValues(0, i0) = (1.0 - x)*(1.0 - y)/4.0;
123  outputValues(1, i0) = (1.0 + x)*(1.0 - y)/4.0;
124  outputValues(2, i0) = (1.0 + x)*(1.0 + y)/4.0;
125  outputValues(3, i0) = (1.0 - x)*(1.0 + y)/4.0;
126  }
127  break;
128 
129  case OPERATOR_GRAD:
130  case OPERATOR_D1:
131  for (int i0 = 0; i0 < dim0; i0++) {
132  x = inputPoints(i0,0);
133  y = inputPoints(i0,1);
134 
135  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
136  outputValues(0, i0, 0) = -(1.0 - y)/4.0;
137  outputValues(0, i0, 1) = -(1.0 - x)/4.0;
138 
139  outputValues(1, i0, 0) = (1.0 - y)/4.0;
140  outputValues(1, i0, 1) = -(1.0 + x)/4.0;
141 
142  outputValues(2, i0, 0) = (1.0 + y)/4.0;
143  outputValues(2, i0, 1) = (1.0 + x)/4.0;
144 
145  outputValues(3, i0, 0) = -(1.0 + y)/4.0;
146  outputValues(3, i0, 1) = (1.0 - x)/4.0;
147  }
148  break;
149 
150  case OPERATOR_CURL:
151  for (int i0 = 0; i0 < dim0; i0++) {
152  x = inputPoints(i0,0);
153  y = inputPoints(i0,1);
154 
155  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
156  outputValues(0, i0, 0) = -(1.0 - x)/4.0;
157  outputValues(0, i0, 1) = (1.0 - y)/4.0;
158 
159  outputValues(1, i0, 0) = -(1.0 + x)/4.0;
160  outputValues(1, i0, 1) = -(1.0 - y)/4.0;
161 
162  outputValues(2, i0, 0) = (1.0 + x)/4.0;
163  outputValues(2, i0, 1) = -(1.0 + y)/4.0;
164 
165  outputValues(3, i0, 0) = (1.0 - x)/4.0;
166  outputValues(3, i0, 1) = (1.0 + y)/4.0;
167  }
168  break;
169 
170  case OPERATOR_DIV:
171  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
172  ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 2D");
173  break;
174 
175  case OPERATOR_D2:
176  for (int i0 = 0; i0 < dim0; i0++) {
177 
178  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality=3)
179  outputValues(0, i0, 0) = 0.0;
180  outputValues(0, i0, 1) = 0.25;
181  outputValues(0, i0, 2) = 0.0;
182 
183  outputValues(1, i0, 0) = 0.0;
184  outputValues(1, i0, 1) = -0.25;
185  outputValues(1, i0, 2) = 0.0;
186 
187  outputValues(2, i0, 0) = 0.0;
188  outputValues(2, i0, 1) = 0.25;
189  outputValues(2, i0, 2) = 0.0;
190 
191  outputValues(3, i0, 0) = 0.0;
192  outputValues(3, i0, 1) = -0.25;
193  outputValues(3, i0, 2) = 0.0;
194  }
195  break;
196 
197  case OPERATOR_D3:
198  case OPERATOR_D4:
199  case OPERATOR_D5:
200  case OPERATOR_D6:
201  case OPERATOR_D7:
202  case OPERATOR_D8:
203  case OPERATOR_D9:
204  case OPERATOR_D10:
205  {
206  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
207  int DkCardinality = Intrepid::getDkCardinality(operatorType,
208  this -> basisCellTopology_.getDimension() );
209  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
210  for (int i0 = 0; i0 < dim0; i0++) {
211  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
212  outputValues(dofOrd, i0, dkOrd) = 0.0;
213  }
214  }
215  }
216  }
217  break;
218 
219  default:
220  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
221  ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): Invalid operator type");
222  }
223 }
224 
225 
226 
227 template<class Scalar, class ArrayScalar>
229  const ArrayScalar & inputPoints,
230  const ArrayScalar & cellVertices,
231  const EOperator operatorType) const {
232  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
233  ">>> ERROR (Basis_HGRAD_QUAD_C1_FEM): FEM Basis calling an FVD member function");
234 }
235 
236 
237 
238 template<class Scalar, class ArrayScalar>
240 #ifdef HAVE_INTREPID_DEBUG
241  // Verify rank of output array.
242  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
243  ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C1_FEM::getDofCoords) rank = 2 required for DofCoords array");
244  // Verify 0th dimension of output array.
245  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
246  ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C1_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
247  // Verify 1st dimension of output array.
248  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
249  ">>> ERROR: (Intrepid::Basis_HGRAD_QUAD_C1_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
250 #endif
251 
252  DofCoords(0,0) = -1.0; DofCoords(0,1) = -1.0;
253  DofCoords(1,0) = 1.0; DofCoords(1,1) = -1.0;
254  DofCoords(2,0) = 1.0; DofCoords(2,1) = 1.0;
255  DofCoords(3,0) = -1.0; DofCoords(3,1) = 1.0;
256 }
257 
258 }// namespace Intrepid
259 #endif
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.