Intrepid
Intrepid_HGRAD_TET_COMP12_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_TET_COMP12_FEMDEF_HPP
2 #define INTREPID_HGRAD_TET_COMP12_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 
52 namespace Intrepid {
53 
54  template<class Scalar, class ArrayScalar>
56  {
57  this -> basisCardinality_ = 10;
58  this -> basisDegree_ = 1;
59  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<11> >() );
60  this -> basisType_ = BASIS_FEM_DEFAULT;
61  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62  this -> basisTagsAreSet_ = false;
63  }
64 
65 
66 template<class Scalar, class ArrayScalar>
68 
69  // Basis-dependent intializations
70  int tagSize = 4; // size of DoF tag
71  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
72  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
73  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
74 
75  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
76  int tags[] = { 0, 0, 0, 1,
77  0, 1, 0, 1,
78  0, 2, 0, 1,
79  0, 3, 0, 1,
80  1, 0, 0, 1,
81  1, 1, 0, 1,
82  1, 2, 0, 1,
83  1, 3, 0, 1,
84  1, 4, 0, 1,
85  1, 5, 0, 1,
86  };
87 
88  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
89  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
90  this -> ordinalToTag_,
91  tags,
92  this -> basisCardinality_,
93  tagSize,
94  posScDim,
95  posScOrd,
96  posDfOrd);
97 }
98 
99 
100 
101 template<class Scalar, class ArrayScalar>
103  const ArrayScalar & inputPoints,
104  const EOperator operatorType) const {
105 
106  // Verify arguments
107 #ifdef HAVE_INTREPID_DEBUG
108  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
109  inputPoints,
110  operatorType,
111  this -> getBaseCellTopology(),
112  this -> getCardinality() );
113 #endif
114 
115  // Number of evaluation points = dim 0 of inputPoints
116  int dim0 = inputPoints.dimension(0);
117 
118  // Temporaries: (x,y,z) coordinates of the evaluation point
119  Scalar r = 0.0;
120  Scalar s = 0.0;
121  Scalar t = 0.0;
122 
123  // Temporary for the auriliary node basis function
124  Scalar aux = 0.0;
125 
126  // Array to store all the subtets containing the given pt
127  Teuchos::Array<int> pt_tets;
128 
129 
130  switch (operatorType) {
131 
132  case OPERATOR_VALUE:
133  for (int i0 = 0; i0 < dim0; i0++) {
134  r = inputPoints(i0, 0);
135  s = inputPoints(i0, 1);
136  t = inputPoints(i0, 2);
137 
138  pt_tets = getLocalSubTetrahedra(r,s,t);
139 
140  // idependent verification shows that a given point will produce
141  // the same shape functions for each tet that contains it, so
142  // we only need to use the first one returned.
143  //for (int pt = 0; pt < pt_tets.size(); ++pt)
144  if (pt_tets[0] != -1) {
145  int subtet = pt_tets[0];
146  aux = 0.0;
147  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
148  switch (subtet) {
149  case 0:
150  outputValues(0, i0) = 1. - 2. * (r + s + t);
151  outputValues(4, i0) = 2. * r;
152  outputValues(6, i0) = 2. * s;
153  outputValues(7, i0) = 2. * t;
154  break;
155  case 1:
156  outputValues(1, i0) = 2. * r - 1.;
157  outputValues(4, i0) = 2. - 2. * (r + s + t);
158  outputValues(5, i0) = 2. * s;
159  outputValues(8, i0) = 2. * t;
160  break;
161  case 2:
162  outputValues(2, i0) = 2. * s - 1.;
163  outputValues(5, i0) = 2. * r;
164  outputValues(6, i0) = 2. - 2. * (r + s + t);
165  outputValues(9, i0) = 2. * t;
166  break;
167  case 3:
168  outputValues(3, i0) = 2. * t - 1.;
169  outputValues(7, i0) = 2. - 2. * (r + s + t);
170  outputValues(8, i0) = 2. * r;
171  outputValues(9, i0) = 2. * s;
172  break;
173  case 4:
174  outputValues(4, i0) = 1. - 2. * (s + t);
175  outputValues(5, i0) = 2. * (r + s) - 1.;
176  outputValues(8, i0) = 2. * (r + t) - 1.;
177  aux = 2. - 4. * r;
178  break;
179  case 5:
180  outputValues(5, i0) = 2. * (r + s) - 1.;
181  outputValues(8, i0) = 2. * (r + t) - 1.;
182  outputValues(9, i0) = 2. * (s + t) - 1.;
183  aux = 4. - 4. * (r + s + t);
184  break;
185  case 6:
186  outputValues(7, i0) = 1. - 2. * (r + s);
187  outputValues(8, i0) = 2. * (r + t) - 1.;
188  outputValues(9, i0) = 2. * (s + t) - 1.;
189  aux = 2. - 4. * t;
190  break;
191  case 7:
192  outputValues(4, i0) = 1. - 2. * (s + t);
193  outputValues(7, i0) = 1. - 2. * (r + s);
194  outputValues(8, i0) = 2. * (r + t) - 1.;
195  aux = 4. * s;
196  break;
197  case 8:
198  outputValues(4, i0) = 1. - 2. * (s + t);
199  outputValues(5, i0) = 2. * (r + s) - 1.;
200  outputValues(6, i0) = 1. - 2. * (r + t);
201  aux = 4. * t;
202  break;
203  case 9:
204  outputValues(5, i0) = 2. * (r + s) - 1.;
205  outputValues(6, i0) = 1. - 2. * (r + t);
206  outputValues(9, i0) = 2. * (s + t) - 1.;
207  aux = 2. - 4. * s;
208  break;
209  case 10:
210  outputValues(6, i0) = 1. - 2. * (r + t);
211  outputValues(7, i0) = 1. - 2. * (r + s);
212  outputValues(9, i0) = 2. * (s + t) - 1.;
213  aux = 4. * r;
214  break;
215  case 11:
216  outputValues(4, i0) = 1. - 2. * (s + t);
217  outputValues(6, i0) = 1. - 2. * (r + t);
218  outputValues(7, i0) = 1. - 2. * (r + s);
219  aux = 4. * (r + s + t) - 2.;
220  break;
221  }
222  outputValues(4, i0) += aux/6.0;
223  outputValues(5, i0) += aux/6.0;
224  outputValues(6, i0) += aux/6.0;
225  outputValues(7, i0) += aux/6.0;
226  outputValues(8, i0) += aux/6.0;
227  outputValues(9, i0) += aux/6.0;
228  }
229  }
230  break;
231 
232  case OPERATOR_GRAD:
233  case OPERATOR_D1:
234  {
235  // initialize to 0.0 since we will be accumulating
236  outputValues.initialize(0.0);
237 
238  FieldContainer<Scalar> Lopt(10,3);
239  for (int pt=0; pt < dim0; ++pt) {
240 
241  r = inputPoints(pt, 0);
242  s = inputPoints(pt, 1);
243  t = inputPoints(pt, 2);
244 
245  Lopt(0,0) = (-17 + 20*r + 20*s + 20*t)/8.;
246  Lopt(0,1) = (-17 + 20*r + 20*s + 20*t)/8.;
247  Lopt(0,2) = (-17 + 20*r + 20*s + 20*t)/8.;
248  Lopt(1,0) = -0.375 + (5*r)/2.;
249  Lopt(1,1) = 0.;
250  Lopt(1,2) = 0.;
251  Lopt(2,0) = 0.;
252  Lopt(2,1) = -0.375 + (5*s)/2.;
253  Lopt(2,2) = 0.;
254  Lopt(3,0) = 0.;
255  Lopt(3,1) = 0.;
256  Lopt(3,2) = -0.375 + (5*t)/2.;
257  Lopt(4,0) = (-35*(-1 + 2*r + s + t))/12.;
258  Lopt(4,1) = (-4 - 35*r + 5*s + 10*t)/12.;
259  Lopt(4,2) = (-4 - 35*r + 10*s + 5*t)/12.;
260  Lopt(5,0) = (-1 + 5*r + 40*s - 5*t)/12.;
261  Lopt(5,1) = (-1 + 40*r + 5*s - 5*t)/12.;
262  Lopt(5,2) = (-5*(-1 + r + s + 2*t))/12.;
263  Lopt(6,0) = (-4 + 5*r - 35*s + 10*t)/12.;
264  Lopt(6,1) = (-35*(-1 + r + 2*s + t))/12.;
265  Lopt(6,2) = (-4 + 10*r - 35*s + 5*t)/12.;
266  Lopt(7,0) = (-4 + 5*r + 10*s - 35*t)/12.;
267  Lopt(7,1) = (-4 + 10*r + 5*s - 35*t)/12.;
268  Lopt(7,2) = (-35*(-1 + r + s + 2*t))/12.;
269  Lopt(8,0) = (-1 + 5*r - 5*s + 40*t)/12.;
270  Lopt(8,1) = (-5*(-1 + r + 2*s + t))/12.;
271  Lopt(8,2) = (-1 + 40*r - 5*s + 5*t)/12.;
272  Lopt(9,0) = (-5*(-1 + 2*r + s + t))/12.;
273  Lopt(9,1) = (-1 - 5*r + 5*s + 40*t)/12.;
274  Lopt(9,2) = (-1 - 5*r + 40*s + 5*t)/12.;
275 
276  for (int node=0; node < 10; ++node) {
277  for (int dim=0; dim < 3; ++dim) {
278  outputValues(node,pt,dim) = Lopt(node,dim);
279  }
280  }
281  }
282  }
283  break;
284 
285  case OPERATOR_CURL:
286  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
287  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
288  break;
289 
290  case OPERATOR_DIV:
291  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
292  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
293  break;
294 
295  case OPERATOR_D2:
296  case OPERATOR_D3:
297  case OPERATOR_D4:
298  case OPERATOR_D5:
299  case OPERATOR_D6:
300  case OPERATOR_D7:
301  case OPERATOR_D8:
302  case OPERATOR_D9:
303  case OPERATOR_D10:
304  {
305  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
306  int DkCardinality = Intrepid::getDkCardinality(operatorType,
307  this -> basisCellTopology_.getDimension() );
308  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
309  for (int i0 = 0; i0 < dim0; i0++) {
310  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
311  outputValues(dofOrd, i0, dkOrd) = 0.0;
312  }
313  }
314  }
315  }
316  break;
317  default:
318  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
319  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Invalid operator type");
320  }
321 }
322 
323 
324 
325 template<class Scalar, class ArrayScalar>
327  const ArrayScalar & inputPoints,
328  const ArrayScalar & cellVertices,
329  const EOperator operatorType) const {
330  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
331  ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): FEM Basis calling an FVD member function");
332 
333 }
334 
335 template<class Scalar, class ArrayScalar>
337 #ifdef HAVE_INTREPID_DEBUG
338  // Verify rank of output array.
339  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
340  ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) rank = 2 required for DofCoords array");
341  // Verify 0th dimension of output array.
342  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
343  ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
344  // Verify 1st dimension of output array.
345  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
346  ">>> ERROR: (Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
347 #endif
348 
349  DofCoords(0,0) = 0.0; DofCoords(0,1) = 0.0; DofCoords(0,2) = 0.0;
350  DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0; DofCoords(1,2) = 0.0;
351  DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0; DofCoords(2,2) = 0.0;
352  DofCoords(3,0) = 0.0; DofCoords(3,1) = 0.0; DofCoords(3,2) = 1.0;
353  DofCoords(4,0) = 0.5; DofCoords(4,1) = 0.0; DofCoords(4,2) = 0.0;
354  DofCoords(5,0) = 0.5; DofCoords(5,1) = 0.5; DofCoords(5,2) = 0.0;
355  DofCoords(6,0) = 0.0; DofCoords(6,1) = 0.5; DofCoords(6,2) = 0.0;
356  DofCoords(7,0) = 0.0; DofCoords(7,1) = 0.0; DofCoords(7,2) = 0.5;
357  DofCoords(8,0) = 0.5; DofCoords(8,1) = 0.0; DofCoords(8,2) = 0.5;
358  DofCoords(9,0) = 0.0; DofCoords(9,1) = 0.5; DofCoords(9,2) = 0.5;
359 }
360 
361 template<class Scalar, class ArrayScalar>
362 Teuchos::Array<int>
364 {
365 
366  Teuchos::Array<int> subTets;
367  int count(0);
368 
369  // local coords
370  Scalar xyz = x + y + z;
371  Scalar xy = x + y;
372  Scalar xz = x + z;
373  Scalar yz = y + z;
374 
375  // cycle through each subdomain and push back if the point lies within
376 
377  // subtet #0 E0 := 0.0 <= r + s + t <= 0.5 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5
378  if ( (0.0 <= xyz && xyz <= 0.5) && (0.0 <= x && x <= 0.5) && (0.0 <= y && y <= 0.5) && (0.0 <= z && z <= 0.5) ) {
379  count++;
380  subTets.push_back(0);
381  }
382 
383  // subtet #1 E1 := 0.5 <= r + s + t <= 1.0 && 0.5 <= r <= 1.0 && 0.0 <= s <= 0.5 && 0.0 <= t <= 0.5
384  if ( (0.5 <= xyz && xyz <= 1.0) && (0.5 <= x && x <= 1.0) && (0.0 <= y && y <= 0.5) && (0.0 <= z && z <= 0.5) ) {
385  count++;
386  subTets.push_back(1);
387  }
388 
389  // subtet #2 E2 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.5 <= s <= 1.0 && 0.0 <= t <= 0.5
390  if ( (0.5 <= xyz && xyz <= 1.0) && (0.0 <= x && x <= 0.5) && (0.5 <= y && y <= 1.0) && (0.0 <= z && z<= 0.5) ) {
391  count++;
392  subTets.push_back(2);
393  }
394 
395  // subtet #3 E3 := 0.5 <= r + s + t <= 1.0 && 0.0 <= r <= 0.5 && 0.0 <= s <= 0.5 && 0.5 <= t <= 1.0
396  if ( (0.5 <= xyz && xyz <= 1.0) && (0.0 <= x && x <= 0.5) && (0.0 <= y && y <= 0.5) && (0.5 <= z && z <= 1.0) ) {
397  count++;
398  subTets.push_back(3);
399  }
400 
401  // subtet #4 E4 := 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= r + t <= 1.0 && 0.0 <= r <= 0.5
402  if ( (0.0 <= yz && yz <= 0.5) && (0.5 <= xy && xy <= 1.0) && (0.5 <= xz && xz <= 1.0) && (0.0 <= x && x <= 0.5) ) {
403  count++;
404  subTets.push_back(4);
405  }
406 
407  // subtet #5 E5 := 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.5 <= r + t <= 1.0 && 0.75 <= r + s + t <= 1.0
408  if ( (0.5 <= xy && xy <= 1.0) && (0.5 <= yz && yz <= 1.0) && (0.5 <= xz && xz <= 1.0) && (0.75 <= xyz && xyz <= 1.0) ) {
409  count++;
410  subTets.push_back(5);
411  }
412 
413  // subtet #6 E6 := 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= t <= 0.5
414  if ( (0.5 <= yz && yz <= 1.0) && (0.0 <= xy && xy <= 0.5) && (0.5 <= xz && xz <= 1.0) && (0.0 <= z && z <= 0.5) ) {
415  count++;
416  subTets.push_back(6);
417  }
418 
419  // subtet #7 E7 := 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5 && 0.5 <= r + t <= 1.0 && 0.0 <= s <= 0.25
420  if ( (0.0 <= yz && yz <= 0.5) && (0.0 <= xy && xy <= 0.5) && (0.5 <= xz && xz <= 1.0) && (0.0 <= y && y <= 0.25) ) {
421  count++;
422  subTets.push_back(7);
423  }
424 
425  // subtet #8 E8 := 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.0 <= t <= 0.25
426  if ( (0.0 <= xz && xz <= 0.5) && (0.0 <= yz && yz <= 0.5) && (0.5 <= xy && xy <= 1.0) && (0.0 <= z && z <= 0.25) ) {
427  count++;
428  subTets.push_back(8);
429  }
430 
431  // subtet #9 E9 := 0.0 <= r + t <= 0.5 && 0.5 <= r + s <= 1.0 && 0.5 <= s + t <= 1.0 && 0.0 <= s <= 0.5
432  if ( (0.0 <= xz && xz <= 0.5) && (0.5 <= xy && xy <= 1.0) && (0.5 <= yz && yz <= 1.0) && (0.0 <= y && y <= 0.5) ) {
433  count++;
434  subTets.push_back(9);
435  }
436 
437  // subtet #10 E10 := 0.0 <= r + t <= 0.5 && 0.5 <= s + t <= 1.0 && 0.0 <= r + s <= 0.5 && 0.0 <= r <= 0.25
438  if ( (0.0 <= xz && xz <= 0.5) && (0.5 <= yz && yz <= 1.0) && (0.0 <= xy && xy <= 0.5) && (0.0 <= x && x <= 0.25) ) {
439  count++;
440  subTets.push_back(10);
441  }
442 
443  // subtet #11 E11 := 0.5 <= r + s + t <= 0.75 && 0.0 <= r + t <= 0.5 && 0.0 <= s + t <= 0.5 && 0.0 <= r + s <= 0.5
444  if ( (0.5 <= xyz && xyz <= 0.75) && (0.0 <= xz && xz <= 0.5) && (0.0 <= yz && yz <= 0.5) && (0.0 <= xy && xy <= 0.5) ) {
445  count++;
446  subTets.push_back(11);
447  }
448 
449  // if the point doesn't lie in the parent domain return -1
450  if (count == 0) {
451  subTets.push_back(-1);
452  }
453 
454  return subTets;
455 }
456 
457 
458 template<class Scalar, class ArrayScalar>
461 {
462  int numPoints = inPts.dimension(0);
463  Intrepid::FieldContainer<Scalar> w(numPoints, 12);
464  w.initialize(0.0);
465  Teuchos::Array< Teuchos::Array<int> > pt_tets;
466 
467  for (int pt = 0; pt < numPoints; ++pt)
468  pt_tets.push_back(getLocalSubTetrahedra(inPts(pt,0), inPts(pt,1), inPts(pt,2)));
469 
470  Teuchos::Array<int> flat;
471  FieldContainer<int> count(12);
472 
473  for (int pt = 0; pt < numPoints; ++pt)
474  for (int i = 0; i < pt_tets[pt].size(); ++i)
475  flat.push_back(pt_tets[pt][i]);
476 
477  for (int i = 0; i < flat.size(); ++i)
478  count(flat[i])++;
479 
480  for (int pt = 0; pt < numPoints; ++pt)
481  for (int i = 0; i < pt_tets[pt].size(); ++i)
482  w(pt, pt_tets[pt][i]) = 1.0/count(pt_tets[pt][i]);
483 
484  return w;
485 }
486 
487 
488 
489 template<class Scalar, class ArrayScalar>
492 {
493  int numPoints = inPts.dimension(0);
494  Intrepid::FieldContainer<Scalar> lambda(numPoints, 4);
495 
496  for (int pt = 0; pt < numPoints; ++pt)
497  {
498  lambda(pt,0) = 1. - inPts(pt,0) - inPts(pt,1) - inPts(pt,2);
499  lambda(pt,1) = inPts(pt,0);
500  lambda(pt,2) = inPts(pt,1);
501  lambda(pt,3) = inPts(pt,2);
502  }
503 
504  return lambda;
505 }
506 
507 template<class Scalar, class ArrayScalar>
508 Scalar
510 {
511  Scalar det = a(0,3) * a(1,2) * a(2,1) * a(3,0)
512  - a(0,2) * a(1,3) * a(2,1) * a(3,0)
513  - a(0,3) * a(1,1) * a(2,2) * a(3,0)
514  + a(0,1) * a(1,3) * a(2,2) * a(3,0)
515  + a(0,2) * a(1,1) * a(2,3) * a(3,0)
516  - a(0,1) * a(1,2) * a(2,3) * a(3,0)
517  - a(0,3) * a(1,2) * a(2,0) * a(3,1)
518  + a(0,2) * a(1,3) * a(2,0) * a(3,1)
519  + a(0,3) * a(1,0) * a(2,2) * a(3,1)
520  - a(0,0) * a(1,3) * a(2,2) * a(3,1)
521  - a(0,2) * a(1,0) * a(2,3) * a(3,1)
522  + a(0,0) * a(1,2) * a(2,3) * a(3,1)
523  + a(0,3) * a(1,1) * a(2,0) * a(3,2)
524  - a(0,1) * a(1,3) * a(2,0) * a(3,2)
525  - a(0,3) * a(1,0) * a(2,1) * a(3,2)
526  + a(0,0) * a(1,3) * a(2,1) * a(3,2)
527  + a(0,1) * a(1,0) * a(2,3) * a(3,2)
528  - a(0,0) * a(1,1) * a(2,3) * a(3,2)
529  - a(0,2) * a(1,1) * a(2,0) * a(3,3)
530  + a(0,1) * a(1,2) * a(2,0) * a(3,3)
531  + a(0,2) * a(1,0) * a(2,1) * a(3,3)
532  - a(0,0) * a(1,2) * a(2,1) * a(3,3)
533  - a(0,1) * a(1,0) * a(2,2) * a(3,3)
534  + a(0,0) * a(1,1) * a(2,2) * a(3,3);
535 
536  return det;
537 }
538 
539 template<class Scalar, class ArrayScalar>
542 {
544  Scalar xj = det44(a);
545 
546  ai(0,0) = (1/xj) * (-a(1,3) * a(2,2) * a(3,1) + a(1,2) * a(2,3) * a(3,1) + a(1,3) * a(2,1) * a(3,2) - a(1,1) * a(2,3) * a(3,2) - a(1,2) * a(2,1) * a(3,3) + a(1,1) * a(2,2) * a(3,3));
547  ai(0,1) = (1/xj) * ( a(0,3) * a(2,2) * a(3,1) - a(0,2) * a(2,3) * a(3,1) - a(0,3) * a(2,1) * a(3,2) + a(0,1) * a(2,3) * a(3,2) + a(0,2) * a(2,1) * a(3,3) - a(0,1) * a(2,2) * a(3,3));
548  ai(0,2) = (1/xj) * (-a(0,3) * a(1,2) * a(3,1) + a(0,2) * a(1,3) * a(3,1) + a(0,3) * a(1,1) * a(3,2) - a(0,1) * a(1,3) * a(3,2) - a(0,2) * a(1,1) * a(3,3) + a(0,1) * a(1,2) * a(3,3));
549  ai(0,3) = (1/xj) * ( a(0,3) * a(1,2) * a(2,1) - a(0,2) * a(1,3) * a(2,1) - a(0,3) * a(1,1) * a(2,2) + a(0,1) * a(1,3) * a(2,2) + a(0,2) * a(1,1) * a(2,3) - a(0,1) * a(1,2) * a(2,3));
550 
551  ai(1,0) = (1/xj) * ( a(1,3) * a(2,2) * a(3,0) - a(1,2) * a(2,3) * a(3,0) - a(1,3) * a(2,0) * a(3,2) + a(1,0) * a(2,3) * a(3,2) + a(1,2) * a(2,0) * a(3,3) - a(1,0) * a(2,2) * a(3,3));
552  ai(1,1) = (1/xj) * (-a(0,3) * a(2,2) * a(3,0) + a(0,2) * a(2,3) * a(3,0) + a(0,3) * a(2,0) * a(3,2) - a(0,0) * a(2,3) * a(3,2) - a(0,2) * a(2,0) * a(3,3) + a(0,0) * a(2,2) * a(3,3));
553  ai(1,2) = (1/xj) * ( a(0,3) * a(1,2) * a(3,0) - a(0,2) * a(1,3) * a(3,0) - a(0,3) * a(1,0) * a(3,2) + a(0,0) * a(1,3) * a(3,2) + a(0,2) * a(1,0) * a(3,3) - a(0,0) * a(1,2) * a(3,3));
554  ai(1,3) = (1/xj) * (-a(0,3) * a(1,2) * a(2,0) + a(0,2) * a(1,3) * a(2,0) + a(0,3) * a(1,0) * a(2,2) - a(0,0) * a(1,3) * a(2,2) - a(0,2) * a(1,0) * a(2,3) + a(0,0) * a(1,2) * a(2,3));
555 
556  ai(2,0) = (1/xj) * (-a(1,3) * a(2,1) * a(3,0) + a(1,1) * a(2,3) * a(3,0) + a(1,3) * a(2,0) * a(3,1) - a(1,0) * a(2,3) * a(3,1) - a(1,1) * a(2,0) * a(3,3) + a(1,0) * a(2,1) * a(3,3));
557  ai(2,1) = (1/xj) * ( a(0,3) * a(2,1) * a(3,0) - a(0,1) * a(2,3) * a(3,0) - a(0,3) * a(2,0) * a(3,1) + a(0,0) * a(2,3) * a(3,1) + a(0,1) * a(2,0) * a(3,3) - a(0,0) * a(2,1) * a(3,3));
558  ai(2,2) = (1/xj) * (-a(0,3) * a(1,1) * a(3,0) + a(0,1) * a(1,3) * a(3,0) + a(0,3) * a(1,0) * a(3,1) - a(0,0) * a(1,3) * a(3,1) - a(0,1) * a(1,0) * a(3,3) + a(0,0) * a(1,1) * a(3,3));
559  ai(2,3) = (1/xj) * ( a(0,3) * a(1,1) * a(2,0) - a(0,1) * a(1,3) * a(2,0) - a(0,3) * a(1,0) * a(2,1) + a(0,0) * a(1,3) * a(2,1) + a(0,1) * a(1,0) * a(2,3) - a(0,0) * a(1,1) * a(2,3));
560 
561  ai(3,0) = (1/xj) * ( a(1,2) * a(2,1) * a(3,0) - a(1,1) * a(2,2) * a(3,0) - a(1,2) * a(2,0) * a(3,1) + a(1,0) * a(2,2) * a(3,1) + a(1,1) * a(2,0) * a(3,2) - a(1,0) * a(2,1) * a(3,2));
562  ai(3,1) = (1/xj) * (-a(0,2) * a(2,1) * a(3,0) + a(0,1) * a(2,2) * a(3,0) + a(0,2) * a(2,0) * a(3,1) - a(0,0) * a(2,2) * a(3,1) - a(0,1) * a(2,0) * a(3,2) + a(0,0) * a(2,1) * a(3,2));
563  ai(3,2) = (1/xj) * ( a(0,2) * a(1,1) * a(3,0) - a(0,1) * a(1,2) * a(3,0) - a(0,2) * a(1,0) * a(3,1) + a(0,0) * a(1,2) * a(3,1) + a(0,1) * a(1,0) * a(3,2) - a(0,0) * a(1,1) * a(3,2));
564  ai(3,3) = (1/xj) * (-a(0,2) * a(1,1) * a(2,0) + a(0,1) * a(1,2) * a(2,0) + a(0,2) * a(1,0) * a(2,1) - a(0,0) * a(1,2) * a(2,1) - a(0,1) * a(1,0) * a(2,2) + a(0,0) * a(1,1) * a(2,2));
565 
566  return ai;
567 }
568 
569 template<class Scalar, class ArrayScalar>
572 {
574  dx.initialize(0.0);
575  // fill in dx
576  dx(0,0,0) = -2.0;
577  dx(0,1,1) = 2.0;
578  dx(0,4,0) = 2.0;
579  dx(0,4,1) = -2.0;
580  dx(0,5,2) = 2.0;
581  dx(0,5,4) = 2.0;
582  dx(0,5,5) = 2.0;
583  dx(0,5,8) = 2.0;
584  dx(0,5,9) = 2.0;
585  dx(0,6,2) = -2.0;
586  dx(0,6,8) = -2.0;
587  dx(0,6,9) = -2.0;
588  dx(0,6,10) = -2.0;
589  dx(0,6,11) = -2.0;
590  dx(0,7,3) = -2.0;
591  dx(0,7,6) = -2.0;
592  dx(0,7,7) = -2.0;
593  dx(0,7,10) = -2.0;
594  dx(0,7,11) = -2.0;
595  dx(0,8,3) = 2.0;
596  dx(0,8,4) = 2.0;
597  dx(0,8,5) = 2.0;
598  dx(0,8,6) = 2.0;
599  dx(0,8,7) = 2.0;
600  dx(0,10,4) = -4.0;
601  dx(0,10,5) = -4.0;
602  dx(0,10,10) = 4.0;
603  dx(0,10,11) = 4.0;
604 
605  // fill in dy
606  dx(1,0,0) = -2.0;
607  dx(1,2,2) = 2.0;
608  dx(1,4,1) = -2.0;
609  dx(1,4,4) = -2.0;
610  dx(1,4,7) = -2.0;
611  dx(1,4,8) = -2.0;
612  dx(1,4,11) = -2.0;
613  dx(1,5,1) = 2.0;
614  dx(1,5,4) = 2.0;
615  dx(1,5,5) = 2.0;
616  dx(1,5,8) = 2.0;
617  dx(1,5,9) = 2.0;
618  dx(1,6,0) = 2.0;
619  dx(1,6,2) = -2.0;
620  dx(1,7,3) = -2.0;
621  dx(1,7,6) = -2.0;
622  dx(1,7,7) = -2.0;
623  dx(1,7,10) = -2.0;
624  dx(1,7,11) = -2.0;
625  dx(1,9,3) = 2.0;
626  dx(1,9,5) = 2.0;
627  dx(1,9,6) = 2.0;
628  dx(1,9,9) = 2.0;
629  dx(1,9,10) = 2.0;
630  dx(1,10,5) = -4.0;
631  dx(1,10,7) = 4.0;
632  dx(1,10,9) = -4.0;
633  dx(1,10,11) = 4.0;
634 
635  // fill in dz
636  dx(2,0,0) = -2.0;
637  dx(2,3,3) = 2.0;
638  dx(2,4,1) = -2.0;
639  dx(2,4,4) = -2.0;
640  dx(2,4,7) = -2.0;
641  dx(2,4,8) = -2.0;
642  dx(2,4,11) = -2.0;
643  dx(2,6,2) = -2.0;
644  dx(2,6,8) = -2.0;
645  dx(2,6,9) = -2.0;
646  dx(2,6,10) = -2.0;
647  dx(2,6,11) = -2.0;
648  dx(2,7,0) = 2.0;
649  dx(2,7,3) = -2.0;
650  dx(2,8,1) = 2.0;
651  dx(2,8,4) = 2.0;
652  dx(2,8,5) = 2.0;
653  dx(2,8,6) = 2.0;
654  dx(2,8,7) = 2.0;
655  dx(2,9,2) = 2.0;
656  dx(2,9,5) = 2.0;
657  dx(2,9,6) = 2.0;
658  dx(2,9,9) = 2.0;
659  dx(2,9,10) = 2.0;
660  dx(2,10,5) = -4.0;
661  dx(2,10,6) = -4.0;
662  dx(2,10,8) = 4.0;
663  dx(2,10,11) = 4.0;
664 
665  return dx;
666 }
667 
668 template<class Scalar, class ArrayScalar>
671 {
673  // set sub-elem jacobians
674  xJ(0) = 1./48.; xJ(1) = 1./48.; xJ(2) = 1./48.; xJ(3) = 1./48.;
675  xJ(4) = 1./96.; xJ(5) = 1./96.; xJ(6) = 1./96.; xJ(7) = 1./96.;
676  xJ(8) = 1./96.; xJ(9) = 1./96.; xJ(10) = 1./96.; xJ(11) = 1./96.;
677 
678  return xJ;
679 }
680 
681 
682 }// namespace Intrepid
683 #endif
Scalar det44(const Intrepid::FieldContainer< Scalar >) const
Returns FieldContainer of Barycentric Coordinates for the input points.
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.
Intrepid::FieldContainer< Scalar > getSubTetGrads() const
Returns FieldContainer of local sub-tet gradients.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
Intrepid::FieldContainer< Scalar > getBarycentricCoords(const ArrayScalar &) const
Returns FieldContainer of Barycentric Coordinates for the input points.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Intrepid::FieldContainer< Scalar > getWeights(const ArrayScalar &) const
Returns FieldContainer of local integration weights.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
Teuchos::Array< int > getLocalSubTetrahedra(Scalar, Scalar, Scalar) const
Returns array of local sub-tetrahdera that the given point resides in.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
Intrepid::FieldContainer< Scalar > inverse44(const Intrepid::FieldContainer< Scalar >) const
Returns FieldContainer of Barycentric Coordinates for the input points.
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.
Intrepid::FieldContainer< Scalar > getSubTetDetF() const
Returns FieldContainer of local sub-tet detF.
int dimension(const int whichDim) const
Returns the specified dimension.