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 "SundanceElementIntegral.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "Teuchos_TimeMonitor.hpp"
00035 
00036 using namespace Sundance;
00037 using namespace Teuchos;
00038 
00039 using std::endl;
00040 
00041 static Time& transCreationTimer() 
00042 {
00043   static RCP<Time> rtn 
00044     = TimeMonitor::getNewTimer("building integral transformation matrices"); 
00045   return *rtn;
00046 }
00047 
00048 
00049 bool& ElementIntegral::alwaysUseCofacets()
00050 {
00051   static bool rtn = true;
00052   return rtn;
00053 }
00054 
00055 ElementIntegral::ElementIntegral(int spatialDim,
00056   const CellType& maxCellType,
00057   int dim, 
00058   const CellType& cellType,
00059   bool isInternalBdry,
00060   const ParametrizedCurve& globalCurve,
00061   const Mesh& mesh,
00062   int verb)
00063   : setupVerb_(verb),
00064     integrationVerb_(0),
00065     transformVerb_(0),
00066     spatialDim_(spatialDim),
00067     dim_(dim),
00068     isInternalBdry_(isInternalBdry),
00069     nFacetCases_(1),
00070     testDerivOrder_(-1), 
00071     nRefDerivTest_(-1),
00072     nNodesTest_(-1),
00073     unkDerivOrder_(-1), 
00074     nRefDerivUnk_(-1),
00075     nNodesUnk_(-1),
00076     nNodes_(-1),
00077     order_(0),
00078     alpha_(),
00079     beta_(),
00080     cellType_(cellType),
00081     maxCellType_(maxCellType),
00082     evalCellType_(cellType),
00083     testBasis_(),
00084     unkBasis_(),
00085     globalCurve_(globalCurve),
00086     mesh_(mesh)
00087 {
00088   Tabs tab0;
00089   SUNDANCE_MSG2(setupVerb(), tab0 << "constructing 0-form ElementIntegral");
00090   /* if we're integrating a derivative along a facet, we need to refer back
00091    * to the maximal cell. */
00092   if (alwaysUseCofacets() || (dim != spatialDim && !isInternalBdry))
00093   {
00094     evalCellType_ = maxCellType_;
00095     nFacetCases_ = numFacets(maxCellType, dim);
00096   }
00097 }
00098 
00099 ElementIntegral::ElementIntegral(int spatialDim,
00100   const CellType& maxCellType,
00101   int dim, 
00102   const CellType& cellType,
00103   const BasisFamily& testBasis,
00104   int alpha,
00105   int testDerivOrder,
00106   bool isInternalBdry,
00107   const ParametrizedCurve& globalCurve,
00108   const Mesh& mesh,
00109   int verb)
00110   : setupVerb_(verb),
00111     integrationVerb_(0),
00112     transformVerb_(0),
00113     spatialDim_(spatialDim),
00114     dim_(dim),
00115     isInternalBdry_(isInternalBdry),
00116     nFacetCases_(1),
00117     testDerivOrder_(testDerivOrder), 
00118     nRefDerivTest_(ipow(spatialDim, testDerivOrder)),
00119     nNodesTest_(testBasis.nReferenceDOFsWithFacets(maxCellType, cellType)),
00120     unkDerivOrder_(-1), 
00121     nRefDerivUnk_(-1),
00122     nNodesUnk_(-1),
00123     nNodes_(nNodesTest_),
00124     order_(1),
00125     alpha_(alpha),
00126     beta_(-1),
00127     cellType_(cellType),
00128     maxCellType_(maxCellType),
00129     evalCellType_(cellType),
00130     testBasis_(testBasis),
00131     unkBasis_(),
00132     globalCurve_(globalCurve),
00133     mesh_(mesh)
00134 {
00135   Tabs tab0(0);
00136   SUNDANCE_MSG2(setupVerb(), tab0 << "constructing 1-form ElementIntegral");
00137   /* if we're integrating a derivative along a facet, we 
00138    * may need to refer back to the maximal cell. */
00139   bool okToRestrictTestToBdry = basisRestrictableToBoundary(testBasis);
00140     
00141   Tabs tab1;
00142   SUNDANCE_MSG2(setupVerb(), tab1 << "dim=" << dim << " spatialDim=" << spatialDim);
00143   if (dim != spatialDim)
00144   {
00145     if (isInternalBdry)
00146     {
00147       TEUCHOS_TEST_FOR_EXCEPT(!okToRestrictTestToBdry);
00148     }
00149     if (alwaysUseCofacets() || testDerivOrder>0)
00150     {
00151       Tabs tab2;
00152       evalCellType_ = maxCellType_;
00153       nFacetCases_ = numFacets(maxCellType, dim);
00154       nNodesTest_ = testBasis.nReferenceDOFsWithFacets(maxCellType, maxCellType);
00155       SUNDANCE_MSG2(setupVerb(), tab2 << "nNodesTest=" << nNodesTest_);
00156       nNodes_ = nNodesTest_;
00157       TEUCHOS_TEST_FOR_EXCEPT(nNodes_ == 0);
00158     }
00159     else
00160     {
00161       TEUCHOS_TEST_FOR_EXCEPT(!okToRestrictTestToBdry);
00162     }
00163   }
00164 
00165   SUNDANCE_MSG2(setupVerb(), tab1 << "nNodes=" << nNodes_);
00166 }
00167 
00168 
00169 
00170 ElementIntegral::ElementIntegral(int spatialDim,
00171   const CellType& maxCellType,
00172   int dim,
00173   const CellType& cellType,
00174   const BasisFamily& testBasis,
00175   int alpha,
00176   int testDerivOrder,
00177   const BasisFamily& unkBasis,
00178   int beta,
00179   int unkDerivOrder,
00180   bool isInternalBdry,
00181   const ParametrizedCurve& globalCurve,
00182   const Mesh& mesh,
00183   int verb)
00184   : setupVerb_(verb),
00185     integrationVerb_(0),
00186     transformVerb_(0),
00187     spatialDim_(spatialDim),
00188     dim_(dim),
00189     isInternalBdry_(isInternalBdry),
00190     nFacetCases_(1),
00191     testDerivOrder_(testDerivOrder), 
00192     nRefDerivTest_(ipow(spatialDim, testDerivOrder)),
00193     nNodesTest_(testBasis.nReferenceDOFsWithFacets(maxCellType, cellType)), 
00194     unkDerivOrder_(unkDerivOrder), 
00195     nRefDerivUnk_(ipow(spatialDim, unkDerivOrder)),
00196     nNodesUnk_(unkBasis.nReferenceDOFsWithFacets(maxCellType, cellType)), 
00197     nNodes_(nNodesTest_*nNodesUnk_),
00198     order_(2),
00199     alpha_(alpha),
00200     beta_(beta),
00201     cellType_(cellType),
00202     maxCellType_(maxCellType),
00203     evalCellType_(cellType),
00204     testBasis_(testBasis),
00205     unkBasis_(unkBasis),
00206     globalCurve_(globalCurve),
00207     mesh_(mesh)
00208 {
00209   Tabs tab0(0);
00210   SUNDANCE_MSG2(setupVerb(), tab0 
00211     << "***** constructing 2-form ElementIntegral");
00212   /* if we're integrating a derivative along a facet, we may need to refer back
00213    * to the maximal cell. */
00214   bool okToRestrictTestToBdry = basisRestrictableToBoundary(testBasis);
00215   bool okToRestrictUnkToBdry = basisRestrictableToBoundary(unkBasis);
00216 
00217     
00218   Tabs tab1;
00219   SUNDANCE_MSG2(setupVerb(), tab1 << "dim=" << dim << " spatialDim=" << spatialDim);
00220   if (dim != spatialDim)
00221   {
00222     if (isInternalBdry)
00223     {
00224       TEUCHOS_TEST_FOR_EXCEPT(!(okToRestrictTestToBdry && okToRestrictUnkToBdry));   
00225     }
00226     if (alwaysUseCofacets() || testDerivOrder>0 || unkDerivOrder>0)
00227     {
00228       Tabs tab2;
00229       evalCellType_ = maxCellType_;
00230       nFacetCases_ = numFacets(maxCellType, dim);
00231       nNodesTest_ = testBasis.nReferenceDOFsWithFacets(maxCellType, maxCellType);
00232       SUNDANCE_MSG2(setupVerb(), tab2 << "nNodesTest=" << nNodesTest_);
00233       nNodesUnk_ = unkBasis.nReferenceDOFsWithFacets(maxCellType, maxCellType);
00234       SUNDANCE_MSG2(setupVerb(), tab2 << "nNodesUnk=" << nNodesUnk_);
00235       nNodes_ = nNodesTest_ * nNodesUnk_;
00236       TEUCHOS_TEST_FOR_EXCEPT(nNodes_ == 0);
00237     }
00238     else
00239     {
00240       TEUCHOS_TEST_FOR_EXCEPT(okToRestrictTestToBdry != okToRestrictUnkToBdry);
00241     }
00242   }
00243 
00244   SUNDANCE_MSG2(setupVerb(), tab1 << "nNodes=" << nNodes_);
00245 }
00246 
00247 
00248 void ElementIntegral::setVerb(
00249   int integrationVerb,
00250   int transformVerb)
00251 {
00252   integrationVerb_ = integrationVerb;
00253   transformVerb_ = transformVerb;
00254 }
00255 
00256 void ElementIntegral::describe(std::ostream& os) const 
00257 {
00258   Tabs tab(0);
00259   bool hasTest = testDerivOrder() >= 0;
00260   bool hasUnk = unkDerivOrder() >= 0;
00261 
00262   if (hasTest)
00263   {
00264     os << tab << "Test functions:" << std::endl;
00265     Tabs tab1;
00266     os << tab1 << "test basis = " << testBasis() << std::endl;
00267     os << tab1 << "test differentiation order = " << testDerivOrder() << std::endl;
00268     os << tab1 << "alpha = " << alpha() << std::endl;
00269     os << tab1 << "num test gradient components=" << nRefDerivTest() << std::endl;
00270   }
00271   if (hasUnk)
00272   {
00273     os << tab << "Unknown functions:" << std::endl;
00274     Tabs tab1;
00275     os << tab1 << "unk basis = " << unkBasis() << std::endl;
00276     os << tab1 << "unk differentiation order = " << unkDerivOrder() << std::endl;
00277     os << tab1 << "beta = " << beta() << std::endl;
00278     os << tab1  << "num unk gradient components=" << nRefDerivUnk() << std::endl;
00279   }
00280   os << tab << "Geometry:" << std::endl;
00281   Tabs tab1;
00282   os << tab1 << "cell type on which integral is defined: " << cellType()
00283      << std::endl;
00284   os << tab1 << "maximal cell type: " << maxCellType() << std::endl;
00285   os << tab1 << "cell type on which quad pts defined: " << evalCellType()
00286      << std::endl;
00287   os << tab << "Number of evaluation cases: " << nFacetCases() << std::endl;
00288 }
00289 
00290 
00291 void ElementIntegral::assertBilinearForm() const 
00292 {
00293   TEUCHOS_TEST_FOR_EXCEPTION(testDerivOrder() < 0 || testDerivOrder() > 1,
00294     std::logic_error,
00295     "Test function derivative order=" << testDerivOrder()
00296     << " must be 0 or 1");
00297   
00298   TEUCHOS_TEST_FOR_EXCEPTION(unkDerivOrder() < 0 || unkDerivOrder() > 1,
00299     std::logic_error,
00300     "Unknown function derivative order=" << unkDerivOrder()
00301     << " must be 0 or 1");
00302 }
00303 
00304 void ElementIntegral::assertLinearForm() const 
00305 {
00306   TEUCHOS_TEST_FOR_EXCEPTION(testDerivOrder() < 0 || testDerivOrder() > 1,
00307     std::logic_error,
00308     "Test function derivative order=" << testDerivOrder()
00309     << " must be 0 or 1");
00310 }
00311 
00312 
00313 void ElementIntegral::getQuad(const QuadratureFamily& quad,
00314   int evalCase, Array<Point>& quadPts, Array<double>& quadWeights) const 
00315 {
00316   Tabs tab(0);
00317 
00318   SUNDANCE_MSG2(setupVerb(), tab << "getting quad points for rule "
00319     << quad);
00320 
00321   if (nFacetCases()==1) 
00322   {
00323     quad.getPoints(cellType(), quadPts, quadWeights);
00324   }
00325   else 
00326   {
00327     int dim = dimension(cellType());
00328     quad.getFacetPoints(maxCellType(), dim, evalCase, 
00329       quadPts, quadWeights);
00330   }
00331 
00332   if (setupVerb() >= 4)
00333   {
00334     Tabs tab1;
00335     SUNDANCE_MSG4(setupVerb(), 
00336       tab1 << "quadrature points on ref cell are:");
00337     printQuad(Out::os(), quadPts, quadWeights);
00338   }
00339 }
00340   
00341 
00342 
00343 Array<double>& ElementIntegral::G(int alpha)
00344 {
00345   static Array<Array<double> > rtn(3);
00346 
00347   return rtn[alpha];
00348 }
00349 
00350 Array<double>& ElementIntegral::G(int alpha, int beta)
00351 {
00352   static Array<Array<double> > rtn(9);
00353 
00354   return rtn[3*alpha + beta];
00355 }
00356 
00357 int& ElementIntegral::transformationMatrixIsValid(int alpha, int beta)
00358 {
00359   static Array<int> rtn(9, false);
00360   return rtn[3*alpha + beta];
00361 }
00362 
00363 int& ElementIntegral::transformationMatrixIsValid(int alpha)
00364 {
00365   static Array<int> rtn(3, false);
00366   return rtn[alpha];
00367 }
00368 
00369 void ElementIntegral::invalidateTransformationMatrices()
00370 {
00371   for (int i=0; i<3; i++)
00372   {
00373     transformationMatrixIsValid(i) = false;
00374     for (int j=0; j<3; j++)
00375     {
00376       transformationMatrixIsValid(i, j) = false;
00377     }
00378   }
00379 }
00380 
00381 
00382 
00383 
00384 int ElementIntegral::ipow(int base, int power) 
00385 {
00386   int rtn = 1;
00387   for (int i=0; i<power; i++) rtn *= base;
00388   return rtn;
00389 }
00390 
00391 void ElementIntegral
00392 ::createTwoFormTransformationMatrix(const CellJacobianBatch& JTrans,
00393   const CellJacobianBatch& JVol) const
00394 {
00395   TimeMonitor timer(transCreationTimer());
00396   Tabs tab;
00397 
00398   int flops = 0;
00399 
00400   int maxDim = JTrans.cellDim();
00401   //int cellDim = JVol.cellDim();
00402 
00403   if (testDerivOrder() == 1 && unkDerivOrder() == 1)
00404   {
00405     Tabs tab2;
00406     if (transformationMatrixIsValid(alpha(), beta())) return;
00407     transformationMatrixIsValid(alpha(), beta()) = true;
00408 
00409     G(alpha(), beta()).resize(JTrans.numCells() * JTrans.cellDim() * JTrans.cellDim());
00410 
00411     double* GPtr = &(G(alpha(),beta())[0]);
00412     int k = 0;
00413 
00414     for (int c=0; c<JTrans.numCells(); c++)
00415     {
00416       static Array<double> invJ;
00417       JTrans.getInvJ(c, invJ);
00418       double detJ = fabs(JVol.detJ()[c]);
00419       for (int gamma=0; gamma<maxDim; gamma++)
00420       {
00421         for (int delta=0; delta<maxDim; delta++, k++)
00422         {
00423           GPtr[k] =  detJ*invJ[alpha() + gamma*maxDim]
00424             * invJ[beta() + maxDim*delta];
00425         }
00426       }
00427     }
00428     flops = 2 * JTrans.numCells() * maxDim * maxDim + JTrans.numCells();
00429   }
00430 
00431   else if (testDerivOrder() == 1 && unkDerivOrder() == 0)
00432   {
00433     if (transformationMatrixIsValid(alpha())) return;
00434     transformationMatrixIsValid(alpha()) = true;
00435 
00436     G(alpha()).resize(JTrans.numCells() * JTrans.cellDim());
00437 
00438     int k = 0;
00439     double* GPtr = &(G(alpha())[0]);
00440 
00441     for (int c=0; c<JTrans.numCells(); c++)
00442     {
00443       static Array<double> invJ;
00444       JTrans.getInvJ(c, invJ);
00445       double detJ = fabs(JVol.detJ()[c]);
00446       for (int gamma=0; gamma<maxDim; gamma++,k++)
00447       {
00448         GPtr[k] = detJ*invJ[alpha() + maxDim * gamma];
00449       }
00450     }
00451     flops = JTrans.numCells() * maxDim + JTrans.numCells();
00452   }
00453 
00454   else 
00455   {
00456     if (transformationMatrixIsValid(beta())) return;
00457     transformationMatrixIsValid(beta()) = true;
00458 
00459     G(beta()).resize(JTrans.numCells() * JTrans.cellDim());
00460 
00461     int k = 0;
00462     double* GPtr = &(G(beta())[0]);
00463 
00464     for (int c=0; c<JTrans.numCells(); c++)
00465     {
00466       static Array<double> invJ;
00467       JTrans.getInvJ(c, invJ);
00468       double detJ = fabs(JVol.detJ()[c]);
00469       for (int gamma=0; gamma<maxDim; gamma++,k++)
00470       {
00471         GPtr[k] = detJ*invJ[beta() + maxDim * gamma];
00472       }
00473     }
00474     flops = JTrans.numCells() * maxDim + JTrans.numCells();
00475   }
00476 
00477   addFlops(flops);
00478 }
00479 
00480 
00481 void ElementIntegral
00482 ::createOneFormTransformationMatrix(const CellJacobianBatch& JTrans,
00483   const CellJacobianBatch& JVol) const 
00484 {
00485   TimeMonitor timer(transCreationTimer());
00486   Tabs tab;
00487   SUNDANCE_MSG2(transformVerb(), 
00488     tab << "ElementIntegral creating linear form trans matrices");
00489 
00490   int maxDim = JTrans.cellDim();
00491 
00492   if (transformationMatrixIsValid(alpha())) return;
00493   transformationMatrixIsValid(alpha()) = true;
00494 
00495   int flops = JTrans.numCells() * maxDim + JTrans.numCells();
00496 
00497   G(alpha()).resize(JTrans.numCells() * JTrans.cellDim());
00498 
00499   int k = 0;
00500   double* GPtr = &(G(alpha())[0]);
00501 
00502   for (int c=0; c<JTrans.numCells(); c++)
00503   {
00504     Array<double> invJ;
00505     JTrans.getInvJ(c, invJ);
00506     double detJ = fabs(JVol.detJ()[c]);
00507     for (int gamma=0; gamma<maxDim; gamma++, k++)
00508     {
00509       GPtr[k] = detJ*invJ[alpha() + maxDim * gamma]; 
00510     }
00511   }
00512   
00513   addFlops(flops);
00514 }
00515
SundanceElementIntegral.cpp
