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>
55  Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::Basis_HGRAD_TET_COMP12_FEM()
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>
67 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::initializeTags() {
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>
102 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues,
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>
326 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
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>
336 void Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const {
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>
363 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getLocalSubTetrahedra(Scalar x, Scalar y, Scalar z) const
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>
459 Intrepid::FieldContainer<Scalar>
460 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getWeights(const ArrayScalar & inPts) const
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>
490 Intrepid::FieldContainer<Scalar>
491 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getBarycentricCoords(const ArrayScalar & inPts) const
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
509 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::det44(const Intrepid::FieldContainer<Scalar> a) const
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>
540 Intrepid::FieldContainer<Scalar>
541 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::inverse44(const Intrepid::FieldContainer<Scalar> a) const
542 {
543  Intrepid::FieldContainer<Scalar> ai(4,4);
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>
570 Intrepid::FieldContainer<Scalar>
571 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getSubTetGrads() const
572 {
573  Intrepid::FieldContainer<Scalar> dx(3,11,12);
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>
669 Intrepid::FieldContainer<Scalar>
670 Basis_HGRAD_TET_COMP12_FEM<Scalar, ArrayScalar>::getSubTetDetF() const
671 {
672  Intrepid::FieldContainer<Scalar> xJ(12);
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
Intrepid::FieldContainer::dimension
int dimension(const int whichDim) const
Returns the specified dimension.
Definition: Intrepid_FieldContainerDef.hpp:658
Intrepid::Basis_HGRAD_TET_COMP12_FEM::inverse44
Intrepid::FieldContainer< Scalar > inverse44(const Intrepid::FieldContainer< Scalar >) const
Returns FieldContainer of Barycentric Coordinates for the input points.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:541
Intrepid::Basis_HGRAD_TET_COMP12_FEM::initializeTags
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:67
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getSubTetDetF
Intrepid::FieldContainer< Scalar > getSubTetDetF() const
Returns FieldContainer of local sub-tet detF.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:670
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:102
Intrepid::Basis_HGRAD_TET_COMP12_FEM::Basis_HGRAD_TET_COMP12_FEM
Basis_HGRAD_TET_COMP12_FEM()
Constructor.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:55
Intrepid::FieldContainer::initialize
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
Definition: Intrepid_FieldContainerDef.hpp:1374
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getSubTetGrads
Intrepid::FieldContainer< Scalar > getSubTetGrads() const
Returns FieldContainer of local sub-tet gradients.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:571
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getLocalSubTetrahedra
Teuchos::Array< int > getLocalSubTetrahedra(Scalar, Scalar, Scalar) const
Returns array of local sub-tetrahdera that the given point resides in.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:363
Intrepid::Basis_HGRAD_TET_COMP12_FEM::det44
Scalar det44(const Intrepid::FieldContainer< Scalar >) const
Returns FieldContainer of Barycentric Coordinates for the input points.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:509
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getWeights
Intrepid::FieldContainer< Scalar > getWeights(const ArrayScalar &) const
Returns FieldContainer of local integration weights.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:460
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getDofCoords
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Tetrahedron.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:336
Intrepid::Basis_HGRAD_TET_COMP12_FEM::getBarycentricCoords
Intrepid::FieldContainer< Scalar > getBarycentricCoords(const ArrayScalar &) const
Returns FieldContainer of Barycentric Coordinates for the input points.
Definition: Intrepid_HGRAD_TET_COMP12_FEMDef.hpp:491
Generated by   doxygen 1.8.5