Intrepid
Intrepid_HCURL_HEX_I1_FEMDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 namespace Intrepid {
50 
51 template<class Scalar, class ArrayScalar>
53  {
54  this -> basisCardinality_ = 12;
55  this -> basisDegree_ = 1;
56  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
57  this -> basisType_ = BASIS_FEM_DEFAULT;
58  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
59  this -> basisTagsAreSet_ = false;
60  }
61 
62 template<class Scalar, class ArrayScalar>
64 
65  // Basis-dependent intializations
66  int tagSize = 4; // size of DoF tag
67  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
68  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
69  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
70 
71  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
72  int tags[] = {
73  1, 0, 0, 1,
74  1, 1, 0, 1,
75  1, 2, 0, 1,
76  1, 3, 0, 1,
77  1, 4, 0, 1,
78  1, 5, 0, 1,
79  1, 6, 0, 1,
80  1, 7, 0, 1,
81  1, 8, 0, 1,
82  1, 9, 0, 1,
83  1, 10, 0, 1,
84  1, 11, 0, 1 };
85 
86  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
87  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
88  this -> ordinalToTag_,
89  tags,
90  this -> basisCardinality_,
91  tagSize,
92  posScDim,
93  posScOrd,
94  posDfOrd);
95 }
96 
97 
98 
99 template<class Scalar, class ArrayScalar>
101  const ArrayScalar & inputPoints,
102  const EOperator operatorType) const {
103 
104 // Verify arguments
105 #ifdef HAVE_INTREPID_DEBUG
106  Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
107  inputPoints,
108  operatorType,
109  this -> getBaseCellTopology(),
110  this -> getCardinality() );
111 #endif
112 
113  // Number of evaluation points = dim 0 of inputPoints
114  int dim0 = inputPoints.dimension(0);
115 
116  // Temporaries: (x,y,z) coordinates of the evaluation point
117  Scalar x = 0.0;
118  Scalar y = 0.0;
119  Scalar z = 0.0;
120 
121  switch (operatorType) {
122  case OPERATOR_VALUE:
123  for (int i0 = 0; i0 < dim0; i0++) {
124  x = inputPoints(i0, 0);
125  y = inputPoints(i0, 1);
126  z = inputPoints(i0, 2);
127 
128  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
129  outputValues(0, i0, 0) = (1.0 - y)*(1.0 - z)/8.0;
130  outputValues(0, i0, 1) = 0.0;
131  outputValues(0, i0, 2) = 0.0;
132 
133  outputValues(1, i0, 0) = 0.0;
134  outputValues(1, i0, 1) = (1.0 + x)*(1.0 - z)/8.0;
135  outputValues(1, i0, 2) = 0.0;
136 
137  outputValues(2, i0, 0) = -(1.0 + y)*(1.0 - z)/8.0;
138  outputValues(2, i0, 1) = 0.0;
139  outputValues(2, i0, 2) = 0.0;
140 
141  outputValues(3, i0, 0) = 0.0;
142  outputValues(3, i0, 1) = -(1.0 - x)*(1.0 - z)/8.0;
143  outputValues(3, i0, 2) = 0.0;
144 
145  outputValues(4, i0, 0) = (1.0 - y)*(1.0 + z)/8.0;
146  outputValues(4, i0, 1) = 0.0;
147  outputValues(4, i0, 2) = 0.0;
148 
149  outputValues(5, i0, 0) = 0.0;
150  outputValues(5, i0, 1) = (1.0 + x)*(1.0 + z)/8.0;
151  outputValues(5, i0, 2) = 0.0;
152 
153  outputValues(6, i0, 0) = -(1.0 + y)*(1.0 + z)/8.0;
154  outputValues(6, i0, 1) = 0.0;
155  outputValues(6, i0, 2) = 0.0;
156 
157  outputValues(7, i0, 0) = 0.0;
158  outputValues(7, i0, 1) = -(1.0 - x)*(1.0 + z)/8.0;
159  outputValues(7, i0, 2) = 0.0;
160 
161  outputValues(8, i0, 0) = 0.0;
162  outputValues(8, i0, 1) = 0.0;
163  outputValues(8, i0, 2) = (1.0 - x)*(1.0 - y)/8.0;
164 
165  outputValues(9, i0, 0) = 0.0;
166  outputValues(9, i0, 1) = 0.0;
167  outputValues(9, i0, 2) = (1.0 + x)*(1.0 - y)/8.0;
168 
169  outputValues(10, i0, 0) = 0.0;
170  outputValues(10, i0, 1) = 0.0;
171  outputValues(10, i0, 2) = (1.0 + x)*(1.0 + y)/8.0;
172 
173  outputValues(11, i0, 0) = 0.0;
174  outputValues(11, i0, 1) = 0.0;
175  outputValues(11, i0, 2) = (1.0 - x)*(1.0 + y)/8.0;
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  z = inputPoints(i0, 2);
184 
185  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
186  outputValues(0, i0, 0) = 0.0;
187  outputValues(0, i0, 1) = -(1.0 - y)/8.0;
188  outputValues(0, i0, 2) = (1.0 - z)/8.0;
189 
190  outputValues(1, i0, 0) = (1.0 + x)/8.0;
191  outputValues(1, i0, 1) = 0.0;
192  outputValues(1, i0, 2) = (1.0 - z)/8.0;
193 
194  outputValues(2, i0, 0) = 0.0;
195  outputValues(2, i0, 1) = (1.0 + y)/8.0;
196  outputValues(2, i0, 2) = (1.0 - z)/8.0;
197 
198  outputValues(3, i0, 0) = -(1.0 - x)/8.0;
199  outputValues(3, i0, 1) = 0.0;
200  outputValues(3, i0, 2) = (1.0 - z)/8.0;
201 
202  outputValues(4, i0, 0) = 0.0;
203  outputValues(4, i0, 1) = (1.0 - y)/8.0;
204  outputValues(4, i0, 2) = (1.0 + z)/8.0;
205 
206  outputValues(5, i0, 0) = -(1.0 + x)/8.0;
207  outputValues(5, i0, 1) = 0.0;
208  outputValues(5, i0, 2) = (1.0 + z)/8.0;
209 
210  outputValues(6, i0, 0) = 0.0;
211  outputValues(6, i0, 1) = -(1.0 + y)/8.0;
212  outputValues(6, i0, 2) = (1.0 + z)/8.0;
213 
214  outputValues(7, i0, 0) = (1.0 - x)/8.0;
215  outputValues(7, i0, 1) = 0.0;
216  outputValues(7, i0, 2) = (1.0 + z)/8.0;
217 
218  outputValues(8, i0, 0) = -(1.0 - x)/8.0;
219  outputValues(8, i0, 1) = (1.0 - y)/8.0;
220  outputValues(8, i0, 2) = 0.0;
221 
222  outputValues(9, i0, 0) = -(1.0 + x)/8.0;
223  outputValues(9, i0, 1) = -(1.0 - y)/8.0;
224  outputValues(9, i0, 2) = 0.0;
225 
226  outputValues(10, i0, 0) = (1.0 + x)/8.0;
227  outputValues(10, i0, 1) = -(1.0 + y)/8.0;
228  outputValues(10, i0, 2) = 0.0;
229 
230  outputValues(11, i0, 0) = (1.0 - x)/8.0;
231  outputValues(11, i0, 1) = (1.0 + y)/8.0;
232  outputValues(11, i0, 2) = 0.0;
233  }
234  break;
235 
236  case OPERATOR_DIV:
237  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
238  ">>> ERROR (Basis_HCURL_HEX_I1_FEM): DIV is invalid operator for HCURL Basis Functions");
239  break;
240 
241  case OPERATOR_GRAD:
242  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
243  ">>> ERROR (Basis_HCURL_HEX_I1_FEM): GRAD is invalid operator for HCURL Basis Functions");
244  break;
245 
246  case OPERATOR_D1:
247  case OPERATOR_D2:
248  case OPERATOR_D3:
249  case OPERATOR_D4:
250  case OPERATOR_D5:
251  case OPERATOR_D6:
252  case OPERATOR_D7:
253  case OPERATOR_D8:
254  case OPERATOR_D9:
255  case OPERATOR_D10:
256  TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
257  (operatorType == OPERATOR_D2) ||
258  (operatorType == OPERATOR_D3) ||
259  (operatorType == OPERATOR_D4) ||
260  (operatorType == OPERATOR_D5) ||
261  (operatorType == OPERATOR_D6) ||
262  (operatorType == OPERATOR_D7) ||
263  (operatorType == OPERATOR_D8) ||
264  (operatorType == OPERATOR_D9) ||
265  (operatorType == OPERATOR_D10) ),
266  std::invalid_argument,
267  ">>> ERROR (Basis_HCURL_HEX_I1_FEM): Invalid operator type");
268  break;
269 
270  default:
271  TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
272  (operatorType != OPERATOR_GRAD) &&
273  (operatorType != OPERATOR_CURL) &&
274  (operatorType != OPERATOR_DIV) &&
275  (operatorType != OPERATOR_D1) &&
276  (operatorType != OPERATOR_D2) &&
277  (operatorType != OPERATOR_D3) &&
278  (operatorType != OPERATOR_D4) &&
279  (operatorType != OPERATOR_D5) &&
280  (operatorType != OPERATOR_D6) &&
281  (operatorType != OPERATOR_D7) &&
282  (operatorType != OPERATOR_D8) &&
283  (operatorType != OPERATOR_D9) &&
284  (operatorType != OPERATOR_D10) ),
285  std::invalid_argument,
286  ">>> ERROR (Basis_HCURL_HEX_I1_FEM): Invalid operator type");
287  }
288 }
289 
290 
291 
292 template<class Scalar, class ArrayScalar>
294  const ArrayScalar & inputPoints,
295  const ArrayScalar & cellVertices,
296  const EOperator operatorType) const {
297  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
298  ">>> ERROR (Basis_HCURL_HEX_I1_FEM): FEM Basis calling an FVD member function");
299  }
300 
301 template<class Scalar, class ArrayScalar>
303 #ifdef HAVE_INTREPID_DEBUG
304  // Verify rank of output array.
305  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
306  ">>> ERROR: (Intrepid::Basis_HCURL_QUAD_I1_FEM::getDofCoords) rank = 2 required for DofCoords array");
307  // Verify 0th dimension of output array.
308  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
309  ">>> ERROR: (Intrepid::Basis_HCURL_QUAD_I1_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
310  // Verify 1st dimension of output array.
311  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
312  ">>> ERROR: (Intrepid::Basis_HCURL_QUAD_I1_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
313 #endif
314 
315  DofCoords(0,0) = 0.0; DofCoords(0,1) = -1.0; DofCoords(0,2) = -1.0;
316  DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0; DofCoords(1,2) = -1.0;
317  DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = -1.0;
318  DofCoords(3,0) = -1.0; DofCoords(3,1) = 0.0; DofCoords(3,2) = -1.0;
319  DofCoords(4,0) = 0.0; DofCoords(4,1) = -1.0; DofCoords(4,2) = 1.0;
320  DofCoords(5,0) = 1.0; DofCoords(5,1) = 0.0; DofCoords(5,2) = 1.0;
321  DofCoords(6,0) = 0.0; DofCoords(6,1) = 1.0; DofCoords(6,2) = 1.0;
322  DofCoords(7,0) = -1.0; DofCoords(7,1) = 0.0; DofCoords(7,2) = 1.0;
323  DofCoords(8,0) = -1.0; DofCoords(8,1) = -1.0; DofCoords(8,2) = 0.0;
324  DofCoords(9,0) = 1.0; DofCoords(9,1) = -1.0; DofCoords(9,2) = 0.0;
325  DofCoords(10,0) = 1.0; DofCoords(10,1) = 1.0; DofCoords(10,2) = 0.0;
326  DofCoords(11,0) = -1.0; DofCoords(11,1) = 1.0; DofCoords(11,2) = 0.0;
327 }
328 
329 }// namespace Intrepid
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.