103 const ArrayScalar & inputPoints,
104 const EOperator operatorType)
const {
107#ifdef HAVE_INTREPID_DEBUG
108 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
111 this -> getBaseCellTopology(),
112 this -> getCardinality() );
116 int dim0 = inputPoints.dimension(0);
127 Teuchos::Array<int> pt_tets;
130 switch (operatorType) {
133 for (
int i0 = 0; i0 < dim0; i0++) {
134 r = inputPoints(i0, 0);
135 s = inputPoints(i0, 1);
136 t = inputPoints(i0, 2);
138 pt_tets = getLocalSubTetrahedra(r,s,t);
144 if (pt_tets[0] != -1) {
145 int subtet = pt_tets[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;
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;
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;
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;
174 outputValues(4, i0) = 1. - 2. * (s + t);
175 outputValues(5, i0) = 2. * (r + s) - 1.;
176 outputValues(8, i0) = 2. * (r + t) - 1.;
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);
186 outputValues(7, i0) = 1. - 2. * (r + s);
187 outputValues(8, i0) = 2. * (r + t) - 1.;
188 outputValues(9, i0) = 2. * (s + t) - 1.;
192 outputValues(4, i0) = 1. - 2. * (s + t);
193 outputValues(7, i0) = 1. - 2. * (r + s);
194 outputValues(8, i0) = 2. * (r + t) - 1.;
198 outputValues(4, i0) = 1. - 2. * (s + t);
199 outputValues(5, i0) = 2. * (r + s) - 1.;
200 outputValues(6, i0) = 1. - 2. * (r + t);
204 outputValues(5, i0) = 2. * (r + s) - 1.;
205 outputValues(6, i0) = 1. - 2. * (r + t);
206 outputValues(9, i0) = 2. * (s + t) - 1.;
210 outputValues(6, i0) = 1. - 2. * (r + t);
211 outputValues(7, i0) = 1. - 2. * (r + s);
212 outputValues(9, i0) = 2. * (s + t) - 1.;
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.;
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;
236 outputValues.initialize(0.0);
238 FieldContainer<Scalar> Lopt(10,3);
239 for (
int pt=0; pt < dim0; ++pt) {
241 r = inputPoints(pt, 0);
242 s = inputPoints(pt, 1);
243 t = inputPoints(pt, 2);
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.;
252 Lopt(2,1) = -0.375 + (5*s)/2.;
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.;
276 for (
int node=0; node < 10; ++node) {
277 for (
int dim=0; dim < 3; ++dim) {
278 outputValues(node,pt,dim) = Lopt(node,dim);
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");
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");
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;
319 ">>> ERROR (Basis_HGRAD_TET_COMP12_FEM): Invalid operator type");
337#ifdef HAVE_INTREPID_DEBUG
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");
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");
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");
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;
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);
544 Scalar xj = det44(a);
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));
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));
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));
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));