00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceLagrange.hpp"
00032 #include "SundanceADReal.hpp"
00033 #include "PlayaExceptions.hpp"
00034 #include "SundanceSpatialDerivSpecifier.hpp"
00035 #include "SundancePoint.hpp"
00036 #include "SundanceObjectWithVerbosity.hpp"
00037 #include "SundanceOut.hpp"
00038
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043
00044
00045 Lagrange::Lagrange(int order)
00046 : order_(order) , doFInfromationCalculated_(false)
00047 {
00048 TEUCHOS_TEST_FOR_EXCEPTION(order < 0, std::runtime_error,
00049 "invalid polynomial order=" << order
00050 << " in Lagrange ctor");
00051 }
00052
00053 bool Lagrange::supportsCellTypePair(
00054 const CellType& maximalCellType,
00055 const CellType& cellType
00056 ) const
00057 {
00058 switch(maximalCellType)
00059 {
00060 case LineCell:
00061 switch(cellType)
00062 {
00063 case LineCell:
00064 case PointCell:
00065 return true;
00066 default:
00067 return false;
00068 }
00069 case TriangleCell:
00070 switch(cellType)
00071 {
00072 case TriangleCell:
00073 case LineCell:
00074 case PointCell:
00075 return true;
00076 default:
00077 return false;
00078 }
00079 case QuadCell:
00080 switch(cellType){
00081 case QuadCell:
00082 case LineCell:
00083 case PointCell:
00084 return true;
00085 default:
00086 return false;
00087 }
00088 case TetCell:
00089 switch(cellType)
00090 {
00091 case TetCell:
00092 case TriangleCell:
00093 case LineCell:
00094 case PointCell:
00095 return true;
00096 default:
00097 return false;
00098 }
00099 case BrickCell:
00100 switch(cellType)
00101 {
00102 case BrickCell:
00103 case QuadCell:
00104 case LineCell:
00105 case PointCell:
00106 return true;
00107 default:
00108 return false;
00109 }
00110 default:
00111 return false;
00112 }
00113 }
00114
00115 void Lagrange::print(std::ostream& os) const
00116 {
00117 os << "Lagrange(" << order_ << ")";
00118 }
00119
00120 int Lagrange::nReferenceDOFsWithoutFacets(
00121 const CellType& maximalCellType,
00122 const CellType& cellType
00123 ) const
00124 {
00125 if (order_==0)
00126 {
00127 if (maximalCellType != cellType) return 0;
00128 return 1;
00129 }
00130
00131 switch(cellType)
00132 {
00133 case PointCell:
00134 return 1;
00135 case LineCell:
00136 return order_-1;
00137 case TriangleCell:
00138 if (order_ < 3) return 0;
00139 if (order_ == 3) return 1;
00140 TEUCHOS_TEST_FOR_EXCEPTION(order_>3, std::runtime_error,
00141 "Lagrange order > 3 not implemented "
00142 "for triangle cells");
00143 return 0;
00144 case QuadCell:
00145 if (order_==1) return 0;
00146 if (order_==2) return 1;
00147 TEUCHOS_TEST_FOR_EXCEPTION(order_>2, std::runtime_error,
00148 "Lagrange order > 2 not implemented "
00149 "for quad cells");
00150 case TetCell:
00151 if (order_<=2) return 0;
00152 TEUCHOS_TEST_FOR_EXCEPTION(order_>2, std::runtime_error,
00153 "Lagrange order > 2 not implemented "
00154 "for tet cells");
00155 case BrickCell:
00156 if (order_<=1) return 0;
00157 if (order_==2) return 1;
00158 TEUCHOS_TEST_FOR_EXCEPTION(order_>2, std::runtime_error,
00159 "Lagrange order > 2 not implemented "
00160 "for brick cells");
00161 default:
00162 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00163 << cellType << " not implemented in Lagrange basis");
00164 return -1;
00165 }
00166 }
00167
00168 void Lagrange::getReferenceDOFs(
00169 const CellType& maximalCellType,
00170 const CellType& cellType,
00171 Array<Array<Array<int> > >& dofs) const
00172 {
00173 typedef Array<int> Aint;
00174
00175 if (order_==0)
00176 {
00177 int dim = dimension(cellType);
00178 dofs.resize(dim+1);
00179 for (int d=0; d<dim; d++)
00180 {
00181 Array<Array<int> > dd;
00182 for (int f=0; f<numFacets(cellType, d); f++) dd.append(Array<int>());
00183 dofs[d] = dd;
00184 }
00185 if (cellType==maximalCellType) dofs[dim] = tuple<Aint>(tuple(0));
00186 else dofs[dim] = tuple<Aint>(Array<int>());
00187 return;
00188 }
00189
00190 switch(cellType)
00191 {
00192 case PointCell:
00193 dofs.resize(1);
00194 dofs[0] = tuple<Aint>(tuple(0));
00195 return;
00196 case LineCell:
00197 dofs.resize(2);
00198 dofs[0] = tuple<Aint>(tuple(0), tuple(1));
00199 dofs[1] = tuple<Aint>(makeRange(2, order()));
00200 return;
00201 case TriangleCell:
00202 {
00203 int n = order()-1;
00204 dofs.resize(3);
00205 dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2));
00206 dofs[1] = tuple<Aint>(makeRange(3,2+n),
00207 makeRange(3+n, 2+2*n),
00208 makeRange(3+2*n, 2+3*n));
00209 if (order()<3)
00210 {
00211 dofs[2] = tuple(Array<int>());
00212 }
00213 else
00214 {
00215 dofs[2] = tuple<Aint>(makeRange(3+3*n, 3+3*n));
00216 }
00217 return;
00218 }
00219 case QuadCell:{
00220 dofs.resize(3);
00221
00222 dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2), tuple(3));
00223 if (order() == 2)
00224 {
00225
00226 dofs[1] = tuple<Aint>(tuple(4), tuple(5), tuple(6), tuple(7));
00227
00228
00229
00230 dofs[2] = tuple<Aint>(tuple(8));
00231 } else {
00232 dofs[1] = tuple(Array<int>(), Array<int>(),
00233 Array<int>(), Array<int>());
00234 dofs[2] = tuple(Array<int>());
00235 }
00236 return;
00237 }
00238 case TetCell:
00239 {
00240 dofs.resize(4);
00241 dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2), tuple(3));
00242 if (order() == 2)
00243 {
00244 dofs[1] = tuple<Aint>(tuple(4), tuple(5), tuple(6),
00245 tuple(7), tuple(8), tuple(9));
00246 }
00247 else
00248 {
00249 dofs[1] = tuple(Array<int>(), Array<int>(),
00250 Array<int>(), Array<int>(),
00251 Array<int>(), Array<int>());
00252 }
00253 dofs[2] = tuple(Array<int>(), Array<int>(),
00254 Array<int>(), Array<int>());
00255 dofs[3] = tuple(Array<int>());
00256 return;
00257 }
00258 case BrickCell:{
00259 dofs.resize(4);
00260
00261 dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2), tuple(3),
00262 tuple(4), tuple(5), tuple(6), tuple(7));
00263 if (order() == 2)
00264 {
00265
00266 dofs[1] = tuple<Aint>( tuple(8), tuple(9) , tuple(10), tuple(11),
00267 tuple(12), tuple(13), tuple(14), tuple(15),
00268 tuple(16), tuple(17), tuple(18), tuple(19));
00269
00270 dofs[2] = tuple<Aint>(tuple(20), tuple(21), tuple(22), tuple(23),
00271 tuple(24), tuple(25) );
00272 dofs[3] = tuple<Aint>(tuple(26));
00273 } else {
00274 dofs[1] = tuple(Array<int>(), Array<int>(), Array<int>(), Array<int>(),
00275 Array<int>(), Array<int>(), Array<int>(), Array<int>(),
00276 Array<int>(), Array<int>(), Array<int>(), Array<int>());
00277 dofs[2] = tuple(Array<int>(), Array<int>(), Array<int>(), Array<int>(),
00278 Array<int>(), Array<int>() );
00279 dofs[3] = tuple(Array<int>());
00280 }
00281 return;
00282 }
00283 default:
00284 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00285 << cellType << " not implemented in Lagrange basis");
00286 }
00287 }
00288
00289
00290
00291 Array<int> Lagrange::makeRange(int low, int high)
00292 {
00293 if (high < low) return Array<int>();
00294
00295 Array<int> rtn(high-low+1);
00296 for (int i=0; i<rtn.length(); i++) rtn[i] = low+i;
00297 return rtn;
00298 }
00299
00300 void Lagrange::refEval(
00301 const CellType& cellType,
00302 const Array<Point>& pts,
00303 const SpatialDerivSpecifier& sds,
00304 Array<Array<Array<double> > >& result,
00305 int verbosity) const
00306 {
00307 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00308 std::runtime_error,
00309 "cannot evaluate spatial derivative " << sds << " on Lagrange basis");
00310 const MultiIndex& deriv = sds.mi();
00311 typedef Array<double> Adouble;
00312 result.resize(1);
00313 result[0].resize(pts.length());
00314
00315 switch(cellType)
00316 {
00317 case PointCell:
00318 result[0] = tuple<Adouble>(tuple(1.0));
00319 return;
00320 case LineCell:
00321 for (int i=0; i<pts.length(); i++)
00322 {
00323 evalOnLine(pts[i], deriv, result[0][i]);
00324 }
00325 return;
00326 case TriangleCell:
00327 for (int i=0; i<pts.length(); i++)
00328 {
00329 evalOnTriangle(pts[i], deriv, result[0][i]);
00330 }
00331 return;
00332 case QuadCell:
00333 for (int i=0; i<pts.length(); i++)
00334 {
00335 evalOnquad(pts[i], deriv, result[0][i]);
00336 }
00337 return;
00338 case TetCell:
00339 for (int i=0; i<pts.length(); i++)
00340 {
00341 evalOnTet(pts[i], deriv, result[0][i]);
00342 }
00343 return;
00344 case BrickCell:
00345 for (int i=0; i<pts.length(); i++)
00346 {
00347 evalOnBrick(pts[i], deriv, result[0][i]);
00348 }
00349 return;
00350 default:
00351 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00352 "Lagrange::refEval() unimplemented for cell type "
00353 << cellType);
00354
00355 }
00356 }
00357
00358
00359
00360 void Lagrange::evalOnLine(const Point& pt,
00361 const MultiIndex& deriv,
00362 Array<double>& result) const
00363 {
00364 ADReal x = ADReal(pt[0], 0, 1);
00365 ADReal one(1.0, 1);
00366
00367 result.resize(order()+1);
00368 Array<ADReal> tmp(result.length());
00369 Array<double> x0(order()+1);
00370
00371 if (order_ == 0)
00372 {
00373 tmp[0] = one;
00374 }
00375 else
00376 {
00377 x0[0] = 0.0;
00378 x0[1] = 1.0;
00379 for (int i=0; i<order_-1; i++)
00380 {
00381 x0[i+2] = (i+1.0)/order_;
00382 }
00383
00384 for (int i=0; i<=order_; i++)
00385 {
00386 tmp[i] = one;
00387 for (int j=0; j<=order_; j++)
00388 {
00389 if (i==j) continue;
00390 tmp[i] *= (x - x0[j])/(x0[i]-x0[j]);
00391 }
00392 }
00393 }
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421 for (int i=0; i<tmp.length(); i++)
00422 {
00423 if (deriv.order()==0) result[i] = tmp[i].value();
00424 else result[i] = tmp[i].gradient()[0];
00425 }
00426 }
00427
00428 void Lagrange::evalOnTriangle(const Point& pt,
00429 const MultiIndex& deriv,
00430 Array<double>& result) const
00431
00432
00433
00434 {
00435 ADReal x = ADReal(pt[0], 0, 2);
00436 ADReal y = ADReal(pt[1], 1, 2);
00437 ADReal one(1.0, 2);
00438
00439 Array<ADReal> tmp;
00440
00441 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00442 << y.value());
00443
00444 switch(order_)
00445 {
00446 case 0:
00447 result.resize(1);
00448 tmp.resize(1);
00449 tmp[0] = one;
00450 break;
00451 case 1:
00452 result.resize(3);
00453 tmp.resize(3);
00454 tmp[0] = 1.0 - x - y;
00455 tmp[1] = x;
00456 tmp[2] = y;
00457 break;
00458 case 2:
00459 result.resize(6);
00460 tmp.resize(6);
00461 tmp[0] = (1.0-x-y)*(1.0-2.0*x-2.0*y);
00462 tmp[1] = 2.0*x*(x-0.5);
00463 tmp[2] = 2.0*y*(y-0.5);
00464 tmp[5] = 4.0*x*(1.0-x-y);
00465 tmp[3] = 4.0*x*y;
00466 tmp[4] = 4.0*y*(1.0-x-y);
00467 break;
00468 case 3:
00469 result.resize(10);
00470 tmp.resize(10);
00471 {
00472 ADReal q1 = 1.0 - x - y;
00473 ADReal q2 = x;
00474 ADReal q3 = y;
00475 tmp[0] = 9.0/2.0 * q1 * (q1 - 2.0/3.0) * (q1 - 1.0/3.0);
00476 tmp[1] = 9.0/2.0 * q2 * (q2 - 2.0/3.0) * (q2 - 1.0/3.0);
00477 tmp[2] = 9.0/2.0 * q3 * (q3 - 2.0/3.0) * (q3 - 1.0/3.0);
00478 tmp[7] = 27.0/2.0 * q1 * q2 * (q1 - 1.0/3.0);
00479 tmp[8] = 27.0/2.0 * q1 * q2 * (q2 - 1.0/3.0);
00480 tmp[3] = 27.0/2.0 * q2 * q3 * (q2 - 1.0/3.0);
00481 tmp[4] = 27.0/2.0 * q2 * q3 * (q3 - 1.0/3.0);
00482 tmp[6] = 27.0/2.0 * q3 * q1 * (q3 - 1.0/3.0);
00483 tmp[5] = 27.0/2.0 * q3 * q1 * (q1 - 1.0/3.0);
00484 tmp[9] = 27.0 * q1 * q2 * q3;
00485 }
00486 break;
00487 default:
00488 #ifndef TRILINOS_7
00489 SUNDANCE_ERROR("Lagrange::evalOnTriangle poly order > 2");
00490 #else
00491 SUNDANCE_ERROR7("Lagrange::evalOnTriangle poly order > 2");
00492 #endif
00493 }
00494
00495 for (int i=0; i<tmp.length(); i++)
00496 {
00497 SUNDANCE_OUT(this->verb() > 3,
00498 "tmp[" << i << "]=" << tmp[i].value()
00499 << " grad=" << tmp[i].gradient());
00500 if (deriv.order()==0) result[i] = tmp[i].value();
00501 else
00502 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00503 }
00504 }
00505
00506 void Lagrange::evalOnquad(const Point& pt,
00507 const MultiIndex& deriv,
00508 Array<double>& result) const {
00509
00510 ADReal x = ADReal(pt[0], 0, 2);
00511 ADReal y = ADReal(pt[1], 1, 2);
00512 ADReal one(1.0, 2);
00513
00514 Array<ADReal> tmp;
00515
00516 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00517 << y.value());
00518
00519 switch(order_)
00520 {
00521 case 0:
00522 result.resize(1);
00523 tmp.resize(1);
00524 tmp[0] = one;
00525 break;
00526 case 1:
00527 result.resize(4);
00528 tmp.resize(4);
00529 tmp[0] = (1.0 - x)*(1.0 - y);
00530 tmp[1] = x*(1.0 - y);
00531 tmp[2] = (1.0 - x)*y;
00532 tmp[3] = x*y;
00533 break;
00534 case 2:
00535
00536 result.resize(9);
00537 tmp.resize(9);
00538 tmp[0] = 9.0 - 9.0*x - 9.0*y + 9.0*x*y - 6.0*(x*x -1)*(y-1)
00539 - 6.0*(x-1)*(y*y-1) + 4.0*(x*x-1)*(y*y-1);
00540 tmp[1] = 6.0 - 3.0*x - 6.0*y + 3.0*x*y - 6.0*(x*x -1)*(y-1)
00541 - 3.0*(x-1)*(y*y-1) + 1.0*(x+1)*(y*y-1) + 4.0*(x*x-1)*(y*y-1);
00542 tmp[2] = 6.0 - 6.0*x - 3.0*y + 3.0*x*y - 3.0*(x*x -1)*(y-1) + 1.0*(x*x -1)*(y+1)
00543 -6.0*(x-1)*(y*y-1) + + 4.0*(x*x-1)*(y*y-1);
00544 tmp[3] = 4.0 - 2.0*x - 2.0*y + 1.0*x*y - 3.0*(x*x -1)*(y-1) + 1.0*(x*x -1)*(y+1)
00545 - 3.0*(x-1)*(y*y-1) + 1.0*(x+1)*(y*y-1) + 4.0*(x*x-1)*(y*y-1);
00546 tmp[4] = -12.0 + 12.0*x + 12.0*y - 12.0*x*y + 12.0*(x*x -1)*(y-1)
00547 + 8.0*(x-1)*(y*y-1) - 8.0*(x*x-1)*(y*y-1);
00548 tmp[5] = -12.0 + 12.0*x + 12.0*y - 12.0*x*y + 8.0*(x*x -1)*(y-1)
00549 + 12.0*(x-1)*(y*y-1) - 8.0*(x*x-1)*(y*y-1);
00550 tmp[6] = -8.0 + 4.0*x + 8.0*y - 4.0*x*y + 8.0*(x*x -1)*(y-1)
00551 + 6.0*(x-1)*(y*y-1) - 2.0*(x+1)*(y*y-1) - 8.0*(x*x-1)*(y*y-1);
00552 tmp[7] = -8.0 + 8.0*x + 4.0*y - 4.0*x*y + 6.0*(x*x -1)*(y-1) - 2.0*(x*x -1)*(y+1)
00553 + 8.0*(x-1)*(y*y-1) - 8.0*(x*x-1)*(y*y-1);
00554 tmp[8] = 16.0 - 16.0*x - 16.0*y + 16.0*x*y - 16.0*(x*x -1)*(y-1)
00555 - 16.0*(x-1)*(y*y-1) + 16.0*(x*x-1)*(y*y-1);
00556 break;
00557 default:{}
00558 }
00559
00560 for (int i=0; i<tmp.length(); i++) {
00561 SUNDANCE_OUT(this->verb() > 3 ,
00562 "tmp[" << i << "]=" << tmp[i].value()
00563 << " grad=" << tmp[i].gradient());
00564 if (deriv.order()==0)
00565 result[i] = tmp[i].value();
00566 else
00567 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00568 }
00569 }
00570
00571 void Lagrange::evalOnTet(const Point& pt,
00572 const MultiIndex& deriv,
00573 Array<double>& result) const
00574 {
00575 ADReal x = ADReal(pt[0], 0, 3);
00576 ADReal y = ADReal(pt[1], 1, 3);
00577 ADReal z = ADReal(pt[2], 2, 3);
00578 ADReal one(1.0, 3);
00579
00580
00581 Array<ADReal> tmp(result.length());
00582
00583 switch(order_)
00584 {
00585 case 0:
00586 tmp.resize(1);
00587 result.resize(1);
00588 tmp[0] = one;
00589 break;
00590 case 1:
00591 result.resize(4);
00592 tmp.resize(4);
00593 tmp[0] = 1.0 - x - y - z;
00594 tmp[1] = x;
00595 tmp[2] = y;
00596 tmp[3] = z;
00597 break;
00598 case 2:
00599 result.resize(10);
00600 tmp.resize(10);
00601 tmp[0] = (1.0-x-y-z)*(1.0-2.0*x-2.0*y-2.0*z);
00602 tmp[1] = 2.0*x*(x-0.5);
00603 tmp[2] = 2.0*y*(y-0.5);
00604 tmp[3] = 2.0*z*(z-0.5);
00605 tmp[9] = 4.0*x*(1.0-x-y-z);
00606 tmp[6] = 4.0*x*y;
00607 tmp[8] = 4.0*y*(1.0-x-y-z);
00608 tmp[7] = 4.0*z*(1.0-x-y-z);
00609 tmp[5] = 4.0*x*z;
00610 tmp[4] = 4.0*y*z;
00611 break;
00612 default:
00613 #ifndef TRILINOS_7
00614 SUNDANCE_ERROR("Lagrange::evalOnTet poly order > 2");
00615 #else
00616 SUNDANCE_ERROR7("Lagrange::evalOnTet poly order > 2");
00617 #endif
00618 }
00619
00620 for (int i=0; i<tmp.length(); i++)
00621 {
00622 if (deriv.order()==0) result[i] = tmp[i].value();
00623 else
00624 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00625 }
00626 }
00627
00628 void Lagrange::evalOnBrick(const Point& pt,
00629 const MultiIndex& deriv,
00630 Array<double>& result) const
00631 {
00632 ADReal x = ADReal(pt[0], 0, 3);
00633 ADReal y = ADReal(pt[1], 1, 3);
00634 ADReal z = ADReal(pt[2], 2, 3);
00635 ADReal one(1.0, 3);
00636
00637
00638 Array<ADReal> tmp(result.length());
00639
00640 switch(order_)
00641 {
00642 case 0:
00643 tmp.resize(1);
00644 result.resize(1);
00645 tmp[0] = one;
00646 break;
00647 case 1:
00648 result.resize(8);
00649 tmp.resize(8);
00650 tmp[0] = (1.0 - x)*(1.0 - y)*(1.0 - z);
00651 tmp[1] = (x)*(1.0 - y)*(1.0 - z);
00652 tmp[2] = (1.0 - x)*(y)*(1.0 - z);
00653 tmp[3] = (x)*(y)*(1.0 - z);
00654 tmp[4] = (1.0 - x)*(1.0 - y)*(z);
00655 tmp[5] = (x)*(1.0 - y)*(z);
00656 tmp[6] = (1.0 - x)*(y)*(z);
00657 tmp[7] = (x)*(y)*(z);
00658 break;
00659 case 2:
00660 result.resize(27);
00661 tmp.resize(27);
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672 tmp[0] = 27.0 - 27.0*x - 27.0*y - 27.0*z + 27.0*x*y + 27.0*x*z + 27.0*y*z - 27.0*x*y*z
00673 + 18.0*(x*x - 1)*(y-1)*(z-1)
00674 + 18.0*(y*y - 1)*(x-1)*(z-1)
00675 + 18.0*(z*z - 1)*(x-1)*(y-1)
00676 - 12.0*(x*x - 1)*(y*y - 1)*(z-1)
00677 - 12.0*(x*x - 1)*(z*z - 1)*(y-1)
00678 - 12.0*(y*y - 1)*(z*z - 1)*(x-1)
00679 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00680
00681 tmp[1] = 18.0 - 9.0*x - 18.0*y - 18.0*z + 9.0*x*y + 9.0*x*z + 18.0*y*z - 9.0*x*y*z
00682 + 18.0*(x*x - 1)*(y-1)*(z-1)
00683 + 9.0*(y*y - 1)*(x-1)*(z-1) - 3.0*(y*y - 1)*(x+1)*(z-1)
00684 + 9.0*(z*z - 1)*(x-1)*(y-1) - 3.0*(z*z - 1)*(x+1)*(y-1)
00685 - 12.0*(x*x - 1)*(y*y - 1)*(z-1)
00686 - 12.0*(x*x - 1)*(z*z - 1)*(y-1)
00687 - 6.0*(y*y - 1)*(z*z - 1)*(x-1) + 2.0*(y*y - 1)*(z*z - 1)*(x+1)
00688 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00689
00690 tmp[2] = 18.0 - 18.0*x - 9.0*y - 18.0*z + 9.0*x*y + 18.0*x*z + 9.0*y*z - 9.0*x*y*z
00691 + 9.0*(x*x - 1)*(y-1)*(z-1) - 3.0*(x*x - 1)*(y+1)*(z-1)
00692 + 18.0*(y*y - 1)*(x-1)*(z-1)
00693 + 9.0*(z*z - 1)*(x-1)*(y-1) - 3.0*(z*z - 1)*(x-1)*(y+1)
00694 - 12.0*(x*x - 1)*(y*y - 1)*(z-1)
00695 - 6.0*(x*x - 1)*(z*z - 1)*(y-1) + 2.0*(x*x - 1)*(z*z - 1)*(y+1)
00696 - 12.0*(y*y - 1)*(z*z - 1)*(x-1)
00697 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00698
00699 tmp[3] = 12.0 - 6.0*x - 6.0*y - 12.0*z + 3.0*x*y + 6.0*x*z + 6.0*y*z - 3.0*x*y*z
00700 + 9.0*(x*x - 1)*(y-1)*(z-1) - 3.0*(x*x - 1)*(y+1)*(z-1)
00701 + 9.0*(y*y - 1)*(x-1)*(z-1) - 3.0*(y*y - 1)*(x+1)*(z-1)
00702 + 4.5*(z*z - 1)*(x-1)*(y-1) - 1.5*(z*z - 1)*(x-1)*(y+1) - 1.5*(z*z - 1)*(x+1)*(y-1) + 0.5*(z*z - 1)*(x+1)*(y+1)
00703 - 12.0*(x*x - 1)*(y*y - 1)*(z-1)
00704 - 6.0*(x*x - 1)*(z*z - 1)*(y-1) + 2.0*(x*x - 1)*(z*z - 1)*(y+1)
00705 - 6.0*(y*y - 1)*(z*z - 1)*(x-1) + 2.0*(y*y - 1)*(z*z - 1)*(x+1)
00706 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00707
00708 tmp[4] = 18.0 - 18.0*x - 18.0*y - 9.0*z + 18.0*x*y + 9.0*x*z + 9.0*y*z - 9.0*x*y*z
00709 + 9.0*(x*x - 1)*(y-1)*(z-1) - 3.0*(x*x - 1)*(y-1)*(z+1)
00710 + 9.0*(y*y - 1)*(x-1)*(z-1) - 3.0*(y*y - 1)*(x-1)*(z+1)
00711 + 18.0*(z*z - 1)*(x-1)*(y-1)
00712 - 6.0*(x*x - 1)*(y*y - 1)*(z-1) + 2.0*(x*x - 1)*(y*y - 1)*(z+1)
00713 - 12.0*(x*x - 1)*(z*z - 1)*(y-1)
00714 - 12.0*(y*y - 1)*(z*z - 1)*(x-1)
00715 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00716
00717 tmp[5] = 12.0 - 6.0*x - 12.0*y - 6.0*z + 6.0*x*y + 3.0*x*z + 6.0*y*z - 3.0*x*y*z
00718 + 9.0*(x*x - 1)*(y-1)*(z-1) - 3.0*(x*x - 1)*(y-1)*(z+1)
00719 + 4.5*(y*y - 1)*(x-1)*(z-1) - 1.5*(y*y - 1)*(x-1)*(z+1) - 1.5*(y*y - 1)*(x+1)*(z-1) + 0.5*(y*y - 1)*(x+1)*(z+1)
00720 + 9.0*(z*z - 1)*(x-1)*(y-1) - 3.0*(z*z - 1)*(x+1)*(y-1)
00721 - 6.0*(x*x - 1)*(y*y - 1)*(z-1) + 2.0*(x*x - 1)*(y*y - 1)*(z+1)
00722 - 12.0*(x*x - 1)*(z*z - 1)*(y-1)
00723 - 6.0*(y*y - 1)*(z*z - 1)*(x-1) + 2.0*(y*y - 1)*(z*z - 1)*(x+1)
00724 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00725
00726 tmp[6] = 12.0 - 12.0*x - 6.0*y - 6.0*z + 6.0*x*y + 6.0*x*z + 3.0*y*z - 3.0*x*y*z
00727 + 4.5*(x*x - 1)*(y-1)*(z-1) - 1.5*(x*x - 1)*(y-1)*(z+1) - 1.5*(x*x - 1)*(y+1)*(z-1) + 0.5*(x*x - 1)*(y+1)*(z+1)
00728 + 9.0*(y*y - 1)*(x-1)*(z-1) - 3.0*(y*y - 1)*(x-1)*(z+1)
00729 + 9.0*(z*z - 1)*(x-1)*(y-1) - 3.0*(z*z - 1)*(x-1)*(y+1)
00730 - 6.0*(x*x - 1)*(y*y - 1)*(z-1) + 2.0*(x*x - 1)*(y*y - 1)*(z+1)
00731 - 6.0*(x*x - 1)*(z*z - 1)*(y-1) + 2.0*(x*x - 1)*(z*z - 1)*(y+1)
00732 - 12.0*(y*y - 1)*(z*z - 1)*(x-1)
00733 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00734
00735 tmp[7] = 8.0 - 4.0*x - 4.0*y - 4.0*z + 2.0*x*y + 2.0*x*z + 2.0*y*z - 1.0*x*y*z
00736 + 4.5*(x*x - 1)*(y-1)*(z-1) - 1.5*(x*x - 1)*(y-1)*(z+1) - 1.5*(x*x - 1)*(y+1)*(z-1) + 0.5*(x*x - 1)*(y+1)*(z+1)
00737 + 4.5*(y*y - 1)*(x-1)*(z-1) - 1.5*(y*y - 1)*(x-1)*(z+1) - 1.5*(y*y - 1)*(x+1)*(z-1) + 0.5*(y*y - 1)*(x+1)*(z+1)
00738 + 4.5*(z*z - 1)*(x-1)*(y-1) - 1.5*(z*z - 1)*(x-1)*(y+1) - 1.5*(z*z - 1)*(x+1)*(y-1) + 0.5*(z*z - 1)*(x+1)*(y+1)
00739 - 6.0*(x*x - 1)*(y*y - 1)*(z-1) + 2.0*(x*x - 1)*(y*y - 1)*(z+1)
00740 - 6.0*(x*x - 1)*(z*z - 1)*(y-1) + 2.0*(x*x - 1)*(z*z - 1)*(y+1)
00741 - 6.0*(y*y - 1)*(z*z - 1)*(x-1) + 2.0*(y*y - 1)*(z*z - 1)*(x+1)
00742 + 8.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00743
00744 tmp[8] = -36.0 + 36.0*x + 36.0*y + 36.0*z - 36.0*x*y - 36.0*x*z - 36.0*y*z + 36.0*x*y*z
00745 - 36.0*(x*x - 1)*(y-1)*(z-1)
00746 - 24.0*(y*y - 1)*(x-1)*(z-1)
00747 - 24.0*(z*z - 1)*(x-1)*(y-1)
00748 + 24.0*(x*x - 1)*(y*y - 1)*(z-1)
00749 + 24.0*(x*x - 1)*(z*z - 1)*(y-1)
00750 + 16.0*(y*y - 1)*(z*z - 1)*(x-1)
00751 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00752
00753 tmp[9] = -36.0 + 36.0*x + 36.0*y + 36.0*z - 36.0*x*y - 36.0*x*z - 36.0*y*z + 36.0*x*y*z
00754 - 24.0*(x*x - 1)*(y-1)*(z-1)
00755 - 36.0*(y*y - 1)*(x-1)*(z-1)
00756 - 24.0*(z*z - 1)*(x-1)*(y-1)
00757 + 24.0*(x*x - 1)*(y*y - 1)*(z-1)
00758 + 16.0*(x*x - 1)*(z*z - 1)*(y-1)
00759 + 24.0*(y*y - 1)*(z*z - 1)*(x-1)
00760 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00761
00762 tmp[10] = -36.0 + 36.0*x + 36.0*y + 36.0*z - 36.0*x*y - 36.0*x*z - 36.0*y*z + 36.0*x*y*z
00763 - 24.0*(x*x - 1)*(y-1)*(z-1)
00764 - 24.0*(y*y - 1)*(x-1)*(z-1)
00765 - 36.0*(z*z - 1)*(x-1)*(y-1)
00766 + 16.0*(x*x - 1)*(y*y - 1)*(z-1)
00767 + 24.0*(x*x - 1)*(z*z - 1)*(y-1)
00768 + 24.0*(y*y - 1)*(z*z - 1)*(x-1)
00769 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00770
00771 tmp[11] = -24.0 + 12.0*x + 24.0*y + 24.0*z - 12.0*x*y - 12.0*x*z - 24.0*y*z + 12.0*x*y*z
00772 - 24.0*(x*x - 1)*(y-1)*(z-1)
00773 - 18.0*(y*y - 1)*(x-1)*(z-1) + 6.0*(y*y - 1)*(x+1)*(z-1)
00774 - 12.0*(z*z - 1)*(x-1)*(y-1) + 4.0*(z*z - 1)*(x+1)*(y-1)
00775 + 24.0*(x*x - 1)*(y*y - 1)*(z-1)
00776 + 16.0*(x*x - 1)*(z*z - 1)*(y-1)
00777 + 12.0*(y*y - 1)*(z*z - 1)*(x-1) - 4.0*(y*y - 1)*(z*z - 1)*(x+1)
00778 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00779
00780 tmp[12] = -24.0 + 12.0*x + 24.0*y + 24.0*z - 12.0*x*y - 12.0*x*z - 24.0*y*z + 12.0*x*y*z
00781 - 24.0*(x*x - 1)*(y-1)*(z-1)
00782 - 12.0*(y*y - 1)*(x-1)*(z-1) + 4.0*(y*y - 1)*(x+1)*(z-1)
00783 - 18.0*(z*z - 1)*(x-1)*(y-1) + 6.0*(z*z - 1)*(x+1)*(y-1)
00784 + 16.0*(x*x - 1)*(y*y - 1)*(z-1)
00785 + 24.0*(x*x - 1)*(z*z - 1)*(y-1)
00786 + 12.0*(y*y - 1)*(z*z - 1)*(x-1) - 4.0*(y*y - 1)*(z*z - 1)*(x+1)
00787 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00788
00789 tmp[13] = -24.0 + 24.0*x + 12.0*y + 24.0*z - 12.0*x*y - 24.0*x*z - 12.0*y*z + 12.0*x*y*z
00790 - 18.0*(x*x - 1)*(y-1)*(z-1) + 6.0*(x*x - 1)*(y+1)*(z-1)
00791 - 24.0*(y*y - 1)*(x-1)*(z-1)
00792 - 12.0*(z*z - 1)*(x-1)*(y-1) + 4.0*(z*z - 1)*(x-1)*(y+1)
00793 + 24.0*(x*x - 1)*(y*y - 1)*(z-1)
00794 + 12.0*(x*x - 1)*(z*z - 1)*(y-1) - 4.0*(x*x - 1)*(z*z - 1)*(y+1)
00795 + 16.0*(y*y - 1)*(z*z - 1)*(x-1)
00796 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00797
00798 tmp[14] = -24.0 + 24.0*x + 12.0*y + 24.0*z - 12.0*x*y - 24.0*x*z - 12.0*y*z + 12.0*x*y*z
00799 - 12.0*(x*x - 1)*(y-1)*(z-1) + 4.0*(x*x - 1)*(y+1)*(z-1)
00800 - 24.0*(y*y - 1)*(x-1)*(z-1)
00801 - 18.0*(z*z - 1)*(x-1)*(y-1) + 6.0*(z*z - 1)*(x-1)*(y+1)
00802 + 16.0*(x*x - 1)*(y*y - 1)*(z-1)
00803 + 12.0*(x*x - 1)*(z*z - 1)*(y-1) - 4.0*(x*x - 1)*(z*z - 1)*(y+1)
00804 + 24.0*(y*y - 1)*(z*z - 1)*(x-1)
00805 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00806
00807 tmp[15] = -16.0 + 8.0*x + 8.0*y + 16.0*z - 4.0*x*y - 8.0*x*z - 8.0*y*z + 4.0*x*y*z
00808 - 12.0*(x*x - 1)*(y-1)*(z-1) + 4.0*(x*x - 1)*(y+1)*(z-1)
00809 - 12.0*(y*y - 1)*(x-1)*(z-1) + 4.0*(y*y - 1)*(x+1)*(z-1)
00810 - 9.0*(z*z - 1)*(x-1)*(y-1) + 3.0*(z*z - 1)*(x-1)*(y+1) + 3.0*(z*z - 1)*(x+1)*(y-1) - 1.0*(z*z - 1)*(x+1)*(y+1)
00811 + 16.0*(x*x - 1)*(y*y - 1)*(z-1)
00812 + 12.0*(x*x - 1)*(z*z - 1)*(y-1) - 4.0*(x*x - 1)*(z*z - 1)*(y+1)
00813 + 12.0*(y*y - 1)*(z*z - 1)*(x-1) - 4.0*(y*y - 1)*(z*z - 1)*(x+1)
00814 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00815
00816 tmp[16] = -24.0 + 24.0*x + 24.0*y + 12.0*z - 24.0*x*y - 12.0*x*z - 12.0*y*z + 12.0*x*y*z
00817 - 18.0*(x*x - 1)*(y-1)*(z-1) + 6.0*(x*x - 1)*(y-1)*(z+1)
00818 - 12.0*(y*y - 1)*(x-1)*(z-1) + 4.0*(y*y - 1)*(x-1)*(z+1)
00819 - 24.0*(z*z - 1)*(x-1)*(y-1)
00820 + 12.0*(x*x - 1)*(y*y - 1)*(z-1) - 4.0*(x*x - 1)*(y*y - 1)*(z+1)
00821 + 24.0*(x*x - 1)*(z*z - 1)*(y-1)
00822 + 16.0*(y*y - 1)*(z*z - 1)*(x-1)
00823 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00824
00825 tmp[17] = -24.0 + 24.0*x + 24.0*y + 12.0*z - 24.0*x*y - 12.0*x*z - 12.0*y*z + 12.0*x*y*z
00826 - 12.0*(x*x - 1)*(y-1)*(z-1) + 4.0*(x*x - 1)*(y-1)*(z+1)
00827 - 18.0*(y*y - 1)*(x-1)*(z-1) + 6.0*(y*y - 1)*(x-1)*(z+1)
00828 - 24.0*(z*z - 1)*(x-1)*(y-1)
00829 + 12.0*(x*x - 1)*(y*y - 1)*(z-1) - 4.0*(x*x - 1)*(y*y - 1)*(z+1)
00830 + 16.0*(x*x - 1)*(z*z - 1)*(y-1)
00831 + 24.0*(y*y - 1)*(z*z - 1)*(x-1)
00832 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00833
00834 tmp[18] = -16.0 + 8.0*x + 16.0*y + 8.0*z - 8.0*x*y - 4.0*x*z - 8.0*y*z + 4.0*x*y*z
00835 - 12.0*(x*x - 1)*(y-1)*(z-1) + 4.0*(x*x - 1)*(y-1)*(z+1)
00836 - 9.0*(y*y - 1)*(x-1)*(z-1) + 3.0*(y*y - 1)*(x-1)*(z+1) + 3.0*(y*y - 1)*(x+1)*(z-1) - 1.0*(y*y - 1)*(x+1)*(z+1)
00837 - 12.0*(z*z - 1)*(x-1)*(y-1) + 4.0*(z*z - 1)*(x+1)*(y-1)
00838 + 12.0*(x*x - 1)*(y*y - 1)*(z-1) - 4.0*(x*x - 1)*(y*y - 1)*(z+1)
00839 + 16.0*(x*x - 1)*(z*z - 1)*(y-1)
00840 + 12.0*(y*y - 1)*(z*z - 1)*(x-1) - 4.0*(y*y - 1)*(z*z - 1)*(x+1)
00841 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00842
00843 tmp[19] = -16.0 + 16.0*x + 8.0*y + 8.0*z - 8.0*x*y - 8.0*x*z - 4.0*y*z + 4.0*x*y*z
00844 - 9.0*(x*x - 1)*(y-1)*(z-1) + 3.0*(x*x - 1)*(y-1)*(z+1) + 3.0*(x*x - 1)*(y+1)*(z-1) - 1.0*(x*x - 1)*(y+1)*(z+1)
00845 - 12.0*(y*y - 1)*(x-1)*(z-1) + 4.0*(y*y - 1)*(x-1)*(z+1)
00846 - 12.0*(z*z - 1)*(x-1)*(y-1) + 4.0*(z*z - 1)*(x-1)*(y+1)
00847 + 12.0*(x*x - 1)*(y*y - 1)*(z-1) - 4.0*(x*x - 1)*(y*y - 1)*(z+1)
00848 + 12.0*(x*x - 1)*(z*z - 1)*(y-1) - 4.0*(x*x - 1)*(z*z - 1)*(y+1)
00849 + 16.0*(y*y - 1)*(z*z - 1)*(x-1)
00850 - 16.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00851
00852 tmp[20] = 48.0 - 48.0*x - 48.0*y - 48.0*z + 48.0*x*y + 48.0*x*z + 48.0*y*z - 48.0*x*y*z
00853 + 48.0*(x*x - 1)*(y-1)*(z-1)
00854 + 48.0*(y*y - 1)*(x-1)*(z-1)
00855 + 32.0*(z*z - 1)*(x-1)*(y-1)
00856 - 48.0*(x*x - 1)*(y*y - 1)*(z-1)
00857 - 32.0*(x*x - 1)*(z*z - 1)*(y-1)
00858 - 32.0*(y*y - 1)*(z*z - 1)*(x-1)
00859 + 32.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00860
00861 tmp[21] = 48.0 - 48.0*x - 48.0*y - 48.0*z + 48.0*x*y + 48.0*x*z + 48.0*y*z - 48.0*x*y*z
00862 + 48.0*(x*x - 1)*(y-1)*(z-1)
00863 + 32.0*(y*y - 1)*(x-1)*(z-1)
00864 + 48.0*(z*z - 1)*(x-1)*(y-1)
00865 - 32.0*(x*x - 1)*(y*y - 1)*(z-1)
00866 - 48.0*(x*x - 1)*(z*z - 1)*(y-1)
00867 - 32.0*(y*y - 1)*(z*z - 1)*(x-1)
00868 + 32.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00869
00870 tmp[22] = 48.0 - 48.0*x - 48.0*y - 48.0*z + 48.0*x*y + 48.0*x*z + 48.0*y*z - 48.0*x*y*z
00871 + 32.0*(x*x - 1)*(y-1)*(z-1)
00872 + 48.0*(y*y - 1)*(x-1)*(z-1)
00873 + 48.0*(z*z - 1)*(x-1)*(y-1)
00874 - 32.0*(x*x - 1)*(y*y - 1)*(z-1)
00875 - 32.0*(x*x - 1)*(z*z - 1)*(y-1)
00876 - 48.0*(y*y - 1)*(z*z - 1)*(x-1)
00877 + 32.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00878
00879 tmp[23] = 32.0 - 16.0*x - 32.0*y - 32.0*z + 16.0*x*y + 16.0*x*z + 32.0*y*z - 16.0*x*y*z
00880 + 32.0*(x*x - 1)*(y-1)*(z-1)
00881 + 24.0*(y*y - 1)*(x-1)*(z-1) - 8.0*(y*y - 1)*(x+1)*(z-1)
00882 + 24.0*(z*z - 1)*(x-1)*(y-1) - 8.0*(z*z - 1)*(x+1)*(y-1)
00883 - 32.0*(x*x - 1)*(y*y - 1)*(z-1)
00884 - 32.0*(x*x - 1)*(z*z - 1)*(y-1)
00885 - 24.0*(y*y - 1)*(z*z - 1)*(x-1) + 8.0*(y*y - 1)*(z*z - 1)*(x+1)
00886 + 32.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00887
00888 tmp[24] = 32.0 - 32.0*x - 16.0*y - 32.0*z + 16.0*x*y + 32.0*x*z + 16.0*y*z - 16.0*x*y*z
00889 + 24.0*(x*x - 1)*(y-1)*(z-1)- 8.0*(x*x - 1)*(y+1)*(z-1)
00890 + 32.0*(y*y - 1)*(x-1)*(z-1)
00891 + 24.0*(z*z - 1)*(x-1)*(y-1) - 8.0*(z*z - 1)*(x-1)*(y+1)
00892 - 32.0*(x*x - 1)*(y*y - 1)*(z-1)
00893 - 24.0*(x*x - 1)*(z*z - 1)*(y-1) + 8.0*(x*x - 1)*(z*z - 1)*(y+1)
00894 - 32.0*(y*y - 1)*(z*z - 1)*(x-1)
00895 + 32.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00896
00897 tmp[25] = 32.0 - 32.0*x - 32.0*y - 16.0*z + 32.0*x*y + 16.0*x*z + 16.0*y*z - 16.0*x*y*z
00898 + 24.0*(x*x - 1)*(y-1)*(z-1) - 8.0*(x*x - 1)*(y-1)*(z+1)
00899 + 24.0*(y*y - 1)*(x-1)*(z-1) - 8.0*(y*y - 1)*(x-1)*(z+1)
00900 + 32.0*(z*z - 1)*(x-1)*(y-1)
00901 - 24.0*(x*x - 1)*(y*y - 1)*(z-1) + 8.0*(x*x - 1)*(y*y - 1)*(z+1)
00902 - 32.0*(x*x - 1)*(z*z - 1)*(y-1)
00903 - 32.0*(y*y - 1)*(z*z - 1)*(x-1)
00904 + 32.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00905
00906 tmp[26] = -64.0 + 64.0*x + 64.0*y + 64.0*z - 64.0*x*y - 64.0*x*z - 64.0*y*z + 64.0*x*y*z
00907 - 64.0*(x*x - 1)*(y-1)*(z-1)
00908 - 64.0*(y*y - 1)*(x-1)*(z-1)
00909 - 64.0*(z*z - 1)*(x-1)*(y-1)
00910 + 64.0*(x*x - 1)*(y*y - 1)*(z-1)
00911 + 64.0*(x*x - 1)*(z*z - 1)*(y-1)
00912 + 64.0*(y*y - 1)*(z*z - 1)*(x-1)
00913 - 64.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
00914
00915 break;
00916 default:
00917 #ifndef TRILINOS_7
00918 SUNDANCE_ERROR("Lagrange::evalOnBrick poly order > 2");
00919 #else
00920 SUNDANCE_ERROR7("Lagrange::evalOnBrick poly order > 2");
00921 #endif
00922 }
00923
00924 for (int i=0; i<tmp.length(); i++)
00925 {
00926 if (deriv.order()==0) result[i] = tmp[i].value();
00927 else
00928 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00929 }
00930 }
00931
00932 void Lagrange::getConstrainsForHNDoF(
00933 const int indexInParent,
00934 const int maxCellDim,
00935 const int maxNrChild,
00936 const int facetDim,
00937 const int facetIndex,
00938 const int nodeIndex,
00939 Array<int>& localDoFs,
00940 Array<int>& parentFacetDim,
00941 Array<int>& parentFacetIndex,
00942 Array<int>& parentFacetNode,
00943 Array<double>& coefs
00944 ) {
00945
00946
00947
00948
00949
00950
00951 Array<Array<Array<double> > > result;
00952 const SpatialDerivSpecifier deriv;
00953 int index = 0;
00954
00955
00956
00957 localDoFs.resize(0);
00958 coefs.resize(0);
00959 parentFacetDim.resize(0);
00960 parentFacetIndex.resize(0);
00961 parentFacetNode.resize(0);
00962
00963 CellType maximalCellType;
00964
00965 switch (maxCellDim){
00966 case 2:{
00967 maximalCellType = QuadCell;
00968
00969 switch (maxNrChild){
00970 case 9:{
00971 Point dofPos( 0.0 , 0.0 );
00972
00973 SUNDANCE_OUT(this->verb() > 3 , ",maxCellDim=" << maxCellDim << ",facetDim="
00974 << facetDim << ",facetIndex=" << facetIndex << ",nodeIndex=" << nodeIndex
00975 << ",dofPos=" << dofPos);
00976 returnDoFPositionOnRef( maxCellDim,
00977 facetDim, facetIndex, nodeIndex, dofPos);
00978 dofPos = dofPos/3.0;
00979 SUNDANCE_OUT(this->verb() > 3 , "dofPos=" << dofPos );
00980
00981
00982 switch (indexInParent){
00983 case 0: { break;}
00984 case 1: {dofPos[0] += 1.0/3.0; break;}
00985 case 2: {dofPos[0] += 2.0/3.0; break;}
00986 case 3: {dofPos[0] += 2.0/3.0; dofPos[1] += 1.0/3.0; break;}
00987 case 5: {dofPos[1] += 1.0/3.0; break;}
00988 case 6: {dofPos[1] += 2.0/3.0; break;}
00989 case 7: {dofPos[0] += 1.0/3.0; dofPos[1] += 2.0/3.0; break;}
00990 case 8: {dofPos[0] += 2.0/3.0; dofPos[1] += 2.0/3.0; break;}
00991 }
00992
00993 SUNDANCE_OUT(this->verb() > 3 , "dofPos=" << dofPos );
00994 Array<Point> tmp_points(1);
00995 tmp_points[0] = dofPos;
00996 refEval( maximalCellType , tmp_points ,
00997 deriv, result , 4);
00998
00999 break;}
01000 case 4:{
01001 break;
01002 }
01003 }
01004 break;}
01005 case 3:{
01006
01007 maximalCellType = BrickCell;
01008
01009 Point dofPos( 0.0 , 0.0 , 0.0);
01010
01011 SUNDANCE_OUT(this->verb() > 3 , ",maxCellDim=" << maxCellDim << ",facetDim="
01012 << facetDim << ",facetIndex=" << facetIndex << ",nodeIndex=" << nodeIndex
01013 << ",dofPos=" << dofPos);
01014 returnDoFPositionOnRef( maxCellDim,
01015 facetDim, facetIndex, nodeIndex, dofPos);
01016 dofPos = dofPos/3.0;
01017 SUNDANCE_OUT(this->verb() > 3 , "dofPos=" << dofPos );
01018
01019
01020 switch (maxNrChild){
01021 case 27:{
01022 switch (indexInParent){
01023 case 0: {dofPos[0] += 0.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 0.0/3.0; break;}
01024 case 1: {dofPos[0] += 1.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 0.0/3.0; break;}
01025 case 2: {dofPos[0] += 2.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 0.0/3.0; break;}
01026 case 3: {dofPos[0] += 0.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 0.0/3.0; break;}
01027 case 4: {dofPos[0] += 1.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 0.0/3.0; break;}
01028 case 5: {dofPos[0] += 2.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 0.0/3.0; break;}
01029 case 6: {dofPos[0] += 0.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 0.0/3.0; break;}
01030 case 7: {dofPos[0] += 1.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 0.0/3.0; break;}
01031 case 8: {dofPos[0] += 2.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 0.0/3.0; break;}
01032 case 9: {dofPos[0] += 0.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 1.0/3.0; break;}
01033 case 10: {dofPos[0] += 1.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 1.0/3.0; break;}
01034 case 11: {dofPos[0] += 2.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 1.0/3.0; break;}
01035 case 12: {dofPos[0] += 0.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 1.0/3.0; break;}
01036
01037 case 13: {dofPos[0] += 1.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 1.0/3.0; break;}
01038 case 14: {dofPos[0] += 2.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 1.0/3.0; break;}
01039 case 15: {dofPos[0] += 0.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 1.0/3.0; break;}
01040 case 16: {dofPos[0] += 1.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 1.0/3.0; break;}
01041 case 17: {dofPos[0] += 2.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 1.0/3.0; break;}
01042 case 18: {dofPos[0] += 0.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 2.0/3.0; break;}
01043 case 19: {dofPos[0] += 1.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 2.0/3.0; break;}
01044 case 20: {dofPos[0] += 2.0/3.0; dofPos[1] += 0.0/3.0; dofPos[2] += 2.0/3.0; break;}
01045 case 21: {dofPos[0] += 0.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 2.0/3.0; break;}
01046 case 22: {dofPos[0] += 1.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 2.0/3.0; break;}
01047 case 23: {dofPos[0] += 2.0/3.0; dofPos[1] += 1.0/3.0; dofPos[2] += 2.0/3.0; break;}
01048 case 24: {dofPos[0] += 0.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 2.0/3.0; break;}
01049 case 25: {dofPos[0] += 1.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 2.0/3.0; break;}
01050 case 26: {dofPos[0] += 2.0/3.0; dofPos[1] += 2.0/3.0; dofPos[2] += 2.0/3.0; break;}
01051 }
01052 break;}
01053 case 8:{
01054 break;}
01055 }
01056
01057 SUNDANCE_OUT(this->verb() > 3 , "dofPos=" << dofPos );
01058 Array<Point> tmp_points(1);
01059 tmp_points[0] = dofPos;
01060 refEval( maximalCellType , tmp_points ,
01061 deriv, result , 4);
01062
01063 break;}
01064 }
01065
01066
01067 Array<double>& res = result[0][0];
01068
01069
01070
01071 if (doFInfromationCalculated_ == false)
01072 {
01073 getDoFsOrdered( maximalCellType, res.size() , facetD_, facetI_, facetN_);
01074 doFInfromationCalculated_ = true;
01075 }
01076
01077
01078 index = 0;
01079 for ( int i=0 ; i < res.size() ; i++){
01080
01081 if ( fabs(res[i]) > 1e-5 ){
01082 SUNDANCE_OUT(this->verb() > 3 , "found DoF i=" << i );
01083 localDoFs.resize(index+1);
01084 coefs.resize(index+1);
01085 parentFacetDim.resize(index+1);
01086 parentFacetIndex.resize(index+1);
01087 parentFacetNode.resize(index+1);
01088 localDoFs[index] = i;
01089 coefs[index] = res[i];
01090 parentFacetDim[index] = facetD_[i];
01091 parentFacetIndex[index] = facetI_[i];
01092 parentFacetNode[index] = facetN_[i];
01093 index++;
01094 }
01095 }
01096
01097 SUNDANCE_OUT(this->verb() > 3 , "localDoFs=" << localDoFs );
01098 SUNDANCE_OUT(this->verb() > 3 , "coefs=" << coefs );
01099 SUNDANCE_OUT(this->verb() > 3 , "parentFacetDim=" << parentFacetDim );
01100 SUNDANCE_OUT(this->verb() > 3 , "parentFacetIndex=" << parentFacetIndex );
01101 SUNDANCE_OUT(this->verb() > 3 , "parentFacetNode=" << parentFacetNode );
01102
01103 }
01104
01105
01106 void Lagrange::returnDoFPositionOnRef(
01107 const int maxCellDim,
01108 const int facetDim,
01109 const int facetIndex,
01110 const int nodeIndex,
01111 Point& pos) const {
01112
01113 switch (facetDim){
01114 case 0:{
01115 switch (maxCellDim){
01116 case 2:{
01117 switch (facetIndex){
01118 case 0: {pos[0] = 0.0; pos[1] = 0.0; break;}
01119 case 1: {pos[0] = 1.0; pos[1] = 0.0; break;}
01120 case 2: {pos[0] = 0.0; pos[1] = 1.0; break;}
01121 case 3: {pos[0] = 1.0; pos[1] = 1.0; break;}
01122 }
01123 break;}
01124 case 3:{
01125 switch (facetIndex){
01126 case 0: {pos[0] = 0.0; pos[1] = 0.0; pos[2] = 0.0; break;}
01127 case 1: {pos[0] = 1.0; pos[1] = 0.0; pos[2] = 0.0; break;}
01128 case 2: {pos[0] = 0.0; pos[1] = 1.0; pos[2] = 0.0; break;}
01129 case 3: {pos[0] = 1.0; pos[1] = 1.0; pos[2] = 0.0; break;}
01130 case 4: {pos[0] = 0.0; pos[1] = 0.0; pos[2] = 1.0; break;}
01131 case 5: {pos[0] = 1.0; pos[1] = 0.0; pos[2] = 1.0; break;}
01132 case 6: {pos[0] = 0.0; pos[1] = 1.0; pos[2] = 1.0; break;}
01133 case 7: {pos[0] = 1.0; pos[1] = 1.0; pos[2] = 1.0; break;}
01134 }
01135 break;}
01136 }
01137 break;}
01138 case 1:{
01139 switch (maxCellDim){
01140 case 2:{
01141
01142 double posOnEdge = (1.0/ (double)(order())) +(double)(nodeIndex) / (double)(order()-1);
01143 switch (facetIndex){
01144 case 0: {pos[0] = 0.0+posOnEdge; pos[1] = 0.0; break;}
01145 case 1: {pos[0] = 0.0; pos[1] = 0.0+posOnEdge; break;}
01146 case 2: {pos[0] = 1.0; pos[1] = 0.0+posOnEdge; break;}
01147 case 3: {pos[0] = 0.0+posOnEdge; pos[1] = 1.0; break;}
01148 }
01149 break;}
01150 case 3:{
01151 double posOnEdge = (1.0/ (double)(order())) +(double)(nodeIndex) / (double)(order()-1);
01152 switch (facetIndex){
01153 case 0: {pos[0] = 0.0+posOnEdge; pos[1] = 0.0; pos[2] = 0.0; break;}
01154 case 1: {pos[0] = 0.0; pos[1] = 0.0+posOnEdge; pos[2] = 0.0; break;}
01155 case 2: {pos[0] = 0.0; pos[1] = 0.0; pos[2] = 0.0+posOnEdge; break;}
01156 case 3: {pos[0] = 1.0; pos[1] = 0.0+posOnEdge; pos[2] = 0.0; break;}
01157 case 4: {pos[0] = 1.0; pos[1] = 0.0; pos[2] = 0.0+posOnEdge; break;}
01158 case 5: {pos[0] = 0.0+posOnEdge; pos[1] = 1.0; pos[2] = 0.0; break;}
01159 case 6: {pos[0] = 0.0; pos[1] = 1.0; pos[2] = 0.0+posOnEdge; break;}
01160 case 7: {pos[0] = 1.0; pos[1] = 1.0; pos[2] = 0.0+posOnEdge; break;}
01161 case 8: {pos[0] = 0.0+posOnEdge; pos[1] = 0.0; pos[2] = 1.0; break;}
01162 case 9: {pos[0] = 0.0; pos[1] = 0.0+posOnEdge; pos[2] = 1.0; break;}
01163 case 10: {pos[0] = 1.0; pos[1] = 0.0+posOnEdge; pos[2] = 1.0; break;}
01164 case 11: {pos[0] = 0.0+posOnEdge; pos[1] = 1.0; pos[2] = 1.0; break;}
01165 }
01166 break;}
01167 }
01168 break;}
01169 case 2:{
01170
01171 int nrDoFoFace = (order()-1)*(order()-1);
01172 double posOnFace1 = (1.0/ (double)(order())) +(double)(nodeIndex % nrDoFoFace) / (double)(order()-1);
01173 double posOnFace2 = (1.0/ (double)(order())) +(double)(nodeIndex / nrDoFoFace) / (double)(order()-1);
01174 switch (facetIndex){
01175 case 0: {pos[0] = 0.0+posOnFace1; pos[1] = 0.0+posOnFace2; pos[2] = 0.0; break;}
01176 case 1: {pos[0] = 0.0+posOnFace1; pos[1] = 0.0; pos[2] = 0.0+posOnFace2; break;}
01177 case 2: {pos[0] = 0.0; pos[1] = 0.0+posOnFace1; pos[2] = 0.0+posOnFace2; break;}
01178 case 3: {pos[0] = 1.0; pos[1] = 0.0+posOnFace1; pos[2] = 0.0+posOnFace2; break;}
01179 case 4: {pos[0] = 0.0+posOnFace1; pos[1] = 1.0; pos[2] = 0.0+posOnFace2; break;}
01180 case 5: {pos[0] = 0.0+posOnFace1; pos[1] = 0.0+posOnFace2; pos[2] = 1.0; break;}
01181 }
01182 break;}
01183 }
01184
01185 }
01186
01187 void Lagrange::getDoFsOrdered(
01188 const CellType maxCellDim,
01189 int nrDoF,
01190 Array<int>& facetD,
01191 Array<int>& facetI,
01192 Array<int>& facetN)
01193 {
01194
01195
01196 Array<Array<Array<int> > > dofs_struct;
01197 getReferenceDOFs(maxCellDim , maxCellDim , dofs_struct);
01198 bool getNext;
01199 facetD.resize(nrDoF); facetI.resize(nrDoF); facetN.resize(nrDoF);
01200 int dimo=0 , indexo=0 , nodo=0;
01201 for ( int i=0 ; i < nrDoF ; i++){
01202 facetD[i] = dimo;
01203 facetI[i] = indexo;
01204 facetN[i] = nodo;
01205
01206
01207 getNext = true;
01208 while ( (getNext) && (dimo < dofs_struct.size() )){
01209
01210
01211 if (dofs_struct[dimo].size() > (indexo+1))
01212 {
01213 if (dofs_struct[dimo][indexo].size() > (nodo+1)){
01214 nodo++; getNext = false;
01215 }else{
01216 nodo = 0;
01217 if (dofs_struct[dimo].size() > (indexo+1)){
01218 indexo++; getNext = false;
01219 }else{
01220 indexo = 0;
01221 if (dofs_struct.size() > (dimo+1)){
01222 dimo++; getNext = false;
01223 }else{
01224
01225 getNext = false;
01226 }
01227 }
01228 }
01229 }
01230 else
01231 {
01232 dimo++; indexo = 0; nodo = 0; getNext = false;
01233 }
01234
01235
01236 if (dofs_struct.size() <= dimo ) {dimo++; getNext = true; continue;}
01237 if (dofs_struct[dimo].size() <= indexo ) {dimo++; getNext = true; continue;}
01238 if (dofs_struct[dimo][indexo].size() <= nodo ) {dimo++; getNext = true; continue;}
01239 }
01240 }
01241 SUNDANCE_OUT(this->verb() > 3, "facetD=" << facetD );
01242 SUNDANCE_OUT(this->verb() > 3, "facetI=" << facetI );
01243 SUNDANCE_OUT(this->verb() > 3, "facetN=" << facetN );
01244
01245 }