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 "SundanceCubicHermite.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 
00043 
00044 bool CubicHermite::supportsCellTypePair(
00045   const CellType& maximalCellType,
00046   const CellType& cellType
00047   ) const
00048 {
00049   switch(maximalCellType)
00050   {
00051     case LineCell:
00052       switch(cellType)
00053       {
00054         case LineCell:
00055         case PointCell:
00056           return true;
00057         default:
00058           return false;
00059       }
00060     case TriangleCell:
00061       switch(cellType)
00062       {
00063         case TriangleCell:
00064         case LineCell:
00065         case PointCell:
00066           return true;
00067         default:
00068           return false;
00069       }
00070 //     case TetCell:
00071 //       switch(cellType)
00072 //       {
00073 //         case TetCell:
00074 //         case TriangleCell:
00075 //         case LineCell:
00076 //         case PointCell:
00077 //           return true;
00078 //         default:
00079 //           return false;
00080 //       }
00081     default:
00082       return false;
00083   }
00084 }
00085 
00086 void CubicHermite::print(std::ostream& os) const 
00087 {
00088   os << "CubicHermite";
00089 }
00090 
00091 int CubicHermite::nReferenceDOFsWithoutFacets(
00092   const CellType& maximalCellType,
00093   const CellType& cellType
00094   ) const
00095 {
00096   switch(maximalCellType)
00097   {
00098     case LineCell:
00099       switch(cellType)
00100       {
00101         case PointCell:
00102           return 2;
00103         case LineCell:
00104           return 0;
00105         default:
00106           TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00107           return -1;
00108       }
00109       break;
00110     case TriangleCell:
00111       switch(cellType)
00112       {
00113         case PointCell:
00114           return 3;
00115         case LineCell:
00116           return 0;
00117         case TriangleCell:
00118           return 1;
00119         default:
00120           TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00121           return -1;
00122       }
00123       break;
00124     default:
00125       TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00126       return -1;
00127   }
00128   
00129 }
00130 
00131 void CubicHermite::getReferenceDOFs(
00132   const CellType& maximalCellType,
00133   const CellType& cellType,
00134   Array<Array<Array<int> > >& dofs) const 
00135 {
00136   typedef Array<int> Aint;
00137   switch(cellType)
00138   {
00139     case PointCell:
00140     {
00141       dofs.resize(1);
00142       if (maximalCellType==LineCell)
00143         dofs[0] = tuple<Aint>(tuple(0,1));
00144       else if (maximalCellType==TriangleCell)
00145         dofs[0] = tuple<Aint>(tuple(0,1,2));
00146       else TEUCHOS_TEST_FOR_EXCEPT(1);
00147       return;
00148     }
00149     break;
00150     case LineCell:
00151     {
00152       dofs.resize(2);
00153       dofs[0].resize(2);
00154       if (maximalCellType==LineCell)
00155       {
00156         dofs[0][0].resize(2);
00157         dofs[0][0][0] = 0;
00158         dofs[0][0][1] = 1;
00159         dofs[0][1].resize(2);
00160         dofs[0][1][0] = 2;
00161         dofs[0][1][1] = 3;
00162       }
00163       else if (maximalCellType==TriangleCell)
00164       {
00165         dofs[0][0].resize(3);
00166         dofs[0][0][0] = 0;
00167         dofs[0][0][1] = 1;
00168         dofs[0][0][2] = 2;
00169         dofs[0][1].resize(3);
00170         dofs[0][1][0] = 3;
00171         dofs[0][1][1] = 4;
00172         dofs[0][1][1] = 5;
00173       }
00174       else
00175       {
00176         TEUCHOS_TEST_FOR_EXCEPT(1);
00177       }
00178       dofs[1].resize(1);
00179       dofs[1][0].resize(0);
00180       return;
00181     }
00182     break;
00183     case TriangleCell:
00184     {
00185       dofs.resize(3);
00186       dofs[0].resize(3);
00187       dofs[0][0] = tuple(0,1,2);
00188       dofs[0][1] = tuple(3,4,5);
00189       dofs[0][2] = tuple(6,7,8);
00190       dofs[1].resize(3);
00191       dofs[1][0].resize(0);
00192       dofs[1][1].resize(0);
00193       dofs[1][2].resize(0);
00194       dofs[2].resize(1);
00195       dofs[2][0].resize(1);
00196       dofs[2][0][0] = 9;
00197     }
00198     break;
00199     default:
00200       TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00201         << cellType << " not implemented in CubicHermite basis");
00202   }
00203 }
00204 
00205 
00206 void CubicHermite::refEval(
00207   const CellType& cellType,
00208   const Array<Point>& pts,
00209   const SpatialDerivSpecifier& sds,
00210   Array<Array<Array<double> > >& result,
00211   int verbosity) const
00212 {
00213   TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()), 
00214     std::runtime_error,
00215     "cannot evaluate spatial derivative " << sds << " on CubicHermite basis");
00216   const MultiIndex& deriv = sds.mi();
00217   typedef Array<double> Adouble;
00218   result.resize(1);
00219   result[0].resize(pts.length());
00220 
00221   switch(cellType)
00222   {
00223     case PointCell:
00224       result[0] = tuple<Adouble>(tuple(1.0));
00225       return;
00226     case LineCell:
00227       for (int i=0; i<pts.length(); i++)
00228       {
00229         evalOnLine(pts[i], deriv, result[0][i]);
00230       }
00231       return;
00232     case TriangleCell:
00233       for (int i=0; i<pts.length(); i++)
00234       {
00235         evalOnTriangle(pts[i], deriv, result[0][i]);
00236       }
00237       return;
00238     default:
00239       TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00240         "CubicHermite::refEval() unimplemented for cell type "
00241         << cellType);
00242 
00243   }
00244 }
00245 
00246 /* ---------- evaluation on different cell types -------------- */
00247 
00248 void CubicHermite::evalOnLine(const Point& pt, 
00249   const MultiIndex& deriv,
00250   Array<double>& result) const
00251 {
00252   result.resize(4);
00253   ADReal x = ADReal(pt[0],0,1);
00254   Array<ADReal> tmp(4);
00255 
00256   tmp[0] = 1 + x * x * ( -3 + 2 * x );
00257   tmp[1] = x * ( 1 + (x - 2 ) * x );
00258   tmp[2] = ( 3 - 2*x ) * x * x;
00259   tmp[3] = (-1+x)*x*x;
00260 
00261   for (int i=0; i<tmp.length(); i++)
00262   {
00263     if (deriv.order()==0) 
00264     {
00265       result[i] = tmp[i].value();
00266     }
00267     else 
00268     {
00269       result[i] = tmp[i].gradient()[0];
00270     }
00271   }  
00272   return;
00273 }
00274 
00275 void CubicHermite::evalOnTriangle(const Point& pt, 
00276   const MultiIndex& deriv,
00277   Array<double>& result) const
00278 
00279 {
00280   result.resize(10);
00281   ADReal x = ADReal(pt[0], 0, 2);
00282   ADReal y = ADReal(pt[1], 1, 2);
00283   ADReal one(1.0, 2);
00284   
00285   Array<ADReal> tmp(10);
00286 
00287   SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00288     << y.value());
00289  
00290   tmp[0] = 1 - 3*x*x + 2 * x*x*x - 13*x*y + 13*x*x*y - 3*y*y + 13 *x*y*y + 2 *y*y*y;
00291   tmp[1] = x - 2 *x*x + x*x*x - 3*x*y + 3*x*x*y + 2*x*y*y;
00292   tmp[2] = y - 3 *x *y + 2* x*x* y - 2* y*y + 3* x*y*y + y*y*y;
00293   tmp[3] = 3 * x*x - 2* x*x*x - 7* x* y + 7* x*x *y + 7* x*y*y;
00294   tmp[4] = -x*x + x*x*x + 2*x *y - 2* x*x* y - 2* x* y*y;
00295   tmp[5] = -x* y + 2* x*x* y + x* y*y;
00296   tmp[6] = -7* x* y + 7* x*x*y + 3* y*y + 7* x* y*y - 2* y*y*y;
00297   tmp[7] = -x *y + x*x* y + 2* x* y*y;
00298   tmp[8] = 2 *x *y - 2* x*x* y - y*y - 2* x* y*y + y*y*y;
00299   tmp[9] = 27* x *y - 27* x*x* y - 27* x* y*y;
00300 
00301   for (int i=0; i<tmp.length(); i++)
00302   {
00303     if (deriv.order()==0) result[i] = tmp[i].value();
00304     else 
00305       result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00306   }
00307 }
00308 
00309 void CubicHermite::evalOnTet(const Point& pt, 
00310   const MultiIndex& deriv,
00311   Array<double>& result) const
00312 {
00313   ADReal x = ADReal(pt[0], 0, 3);
00314   ADReal y = ADReal(pt[1], 1, 3);
00315   ADReal z = ADReal(pt[2], 2, 3);
00316   ADReal one(1.0, 3);
00317   
00318   TEUCHOS_TEST_FOR_EXCEPT(true);
00319 }
00320 
00321 void CubicHermite::preApplyTransformation( const CellType &maxCellType ,
00322              const Mesh &mesh, 
00323              const Array<int> &cellLIDs,
00324              const CellJacobianBatch& JVol,
00325              RCP<Array<double> >& A ) const
00326   {
00327     switch(maxCellType)
00328       {
00329       case TriangleCell:
00330   preApplyTransformationTriangle( mesh , cellLIDs, JVol , A );
00331   break;
00332       default:
00333   TEUCHOS_TEST_FOR_EXCEPT(1);
00334   break;
00335       }
00336     return;
00337   }
00338 
00339 void CubicHermite::postApplyTransformation( const CellType &maxCellType ,
00340               const Mesh &mesh, 
00341               const Array<int> &cellLIDs,
00342               const CellJacobianBatch& JVol,
00343               RCP<Array<double> >& A ) const
00344   {
00345     switch(maxCellType)
00346       {
00347       case TriangleCell:
00348   postApplyTransformationTriangle( mesh , cellLIDs, JVol , A );
00349   break;
00350       default:
00351   TEUCHOS_TEST_FOR_EXCEPT(1);
00352   break;
00353       }
00354     return;
00355   }
00356 
00357 void CubicHermite::preApplyTransformationTranspose( const CellType &maxCellType ,
00358                 const Mesh &mesh, 
00359                 const Array<int> &cellLIDs,
00360                 const CellJacobianBatch& JVol ,
00361                 Array<double> & A ) const
00362 {
00363   switch(maxCellType)
00364     {
00365     case TriangleCell:
00366       preApplyTransformationTransposeTriangle( mesh , cellLIDs, JVol , A );
00367       break;
00368     default:
00369       TEUCHOS_TEST_FOR_EXCEPT(1);
00370       break;
00371     }
00372   return;
00373 }
00374 
00375 
00376 void CubicHermite::preApplyTransformationTriangle( const Mesh &mesh, 
00377                const Array<int> &cellLIDs,
00378                const CellJacobianBatch& JVol,
00379                RCP<Array<double> >& A) const
00380 {
00381   Array<double> &Anoptr = *A;
00382   
00383   Array<double> cellVertH;
00384   getVertexH( mesh , cellLIDs , cellVertH );
00385   
00386   
00387   // this applies M from the left on each cell
00388   // A for each cell has 10 rows because it is Hermite
00389   // so, this gives the number of columns per cell
00390   const int numCols = Anoptr.size() / JVol.numCells() / 10;
00391   
00392   for (int i=0;i<JVol.numCells();i++) 
00393     {
00394       const int cell_start = i * numCols * 10;
00395       const double *invJ = JVol.jVals(i);
00396       
00397       for (int j=0;j<numCols;j++) 
00398   {
00399     const int col_start = cell_start + 10 * j;
00400     const double a1 = Anoptr[col_start + 1];
00401     const double a2 = Anoptr[col_start + 2];
00402     const double a4 = Anoptr[col_start + 4];
00403     const double a5 = Anoptr[col_start + 5];
00404     const double a7 = Anoptr[col_start + 7];
00405     const double a8 = Anoptr[col_start + 8];
00406     Anoptr[col_start+1] = (invJ[0]*a1 + invJ[1]*a2)/cellVertH[3*i];
00407     Anoptr[col_start+2] = (invJ[2]*a1 + invJ[3]*a2)/cellVertH[3*i];
00408     Anoptr[col_start+4] = (invJ[0]*a4 + invJ[1]*a5)/cellVertH[3*i+1];
00409     Anoptr[col_start+5] = (invJ[2]*a4 + invJ[3]*a5)/cellVertH[3*i+1];
00410     Anoptr[col_start+7] = (invJ[0]*a7 + invJ[1]*a8)/cellVertH[3*i+2];
00411     Anoptr[col_start+8] = (invJ[2]*a7 + invJ[3]*a8)/cellVertH[3*i+2];
00412   }
00413     }
00414   
00415 }
00416 
00417 void CubicHermite::postApplyTransformationTriangle( const Mesh &mesh,
00418                 const Array<int> &cellLIDs , 
00419                 const CellJacobianBatch& JVol,
00420                 RCP<Array<double> >& A ) const
00421 {
00422   Array<double> &Anoptr = *A;
00423   
00424   Array<double> cellVertH;
00425   getVertexH( mesh , cellLIDs , cellVertH );
00426   
00427   
00428   const int numRows = Anoptr.size() / 10 / JVol.numCells();
00429   
00430   for (int i=0;i<JVol.numCells();i++) 
00431     {
00432       const double *invJ = JVol.jVals(i);
00433       
00434       const int cell_start = i * numRows * 10;
00435       // handle columns 1 and 2
00436       for (int j=0;j<numRows;j++) 
00437   {
00438     const double a = Anoptr[cell_start + numRows + j];
00439     const double b = Anoptr[cell_start + 2*numRows + j];
00440     Anoptr[cell_start + numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i];
00441     Anoptr[cell_start + 2*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i];
00442   }
00443       
00444       // handle columns 4 and 5
00445       for (int j=0;j<numRows;j++) 
00446   {
00447     const double a = Anoptr[cell_start + 4*numRows + j];
00448     const double b = Anoptr[cell_start + 5*numRows + j];
00449     Anoptr[cell_start + 4*numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i+1];
00450     Anoptr[cell_start + 5*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i+1];
00451   }
00452       
00453       // handle columns 7 and 8
00454       for (int j=0;j<numRows;j++) 
00455   {
00456     const double a = Anoptr[cell_start + 7*numRows + j];
00457     const double b = Anoptr[cell_start + 8*numRows + j];
00458     Anoptr[cell_start + 7*numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i+2];
00459     Anoptr[cell_start + 8*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i+2];
00460   }
00461       
00462     }
00463   
00464   return;
00465 }
00466 
00467 void CubicHermite::preApplyTransformationTransposeTriangle( const Mesh &mesh,
00468                   const Array<int> &cellLIDs,
00469                   const CellJacobianBatch& JVol,
00470                   Array<double>& A ) const
00471 {
00472   // this applies M from the left on each cell
00473   // A for each cell has 10 rows because it is Hermite
00474   // so, this gives the number of columns per cell
00475   const int numCols = A.size() / JVol.numCells() / 10;
00476   
00477   
00478   Array<double> cellVertH;
00479   getVertexH( mesh , cellLIDs , cellVertH );
00480   
00481   for (int i=0;i<JVol.numCells();i++) 
00482     {
00483       const int cell_start = i * numCols * 10;
00484       const double *invJ = JVol.jVals(i);
00485       
00486       for (int j=0;j<numCols;j++) 
00487   {
00488     const int col_start = cell_start + 10 * j;
00489     const double a1 = A[col_start + 1];
00490     const double a2 = A[col_start + 2];
00491     const double a4 = A[col_start + 4];
00492     const double a5 = A[col_start + 5];
00493     const double a7 = A[col_start + 7];
00494     const double a8 = A[col_start + 8];
00495     A[col_start+1] = (invJ[0]*a1 + invJ[2]*a2)/cellVertH[3*i];
00496     A[col_start+2] = (invJ[1]*a1 + invJ[3]*a2)/cellVertH[3*i];
00497     A[col_start+4] = (invJ[0]*a4 + invJ[2]*a5)/cellVertH[3*i+1];
00498     A[col_start+5] = (invJ[1]*a4 + invJ[3]*a5)/cellVertH[3*i+1];
00499     A[col_start+7] = (invJ[0]*a7 + invJ[2]*a8)/cellVertH[3*i+2];
00500     A[col_start+8] = (invJ[1]*a7 + invJ[3]*a8)/cellVertH[3*i+2];
00501   }
00502     }
00503   
00504 }
00505 
00506 // cellVertexH is 3 * cellLIDs.size and has the appropriate h for each
00507 // vertex of each cell
00508 void CubicHermite::getVertexH( const Mesh &mesh, 
00509              const Array<int> &cellLIDs ,
00510              Array<double> &cellVertexH ) const
00511 {
00512   cellVertexH.resize(3*cellLIDs.size());
00513   const int N = mesh.numCells(2);
00514   const double h = 1.0 / sqrtl( N );
00515   for (int i=0;i<cellVertexH.size();i++) 
00516     {
00517       cellVertexH[i] = h;
00518     }
00519 //     // now I get the vertices that are included in the batch of cells
00520 
00521 
00522 //     map<int,double> vert_h;
00523 
00524 //     Array<int> cellVertLIDs(3);
00525 //     Array<int> vertOrient(3);
00526 //     const int cellDim = mesh.spatialDim();
00527 //     // set up map structure
00528 //     for (int i=0;i<cellLIDs.size();i++) 
00529 //       {
00530 //  const int cellLID = cellLIDs[i];
00531 //  mesh.getFacetArray( cellDim , cellLID , 0 , cellVertLIDs , vertOrient );
00532 //  for (int j=0;j<3;j++)
00533 //    {
00534 //      map<int,double>::iterator vert_h_cur = vert_h.find(cellVertLIDs[j]);
00535 //      if (vert_h_cur == vert_h.end())
00536 //        {
00537 //    vert_h[cellVertLIDs[j]] = 0.0;
00538 //        }
00539 //    }
00540 //       }
00541     
00542 //     // fill in numbers
00543 //     for (map<int,double>::iterator it=vert_h.begin();it!=vert_h.end();++it)
00544 //       {
00545 //  const int vert = it->first;
00546 //  it->second = 0.0;
00547 //  // loop over maximal cofacets of each vertex
00548 //  Array<int> incident_cells;
00549 //  Array<double> cell_diams;
00550 //  mesh.getCofacets( 0 , vert , cellDim , incident_cells );
00551 //  mesh.getCellDiameters( 2 , incident_cells , cell_diams );
00552 //  for (int i=0;i<cell_diams.size();i++) 
00553 //    {
00554 //      it->second += cell_diams[i];
00555 //    }
00556 //  it->second /= cell_diams.size();
00557 
00558 //       }
00559 
00560 //     cellVertexH.resize( 3 * cellLIDs.size() );
00561 
00562 //     for (int i=0;i<cellLIDs.size();i++) 
00563 //       {
00564 //  const int cellLID = cellLIDs[i];
00565 //  mesh.getFacetArray( cellDim , cellLID , 0 , cellVertLIDs , vertOrient );
00566 //  for (int j=0;j<3;j++) 
00567 //    {
00568 //      cellVertexH[3*i+j] = 1./vert_h[cellVertLIDs[j]];
00569 //    }
00570 //       }
00571 
00572 
00573 
00574 
00575     return;
00576 
00577 }
SundanceCubicHermite.cpp
