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 }// namespace Intrepid
321 #endif
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
FEM basis evaluation on a reference Tetrahedron cell.