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>
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>
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>
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>
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>
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 
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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...
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
Evaluation of a FEM basis on a reference Hexahedron cell.
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.