Intrepid
Intrepid_HGRAD_HEX_C1_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_HEX_C1_FEMDEF_HPP
2 #define INTREPID_HGRAD_HEX_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 
54 template<class Scalar, class ArrayScalar>
55 Basis_HGRAD_HEX_C1_FEM<Scalar,ArrayScalar>::Basis_HGRAD_HEX_C1_FEM()
56  {
57  this -> basisCardinality_ = 8;
58  this -> basisDegree_ = 1;
59  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
60  this -> basisType_ = BASIS_FEM_DEFAULT;
61  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62  this -> basisTagsAreSet_ = false;
63  }
64 
65 
66 
67 template<class Scalar, class ArrayScalar>
68 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::initializeTags() {
69 
70  // Basis-dependent intializations
71  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
72  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
73  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
74  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
75 
76  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
77  int tags[] = { 0, 0, 0, 1,
78  0, 1, 0, 1,
79  0, 2, 0, 1,
80  0, 3, 0, 1,
81  0, 4, 0, 1,
82  0, 5, 0, 1,
83  0, 6, 0, 1,
84  0, 7, 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>
100 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues,
101  const ArrayScalar & inputPoints,
102  const EOperator operatorType) const {
103 
104  // Verify arguments
105 #ifdef HAVE_INTREPID_DEBUG
106  Intrepid::getValues_HGRAD_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 
123  case OPERATOR_VALUE:
124  for (int i0 = 0; i0 < dim0; i0++) {
125  x = inputPoints(i0, 0);
126  y = inputPoints(i0, 1);
127  z = inputPoints(i0, 2);
128 
129  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
130  outputValues(0, i0) = (1.0 - x)*(1.0 - y)*(1.0 - z)/8.0;
131  outputValues(1, i0) = (1.0 + x)*(1.0 - y)*(1.0 - z)/8.0;
132  outputValues(2, i0) = (1.0 + x)*(1.0 + y)*(1.0 - z)/8.0;
133  outputValues(3, i0) = (1.0 - x)*(1.0 + y)*(1.0 - z)/8.0;
134 
135  outputValues(4, i0) = (1.0 - x)*(1.0 - y)*(1.0 + z)/8.0;
136  outputValues(5, i0) = (1.0 + x)*(1.0 - y)*(1.0 + z)/8.0;
137  outputValues(6, i0) = (1.0 + x)*(1.0 + y)*(1.0 + z)/8.0;
138  outputValues(7, i0) = (1.0 - x)*(1.0 + y)*(1.0 + z)/8.0;
139  }
140  break;
141 
142  case OPERATOR_GRAD:
143  case OPERATOR_D1:
144  for (int i0 = 0; i0 < dim0; i0++) {
145  x = inputPoints(i0,0);
146  y = inputPoints(i0,1);
147  z = inputPoints(i0,2);
148 
149  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
150  outputValues(0, i0, 0) = -(1.0 - y)*(1.0 - z)/8.0;
151  outputValues(0, i0, 1) = -(1.0 - x)*(1.0 - z)/8.0;
152  outputValues(0, i0, 2) = -(1.0 - x)*(1.0 - y)/8.0;
153 
154  outputValues(1, i0, 0) = (1.0 - y)*(1.0 - z)/8.0;
155  outputValues(1, i0, 1) = -(1.0 + x)*(1.0 - z)/8.0;
156  outputValues(1, i0, 2) = -(1.0 + x)*(1.0 - y)/8.0;
157 
158  outputValues(2, i0, 0) = (1.0 + y)*(1.0 - z)/8.0;
159  outputValues(2, i0, 1) = (1.0 + x)*(1.0 - z)/8.0;
160  outputValues(2, i0, 2) = -(1.0 + x)*(1.0 + y)/8.0;
161 
162  outputValues(3, i0, 0) = -(1.0 + y)*(1.0 - z)/8.0;
163  outputValues(3, i0, 1) = (1.0 - x)*(1.0 - z)/8.0;
164  outputValues(3, i0, 2) = -(1.0 - x)*(1.0 + y)/8.0;
165 
166  outputValues(4, i0, 0) = -(1.0 - y)*(1.0 + z)/8.0;
167  outputValues(4, i0, 1) = -(1.0 - x)*(1.0 + z)/8.0;
168  outputValues(4, i0, 2) = (1.0 - x)*(1.0 - y)/8.0;
169 
170  outputValues(5, i0, 0) = (1.0 - y)*(1.0 + z)/8.0;
171  outputValues(5, i0, 1) = -(1.0 + x)*(1.0 + z)/8.0;
172  outputValues(5, i0, 2) = (1.0 + x)*(1.0 - y)/8.0;
173 
174  outputValues(6, i0, 0) = (1.0 + y)*(1.0 + z)/8.0;
175  outputValues(6, i0, 1) = (1.0 + x)*(1.0 + z)/8.0;
176  outputValues(6, i0, 2) = (1.0 + x)*(1.0 + y)/8.0;
177 
178  outputValues(7, i0, 0) = -(1.0 + y)*(1.0 + z)/8.0;
179  outputValues(7, i0, 1) = (1.0 - x)*(1.0 + z)/8.0;
180  outputValues(7, i0, 2) = (1.0 - x)*(1.0 + y)/8.0;
181  }
182  break;
183 
184  case OPERATOR_CURL:
185  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
186  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
187  break;
188 
189  case OPERATOR_DIV:
190  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
191  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
192  break;
193 
194  case OPERATOR_D2:
195  for (int i0 = 0; i0 < dim0; i0++) {
196  x = inputPoints(i0,0);
197  y = inputPoints(i0,1);
198  z = inputPoints(i0,2);
199 
200  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6)
201  outputValues(0, i0, 0) = 0.0; // {2, 0, 0}
202  outputValues(0, i0, 1) = (1.0 - z)/8.0; // {1, 1, 0}
203  outputValues(0, i0, 2) = (1.0 - y)/8.0; // {1, 0, 1}
204  outputValues(0, i0, 3) = 0.0; // {0, 2, 0}
205  outputValues(0, i0, 4) = (1.0 - x)/8.0; // {0, 1, 1}
206  outputValues(0, i0, 5) = 0.0; // {0, 0, 2}
207 
208  outputValues(1, i0, 0) = 0.0; // {2, 0, 0}
209  outputValues(1, i0, 1) = -(1.0 - z)/8.0; // {1, 1, 0}
210  outputValues(1, i0, 2) = -(1.0 - y)/8.0; // {1, 0, 1}
211  outputValues(1, i0, 3) = 0.0; // {0, 2, 0}
212  outputValues(1, i0, 4) = (1.0 + x)/8.0; // {0, 1, 1}
213  outputValues(1, i0, 5) = 0.0; // {0, 0, 2}
214 
215  outputValues(2, i0, 0) = 0.0; // {2, 0, 0}
216  outputValues(2, i0, 1) = (1.0 - z)/8.0; // {1, 1, 0}
217  outputValues(2, i0, 2) = -(1.0 + y)/8.0; // {1, 0, 1}
218  outputValues(2, i0, 3) = 0.0; // {0, 2, 0}
219  outputValues(2, i0, 4) = -(1.0 + x)/8.0; // {0, 1, 1}
220  outputValues(2, i0, 5) = 0.0; // {0, 0, 2}
221 
222  outputValues(3, i0, 0) = 0.0; // {2, 0, 0}
223  outputValues(3, i0, 1) = -(1.0 - z)/8.0; // {1, 1, 0}
224  outputValues(3, i0, 2) = (1.0 + y)/8.0; // {1, 0, 1}
225  outputValues(3, i0, 3) = 0.0; // {0, 2, 0}
226  outputValues(3, i0, 4) = -(1.0 - x)/8.0; // {0, 1, 1}
227  outputValues(3, i0, 5) = 0.0; // {0, 0, 2}
228 
229 
230  outputValues(4, i0, 0) = 0.0; // {2, 0, 0}
231  outputValues(4, i0, 1) = (1.0 + z)/8.0; // {1, 1, 0}
232  outputValues(4, i0, 2) = -(1.0 - y)/8.0; // {1, 0, 1}
233  outputValues(4, i0, 3) = 0.0; // {0, 2, 0}
234  outputValues(4, i0, 4) = -(1.0 - x)/8.0; // {0, 1, 1}
235  outputValues(4, i0, 5) = 0.0; // {0, 0, 2}
236 
237  outputValues(5, i0, 0) = 0.0; // {2, 0, 0}
238  outputValues(5, i0, 1) = -(1.0 + z)/8.0; // {1, 1, 0}
239  outputValues(5, i0, 2) = (1.0 - y)/8.0; // {1, 0, 1}
240  outputValues(5, i0, 3) = 0.0; // {0, 2, 0}
241  outputValues(5, i0, 4) = -(1.0 + x)/8.0; // {0, 1, 1}
242  outputValues(5, i0, 5) = 0.0; // {0, 0, 2}
243 
244  outputValues(6, i0, 0) = 0.0; // {2, 0, 0}
245  outputValues(6, i0, 1) = (1.0 + z)/8.0; // {1, 1, 0}
246  outputValues(6, i0, 2) = (1.0 + y)/8.0; // {1, 0, 1}
247  outputValues(6, i0, 3) = 0.0; // {0, 2, 0}
248  outputValues(6, i0, 4) = (1.0 + x)/8.0; // {0, 1, 1}
249  outputValues(6, i0, 5) = 0.0; // {0, 0, 2}
250 
251  outputValues(7, i0, 0) = 0.0; // {2, 0, 0}
252  outputValues(7, i0, 1) = -(1.0 + z)/8.0; // {1, 1, 0}
253  outputValues(7, i0, 2) = -(1.0 + y)/8.0; // {1, 0, 1}
254  outputValues(7, i0, 3) = 0.0; // {0, 2, 0}
255  outputValues(7, i0, 4) = (1.0 - x)/8.0; // {0, 1, 1}
256  outputValues(7, i0, 5) = 0.0; // {0, 0, 2}
257  }
258  break;
259 
260  case OPERATOR_D3:
261  case OPERATOR_D4:
262  case OPERATOR_D5:
263  case OPERATOR_D6:
264  case OPERATOR_D7:
265  case OPERATOR_D8:
266  case OPERATOR_D9:
267  case OPERATOR_D10:
268  {
269  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
270  int DkCardinality = Intrepid::getDkCardinality(operatorType,
271  this -> basisCellTopology_.getDimension() );
272  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
273  for (int i0 = 0; i0 < dim0; i0++) {
274  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
275  outputValues(dofOrd, i0, dkOrd) = 0.0;
276  }
277  }
278  }
279  }
280  break;
281 
282  default:
283  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
284  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): Invalid operator type");
285  }
286 }
287 
288 
289 
290 template<class Scalar, class ArrayScalar>
291 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
292  const ArrayScalar & inputPoints,
293  const ArrayScalar & cellVertices,
294  const EOperator operatorType) const {
295  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
296  ">>> ERROR (Basis_HGRAD_HEX_C1_FEM): FEM Basis calling an FVD member function");
297  }
298 
299 template<class Scalar, class ArrayScalar>
300 void Basis_HGRAD_HEX_C1_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const {
301 #ifdef HAVE_INTREPID_DEBUG
302  // Verify rank of output array.
303  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
304  ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords) rank = 2 required for DofCoords array");
305  // Verify 0th dimension of output array.
306  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
307  ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
308  // Verify 1st dimension of output array.
309  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
310  ">>> ERROR: (Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
311 #endif
312 
313  DofCoords(0,0) = -1.0; DofCoords(0,1) = -1.0; DofCoords(0,2) = -1.0;
314  DofCoords(1,0) = 1.0; DofCoords(1,1) = -1.0; DofCoords(1,2) = -1.0;
315  DofCoords(2,0) = 1.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = -1.0;
316  DofCoords(3,0) = -1.0; DofCoords(3,1) = 1.0; DofCoords(3,2) = -1.0;
317  DofCoords(4,0) = -1.0; DofCoords(4,1) = -1.0; DofCoords(4,2) = 1.0;
318  DofCoords(5,0) = 1.0; DofCoords(5,1) = -1.0; DofCoords(5,2) = 1.0;
319  DofCoords(6,0) = 1.0; DofCoords(6,1) = 1.0; DofCoords(6,2) = 1.0;
320  DofCoords(7,0) = -1.0; DofCoords(7,1) = 1.0; DofCoords(7,2) = 1.0;
321 }
322 
323 }// namespace Intrepid
324 
325 #endif
326 
Intrepid::Basis_HGRAD_HEX_C1_FEM::Basis_HGRAD_HEX_C1_FEM
Basis_HGRAD_HEX_C1_FEM()
Constructor.
Definition: Intrepid_HGRAD_HEX_C1_FEMDef.hpp:55
Intrepid::Basis_HGRAD_HEX_C1_FEM::initializeTags
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Definition: Intrepid_HGRAD_HEX_C1_FEMDef.hpp:68
Intrepid::Basis_HGRAD_HEX_C1_FEM::getDofCoords
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.
Definition: Intrepid_HGRAD_HEX_C1_FEMDef.hpp:300
Intrepid::Basis_HGRAD_HEX_C1_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
Definition: Intrepid_HGRAD_HEX_C1_FEMDef.hpp:100
Generated by   doxygen 1.8.5