Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
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; // -Wall
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:{ // TODO: ask Kevin about Q0 (P0)
00220          dofs.resize(3);
00221        // dofs[0] are DOFs at Points
00222        dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2), tuple(3));
00223        if (order() == 2)
00224        {
00225          // dofs[1] are the DOFs on the line
00226          dofs[1] = tuple<Aint>(tuple(4), tuple(5), tuple(6), tuple(7));
00227          //dofs[2] are DOFs inside the Cell
00228 
00229          //dofs[2] = tuple(Array<int>());
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        // dofs[0] are DOFs at Points
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          // dofs[1] are the DOFs on the line
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          //dofs[2] are DOFs on the faces
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 /* ---------- evaluation on different cell types -------------- */
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   switch(order_)
00396     {
00397     case 0:
00398       tmp[0] = one;
00399       break;
00400     case 1:
00401       tmp[0] = 1.0 - x;
00402       tmp[1] = x;
00403       break;
00404     case 2:
00405       tmp[0] = 2.0*(1.0-x)*(0.5-x);
00406       tmp[1] = 2.0*x*(x-0.5);
00407       tmp[2] = 4.0*x*(1.0-x);
00408       break;
00409     case 3:
00410       tmp[0] = 9.0/2.0 * (1.0 - x) * (1.0/3.0 - x) * (2.0/3.0 - x);
00411       tmp[1] = 9.0/2.0 * x * (x - 1.0/3.0) * (x - 2.0/3.0);
00412       tmp[2] = 27.0/2.0 * x * (1.0 - x) * (2.0/3.0 - x);
00413       tmp[3] = 27.0/2.0 * x * (1.0 - x) * (x - 1.0/3.0);
00414 
00415       break;
00416     default:
00417       SUNDANCE_ERROR("Lagrange::evalOnLine polynomial order > 2 has not been"
00418                      " implemented");
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; //These are the 4 basis functions for the first order
00533       break;
00534     case 2: //The basis functions calculated with Octave / Matlab
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: // tri linear basis functions for bricks
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: // second order function in 3D
00660           result.resize(27);
00661           tmp.resize(27);
00662                    // 0.0 + 0.0*x + 0.0*y + 0.0*z + 0.0*x*y + 0.0*x*z + 0.0*y*z + 0.0*x*y*z
00663                    //+ 0.0*(x*x - 1)*(y-1)*(z-1) + 0.0*(x*x - 1)*(y-1)*(z+1) + 0.0*(x*x - 1)*(y+1)*(z-1) + 0.0*(x*x - 1)*(y+1)*(z+1)
00664                    //+ 0.0*(y*y - 1)*(x-1)*(z-1) + 0.0*(y*y - 1)*(x-1)*(z+1) + 0.0*(y*y - 1)*(x+1)*(z-1) + 0.0*(y*y - 1)*(x+1)*(z+1)
00665                    //+ 0.0*(z*z - 1)*(x-1)*(y-1) + 0.0*(z*z - 1)*(x-1)*(y+1) + 0.0*(z*z - 1)*(x+1)*(y-1) + 0.0*(z*z - 1)*(x+1)*(y+1)
00666                    //+ 0.0*(x*x - 1)*(y*y - 1)*(z-1) + 0.0*(x*x - 1)*(y*y - 1)*(z+1)
00667                    //+ 0.0*(x*x - 1)*(z*z - 1)*(y-1) + 0.0*(x*x - 1)*(z*z - 1)*(y+1)
00668                    //+ 0.0*(y*y - 1)*(z*z - 1)*(x-1) + 0.0*(y*y - 1)*(z*z - 1)*(x+1)
00669                    //+ 0.0*(x*x - 1)*(y*y - 1)*(z*z - 1);
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   // the index in the Parent is very important
00947   // - in case of Lagrange first get the position of the evaluation
00948   // - then evaluate the refElement at that point
00949   // - see which local DoFs are not zero and return them
00950 
00951   Array<Array<Array<double> > >    result;
00952   const SpatialDerivSpecifier      deriv;
00953   int                              index = 0;
00954 
00955   //setVerb(6);
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:{ // 2D
00967     maximalCellType = QuadCell;
00968 
00969     switch (maxNrChild){
00970     case 9:{ //trisection on quads
00971       Point dofPos( 0.0 , 0.0 );
00972             // get the position of the DoF on the refCell
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       // shitf the position of the DOFs according to the position of the child in the parent cell
00982       switch (indexInParent){
00983         case 0: {/* nothing to add */  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:{ //todo: bisection on quads
01001        break;
01002        }
01003     }
01004   break;} // ------ end 2D -------
01005   case 3:{ //  ----- 3D ------------
01006 
01007     maximalCellType = BrickCell;
01008 
01009     Point dofPos( 0.0 , 0.0 , 0.0);
01010         // get the position of the DoF on the refCell
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     // shitf the position of the DOFs according to the position of the child in the parent cell
01020     switch (maxNrChild){
01021     case 27:{ // tri section on bricks
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         // 13 should never occur !!!!
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:{ // todo: bisection on bricks
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;} //end 3D
01064   }
01065 
01066     // evalute the results from the bi or trisection
01067   Array<double>& res = result[0][0];
01068   //SUNDANCE_OUT(this->verb() > 3 , "res=" << res );
01069 
01070   // get those data only once (since it is rather costly)
01071   if (doFInfromationCalculated_ == false)
01072   {
01073     getDoFsOrdered( maximalCellType, res.size() , facetD_, facetI_, facetN_);
01074     doFInfromationCalculated_ = true;
01075   }
01076 
01077   // here take the "result" array and get the nonzero elements
01078   index = 0;
01079   for ( int i=0 ; i < res.size() ; i++){
01080     // take only those DoFs which are non zero
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 /* -------------- Get the position of one DoF ----------------------- */
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:{  // Point DoFs
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:{  // Edge DoFs
01139     switch (maxCellDim){
01140     case 2:{
01141       // here we calculate the position of the DoF on the edge
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:{ // Edge in the brick test this
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:{ //face DoFs only in 3D
01170     // DoFs position on the face
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   // this could be done once at initialization
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     // found the next DoF and change the variables
01206     // nodo,indexo,dimo
01207     getNext = true;
01208     while ( (getNext) && (dimo < dofs_struct.size() )){
01209       //SUNDANCE_OUT(this->verb() > 3, "dimo=" << dimo << " ,indexo:" << indexo << ",nodo:" << nodo <<
01210       //    ",maxCellDim:" << maxCellDim << " ,nrDoF:" << nrDoF);
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               // actually we should not arrived here
01225               getNext = false;
01226             }
01227           }
01228          }
01229         }
01230         else
01231         {
01232                dimo++; indexo = 0; nodo = 0; getNext = false;
01233       }
01234       //SUNDANCE_OUT(this->verb() > 3, "dofs_struct.size()=" << dofs_struct.size() << " ,dofs_struct[dimo].size():"
01235       //    << dofs_struct[dimo].size() << ",dofs_struct[dimo][indexo].size():" << dofs_struct[dimo][indexo].size());
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 }
SundanceLagrange.cpp
