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 "SundanceLegendre.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 Teuchos;
00041 
00042 // ctor
00043 Legendre::Legendre(int order):order_(order)
00044 {
00045 
00046   // set the nr DOFs
00047 
00048   if (order_ > 1) nrDOF_edge_ = order_ - 1;
00049   else nrDOF_edge_ = 0;
00050 
00051   if (order_ > 3) nrDOF_face_ = ((order_-2)*(order_-3))/2;
00052   else nrDOF_face_ = 0;
00053 
00054   nrDOF_brick_ = 0;
00055 
00056   //SUNDANCE_OUT( true , "Legendre ctor nrDOF_edge_:" << nrDOF_edge_ << " , nrDOF_face_:" << nrDOF_face_ << " , nrDOF_brick_:"<<nrDOF_brick_);
00057 }
00058 
00059 bool Legendre::supportsCellTypePair(
00060   const CellType& maximalCellType,
00061   const CellType& cellType
00062   ) const
00063 {
00064   switch(maximalCellType)
00065   {
00066     case LineCell:
00067       switch(cellType)
00068       {
00069         case LineCell:
00070         case PointCell:
00071           return true;
00072         default:
00073           return false;
00074       }
00075     case QuadCell:
00076       switch(cellType)
00077       {
00078         case QuadCell:
00079         case LineCell:
00080         case PointCell:
00081           return true;
00082         default:
00083           return false;
00084       }
00085      case BrickCell:
00086        switch(cellType)
00087        {
00088          case BrickCell:
00089          case QuadCell:
00090          case LineCell:
00091          case PointCell:
00092            return true;
00093          default:
00094            return false;
00095       }
00096     default:
00097       return false;
00098   }
00099 }
00100 
00101 void Legendre::print(std::ostream& os) const
00102 {
00103   os << "Legendre";
00104 }
00105 
00106 int Legendre::nReferenceDOFsWithoutFacets(
00107   const CellType& maximalCellType,
00108   const CellType& cellType
00109   ) const
00110 {
00111       switch(cellType)
00112       {
00113         case PointCell:
00114           return 1;
00115         case LineCell:
00116           return nrDOF_edge_;
00117         case QuadCell:
00118             return nrDOF_face_;
00119         case BrickCell:
00120             return nrDOF_brick_;
00121         default:
00122             TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00123             return -1;
00124       }
00125 }
00126 
00127 Array<int> Legendre::makeRange(int low, int high)
00128 {
00129   if (high < low) return Array<int>();
00130 
00131   Array<int> rtn(high-low+1);
00132   for (int i=0; i<rtn.length(); i++) rtn[i] = low+i;
00133   return rtn;
00134 }
00135 
00136 void Legendre::getReferenceDOFs(
00137   const CellType& maximalCellType,
00138   const CellType& cellType,
00139   Array<Array<Array<int> > >& dofs) const 
00140 {
00141 
00142   typedef Array<int> Aint;
00143 
00144   //switch(order_)
00145 
00146   switch(cellType)
00147   {
00148     case PointCell:
00149     {
00150         dofs.resize(1);
00151         dofs[0] = tuple<Aint>(tuple(0));
00152         return;
00153     }
00154     break;
00155     case LineCell:
00156     {
00157         dofs.resize(2);
00158         dofs[0] = tuple<Aint>(tuple(0), tuple(1));
00159         if (nrDOF_edge_>0)
00160         {
00161           dofs[1] = tuple<Aint>(makeRange(2, 2+nrDOF_edge_-1));
00162         }
00163         else
00164         {
00165           dofs[1] = tuple(Array<int>());
00166         }
00167       return;
00168     }
00169     break;
00170     case QuadCell:
00171     {
00172       int offs = 0;
00173         dofs.resize(3);
00174         // dofs[0] are DOFs at Points
00175         dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2), tuple(3));
00176         offs = 4;
00177         if (nrDOF_edge_>0)
00178         {
00179           dofs[1].resize(4);
00180           dofs[1][0] = makeRange(offs, offs+nrDOF_edge_-1);  offs += nrDOF_edge_;
00181           dofs[1][1] = makeRange(offs, offs+nrDOF_edge_-1);  offs += nrDOF_edge_;
00182           dofs[1][2] = makeRange(offs, offs+nrDOF_edge_-1);  offs += nrDOF_edge_;
00183           dofs[1][3] = makeRange(offs, offs+nrDOF_edge_-1);  offs += nrDOF_edge_;
00184         }
00185         else
00186         {
00187           dofs[1] = tuple(Array<int>(), Array<int>(),
00188                    Array<int>(), Array<int>());
00189         }
00190 
00191         if (nrDOF_edge_>0)
00192         {
00193           dofs[2].resize(1);
00194           dofs[2][0] = makeRange(offs, offs+nrDOF_face_-1);  offs += nrDOF_face_;
00195         }
00196         else
00197         {
00198           dofs[2] = tuple(Array<int>());
00199         }
00200         //SUNDANCE_OUT( true , "Legendre::getReferenceDOFs offs:" << offs );
00201     }
00202     break;
00203     default:
00204       TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00205         << cellType << " not implemented in Legendre basis");
00206   }
00207 }
00208 
00209 
00210 void Legendre::refEval(
00211   const CellType& cellType,
00212   const Array<Point>& pts,
00213   const SpatialDerivSpecifier& sds,
00214   Array<Array<Array<double> > >& result,
00215   int verbosity) const
00216 {
00217   TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()), 
00218     std::runtime_error,
00219     "cannot evaluate spatial derivative " << sds << " on Legendre basis");
00220   const MultiIndex& deriv = sds.mi();
00221   typedef Array<double> Adouble;
00222   result.resize(1);
00223   result[0].resize(pts.length());
00224 
00225   switch(cellType)
00226   {
00227     case PointCell:
00228       result[0] = tuple<Adouble>(tuple(1.0));
00229       return;
00230     case LineCell:
00231       for (int i=0; i<pts.length(); i++)
00232       {
00233         evalOnLine(pts[i], deriv, result[0][i]);
00234       }
00235       return;
00236     case QuadCell:
00237       for (int i=0; i<pts.length(); i++)
00238       {
00239         evalOnQuad(pts[i], deriv, result[0][i]);
00240       }
00241       return;
00242     case BrickCell:
00243       for (int i=0; i<pts.length(); i++)
00244       {
00245         evalOnBrick(pts[i], deriv, result[0][i]);
00246       }
00247       return;
00248     default:
00249       TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00250         "Legendre::refEval() unimplemented for cell type "
00251         << cellType);
00252 
00253   }
00254 }
00255 
00256 /* ---------- evaluation on different cell types -------------- */
00257 
00258 void Legendre::evalOnLine(const Point& pt,
00259   const MultiIndex& deriv,
00260   Array<double>& result) const
00261 {
00262   result.resize(2+nrDOF_edge_);
00263   ADReal x = ADReal(pt[0],0,1);
00264 
00265   Array<ADReal> refAll(7);
00266 
00267   refAll[0] = 1-x;
00268   refAll[1] = x;
00269   refAll[2] = 2.44948974278318 * ( (2*x-1)*(2*x-1) - 1 ) / 4;
00270   refAll[3] = 3.16227766016838 * ( (2*x-1)*(2*x-1) - 1 ) * (2*x-1) / 4;
00271   refAll[4] = 3.74165738677394 * ( 5*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 6*(2*x-1)*(2*x-1) + 1) / 16;
00272   refAll[5] = 4.24264068711929 * (2*x-1) * (7*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 10*(2*x-1)*(2*x-1) + 3) / 16;
00273   refAll[6] = 4.69041575982343 * (21*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) -
00274                               35*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) + 15*(2*x-1)*(2*x-1) - 1) / 32;
00275 
00276 
00277   for (int i=0; i<result.length(); i++)
00278   {
00279     if (deriv.order()==0) 
00280     {
00281       result[i] = refAll[i].value();
00282     }
00283     else 
00284     {
00285       result[i] = refAll[i].gradient()[0];
00286     }
00287   }  
00288   //SUNDANCE_OUT( true , "Legendre::evalOnLine result.length():" << result.length() );
00289   return;
00290 }
00291 
00292 void Legendre::evalOnQuad(const Point& pt,
00293   const MultiIndex& deriv,
00294   Array<double>& result) const
00295 
00296 {
00297   result.resize( 4 + 4*nrDOF_edge_ + nrDOF_face_);
00298   ADReal x = ADReal(pt[0], 0, 2);
00299   ADReal y = ADReal(pt[1], 1, 2);
00300   ADReal one(1.0, 2);
00301   
00302   Array<ADReal> refAllx(7);
00303   Array<ADReal> refAlly(7);
00304 
00305   refAllx[0] = 1-x;
00306   refAllx[1] = x;
00307   refAllx[2] = 2.44948974278318 * ( (2*x-1)*(2*x-1) - 1 ) / 4;
00308   refAllx[3] = 3.16227766016838 * ( (2*x-1)*(2*x-1) - 1 ) * (2*x-1) / 4;
00309   refAllx[4] = 3.74165738677394 * ( 5*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 6*(2*x-1)*(2*x-1) + 1) / 16;
00310   refAllx[5] = 4.24264068711929 * (2*x-1) * (7*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 10*(2*x-1)*(2*x-1) + 3) / 16;
00311   refAllx[6] = 4.69041575982343 * (21*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) -
00312                               35*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) + 15*(2*x-1)*(2*x-1) - 1) / 32;
00313 
00314   refAlly[0] = 1-y;
00315   refAlly[1] = y;
00316   refAlly[2] = 2.44948974278318 * ( (2*y-1)*(2*y-1) - 1 ) / 4;
00317   refAlly[3] = 3.16227766016838 * ( (2*y-1)*(2*y-1) - 1 ) * (2*y-1) / 4;
00318   refAlly[4] = 3.74165738677394 * ( 5*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) - 6*(2*y-1)*(2*y-1) + 1) / 16;
00319   refAlly[5] = 4.24264068711929 * (2*y-1) * (7*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) - 10*(2*y-1)*(2*y-1) + 3) / 16;
00320   refAlly[6] = 4.69041575982343 * (21*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) -
00321                               35*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) + 15*(2*y-1)*(2*y-1) - 1) / 32;
00322 
00323   SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00324     << y.value());
00325 
00326   int sideIndex[4][2] = { {0,0} , {1,0} , {0,1} , {1,1}};
00327   int edgeIndex[4]    = { 0 , 1 , 1 , 0};
00328   int edgeMultI[4]    = { 0 , 0 , 1 , 1};
00329   int offs = 0;
00330   Array<ADReal> tmp(4 + 4*nrDOF_edge_ + nrDOF_face_);
00331 
00332   // loop over vertexes
00333   for (int i=0; i < 4 ; i++, offs++){
00334     tmp[offs] = refAllx[sideIndex[i][0]] * refAlly[sideIndex[i][1]];
00335   }
00336 
00337   // loop over edges
00338   for (int i=0; i < 4 ; i++){
00339     // loop over each DOF on the edge
00340     if (edgeIndex[i] == 0){
00341       for (int j = 0 ; j < nrDOF_edge_ ; j++ , offs++){
00342         tmp[offs] = refAllx[2+j] * refAlly[edgeMultI[i]];
00343       }
00344     }
00345     else
00346     {
00347       for (int j = 0 ; j < nrDOF_edge_ ; j++ , offs++){
00348         tmp[offs] = refAllx[edgeMultI[i]] * refAlly[2+j];
00349       }
00350     }
00351   }
00352 
00353   // loop over all internal DOFs
00354   if ( nrDOF_face_ > 0 ){
00355     // loop for each hierarchical layer
00356     for (int hierarch = 0 ; hierarch < order_ - 3 ; hierarch++)
00357     {
00358       //SUNDANCE_OUT( true , "Legendre::evalOnQuad hierarch:" << hierarch );
00359       // for each layer add the basis function
00360       for (int i=0 ; i < hierarch+1 ; i++ , offs++)
00361       {
00362         //SUNDANCE_OUT( true , "Legendre::evalOnQuad offs:" << offs << " 2+i:" << 2+i << " , 2+(hierarch-1-i):" << 2+(hierarch-i));
00363         tmp[offs] = refAllx[2+i] * refAlly[2+(hierarch-i)];
00364       }
00365     }
00366   }
00367 
00368   // compute the results
00369   for (int i=0; i<result.length(); i++)
00370   {
00371     if (deriv.order()==0) result[i] = tmp[i].value();
00372     else 
00373       result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00374   }
00375   //SUNDANCE_OUT( true , "Legendre::evalOnQuad result.length():" << result.length() );
00376 }
00377 
00378 void Legendre::evalOnBrick(const Point& pt,
00379   const MultiIndex& deriv,
00380   Array<double>& result) const
00381 {
00382   ADReal x = ADReal(pt[0], 0, 3);
00383   ADReal y = ADReal(pt[1], 1, 3);
00384   ADReal z = ADReal(pt[2], 2, 3);
00385   ADReal one(1.0, 3);
00386   
00387   TEUCHOS_TEST_FOR_EXCEPT(true);
00388 }
SundanceLegendre.cpp
