Intrepid
Intrepid_HGRAD_TET_C2_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_TET_C2_FEMDEF_HPP
2 #define INTREPID_HGRAD_TET_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_ = 10;
57  this -> basisDegree_ = 2;
58  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<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
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  1, 0, 0, 1,
80  1, 1, 0, 1,
81  1, 2, 0, 1,
82  1, 3, 0, 1,
83  1, 4, 0, 1,
84  1, 5, 0, 1,
85  };
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,z) coordinates of the evaluation point
118  Scalar x = 0.0;
119  Scalar y = 0.0;
120  Scalar z = 0.0;
121 
122  switch (operatorType) {
123 
124  case OPERATOR_VALUE:
125  for (int i0 = 0; i0 < dim0; i0++) {
126  x = inputPoints(i0, 0);
127  y = inputPoints(i0, 1);
128  z = inputPoints(i0, 2);
129 
130  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
131  outputValues(0, i0) = (-1. + x + y + z)*(-1. + 2.*x + 2.*y + 2.*z);
132  outputValues(1, i0) = x*(-1. + 2.*x);
133  outputValues(2, i0) = y*(-1. + 2.*y);
134  outputValues(3, i0) = z*(-1. + 2.*z);
135 
136  outputValues(4, i0) = -4.*x*(-1. + x + y + z);
137  outputValues(5, i0) = 4.*x*y;
138  outputValues(6, i0) = -4.*y*(-1. + x + y + z);
139  outputValues(7, i0) = -4.*z*(-1. + x + y + z);
140  outputValues(8, i0) = 4.*x*z;
141  outputValues(9, i0) = 4.*y*z;
142 
143  }
144  break;
145 
146  case OPERATOR_GRAD:
147  case OPERATOR_D1:
148  for (int i0 = 0; i0 < dim0; i0++) {
149  x = inputPoints(i0, 0);
150  y = inputPoints(i0, 1);
151  z = inputPoints(i0, 2);
152 
153  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
154  outputValues(0, i0, 0) = -3.+ 4.*x + 4.*y + 4.*z;
155  outputValues(0, i0, 1) = -3.+ 4.*x + 4.*y + 4.*z;
156  outputValues(0, i0, 2) = -3.+ 4.*x + 4.*y + 4.*z;
157 
158  outputValues(1, i0, 0) = -1.+ 4.*x;
159  outputValues(1, i0, 1) = 0.;
160  outputValues(1, i0, 2) = 0.;
161 
162  outputValues(2, i0, 0) = 0.;
163  outputValues(2, i0, 1) = -1.+ 4.*y;
164  outputValues(2, i0, 2) = 0.;
165 
166  outputValues(3, i0, 0) = 0.;
167  outputValues(3, i0, 1) = 0.;
168  outputValues(3, i0, 2) = -1.+ 4.*z;
169 
170 
171  outputValues(4, i0, 0) = -4.*(-1.+ 2*x + y + z);
172  outputValues(4, i0, 1) = -4.*x;
173  outputValues(4, i0, 2) = -4.*x;
174 
175  outputValues(5, i0, 0) = 4.*y;
176  outputValues(5, i0, 1) = 4.*x;
177  outputValues(5, i0, 2) = 0.;
178 
179  outputValues(6, i0, 0) = -4.*y;
180  outputValues(6, i0, 1) = -4.*(-1.+ x + 2*y + z);
181  outputValues(6, i0, 2) = -4.*y;
182 
183  outputValues(7, i0, 0) = -4.*z;
184  outputValues(7, i0, 1) = -4.*z;
185  outputValues(7, i0, 2) = -4.*(-1.+ x + y + 2*z);
186 
187  outputValues(8, i0, 0) = 4.*z;
188  outputValues(8, i0, 1) = 0.;
189  outputValues(8, i0, 2) = 4.*x;
190 
191  outputValues(9, i0, 0) = 0.;
192  outputValues(9, i0, 1) = 4.*z;
193  outputValues(9, i0, 2) = 4.*y;
194 
195  }
196  break;
197 
198  case OPERATOR_CURL:
199  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
200  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
201  break;
202 
203  case OPERATOR_DIV:
204  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
205  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
206  break;
207 
208  case OPERATOR_D2:
209  for (int i0 = 0; i0 < dim0; i0++) {
210 
211  outputValues(0, i0, 0) = 4.;
212  outputValues(0, i0, 1) = 4.;
213  outputValues(0, i0, 2) = 4.;
214  outputValues(0, i0, 3) = 4.;
215  outputValues(0, i0, 4) = 4.;
216  outputValues(0, i0, 5) = 4.;
217 
218  outputValues(1, i0, 0) = 4.;
219  outputValues(1, i0, 1) = 0.;
220  outputValues(1, i0, 2) = 0.;
221  outputValues(1, i0, 3) = 0.;
222  outputValues(1, i0, 4) = 0.;
223  outputValues(1, i0, 5) = 0.;
224 
225  outputValues(2, i0, 0) = 0.;
226  outputValues(2, i0, 1) = 0.;
227  outputValues(2, i0, 2) = 0.;
228  outputValues(2, i0, 3) = 4.;
229  outputValues(2, i0, 4) = 0.;
230  outputValues(2, i0, 5) = 0.;
231 
232  outputValues(3, i0, 0) = 0.;
233  outputValues(3, i0, 1) = 0.;
234  outputValues(3, i0, 2) = 0.;
235  outputValues(3, i0, 3) = 0.;
236  outputValues(3, i0, 4) = 0.;
237  outputValues(3, i0, 5) = 4.;
238 
239  outputValues(4, i0, 0) = -8.;
240  outputValues(4, i0, 1) = -4.;
241  outputValues(4, i0, 2) = -4.;
242  outputValues(4, i0, 3) = 0.;
243  outputValues(4, i0, 4) = 0.;
244  outputValues(4, i0, 5) = 0.;
245 
246  outputValues(5, i0, 0) = 0.;
247  outputValues(5, i0, 1) = 4.;
248  outputValues(5, i0, 2) = 0.;
249  outputValues(5, i0, 3) = 0.;
250  outputValues(5, i0, 4) = 0.;
251  outputValues(5, i0, 5) = 0.;
252 
253  outputValues(6, i0, 0) = 0.;
254  outputValues(6, i0, 1) = -4.;
255  outputValues(6, i0, 2) = 0.;
256  outputValues(6, i0, 3) = -8.;
257  outputValues(6, i0, 4) = -4.;
258  outputValues(6, i0, 5) = 0;
259 
260  outputValues(7, i0, 0) = 0.;
261  outputValues(7, i0, 1) = 0.;
262  outputValues(7, i0, 2) = -4.;
263  outputValues(7, i0, 3) = 0.;
264  outputValues(7, i0, 4) = -4.;
265  outputValues(7, i0, 5) = -8.;
266 
267  outputValues(8, i0, 0) = 0.;
268  outputValues(8, i0, 1) = 0.;
269  outputValues(8, i0, 2) = 4.;
270  outputValues(8, i0, 3) = 0.;
271  outputValues(8, i0, 4) = 0.;
272  outputValues(8, i0, 5) = 0.;
273 
274  outputValues(9, i0, 0) = 0.;
275  outputValues(9, i0, 1) = 0.;
276  outputValues(9, i0, 2) = 0.;
277  outputValues(9, i0, 3) = 0.;
278  outputValues(9, i0, 4) = 4.;
279  outputValues(9, i0, 5) = 0.;
280  }
281  break;
282 
283  case OPERATOR_D3:
284  case OPERATOR_D4:
285  case OPERATOR_D5:
286  case OPERATOR_D6:
287  case OPERATOR_D7:
288  case OPERATOR_D8:
289  case OPERATOR_D9:
290  case OPERATOR_D10:
291  {
292  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
293  int DkCardinality = Intrepid::getDkCardinality(operatorType,
294  this -> basisCellTopology_.getDimension() );
295  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
296  for (int i0 = 0; i0 < dim0; i0++) {
297  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
298  outputValues(dofOrd, i0, dkOrd) = 0.0;
299  }
300  }
301  }
302  }
303  break;
304  default:
305  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
306  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): Invalid operator type");
307  }
308 }
309 
310 
311 
312 template<class Scalar, class ArrayScalar>
314  const ArrayScalar & inputPoints,
315  const ArrayScalar & cellVertices,
316  const EOperator operatorType) const {
317  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
318  ">>> ERROR (Basis_HGRAD_TET_C2_FEM): FEM Basis calling an FVD member function");
319 }
320 
321 
322 
323 template<class Scalar, class ArrayScalar>
325 #ifdef HAVE_INTREPID_DEBUG
326  // Verify rank of output array.
327  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
328  ">>> ERROR: (Intrepid::Basis_HGRAD_TET_C2_FEM::getDofCoords) rank = 2 required for DofCoords array");
329  // Verify 0th dimension of output array.
330  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
331  ">>> ERROR: (Intrepid::Basis_HGRAD_TET_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
332  // Verify 1st dimension of output array.
333  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
334  ">>> ERROR: (Intrepid::Basis_HGRAD_TET_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
335 #endif
336 
337  DofCoords(0,0) = 0.0; DofCoords(0,1) = 0.0; DofCoords(0,2) = 0.0;
338  DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0; DofCoords(1,2) = 0.0;
339  DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = 0.0;
340  DofCoords(3,0) = 0.0; DofCoords(3,1) = 0.0; DofCoords(3,2) = 1.0;
341 
342  DofCoords(4,0) = 0.5; DofCoords(4,1) = 0.0; DofCoords(4,2) = 0.0;
343  DofCoords(5,0) = 0.5; DofCoords(5,1) = 0.5; DofCoords(5,2) = 0.0;
344  DofCoords(6,0) = 0.0; DofCoords(6,1) = 0.5; DofCoords(6,2) = 0.0;
345  DofCoords(7,0) = 0.0; DofCoords(7,1) = 0.0; DofCoords(7,2) = 0.5;
346  DofCoords(8,0) = 0.5; DofCoords(8,1) = 0.0; DofCoords(8,2) = 0.5;
347  DofCoords(9,0) = 0.0; DofCoords(9,1) = 0.5; DofCoords(9,2) = 0.5;
348 }
349 
350 }// namespace Intrepid
351 #endif
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Tetrahedron.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.