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
52namespace 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
66template<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
101template<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
325template<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
335template<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
361template<class Scalar, class ArrayScalar>
362Teuchos::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
458template<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
489template<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
507template<class Scalar, class ArrayScalar>
508Scalar
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
539template<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
569template<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
668template<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
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
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.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Tetrahedron.
Intrepid::FieldContainer< Scalar > getSubTetDetF() const
Returns FieldContainer of local sub-tet detF.
Intrepid::FieldContainer< Scalar > getBarycentricCoords(const ArrayScalar &) const
Returns FieldContainer of Barycentric Coordinates for the input points.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Scalar det44(const Intrepid::FieldContainer< Scalar >) const
Returns FieldContainer of Barycentric Coordinates for the input points.
Intrepid::FieldContainer< Scalar > getWeights(const ArrayScalar &) const
Returns FieldContainer of local integration weights.
Intrepid::FieldContainer< Scalar > inverse44(const Intrepid::FieldContainer< Scalar >) const
Returns FieldContainer of Barycentric Coordinates for the input points.
Intrepid::FieldContainer< Scalar > getSubTetGrads() const
Returns FieldContainer of local sub-tet gradients.
Teuchos::Array< int > getLocalSubTetrahedra(Scalar, Scalar, Scalar) const
Returns array of local sub-tetrahdera that the given point resides in.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.