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 "SundanceCurveEvalMediator.hpp"
00032 #include "SundanceCoordExpr.hpp"
00033 #include "SundanceTempStack.hpp"
00034 #include "SundanceCellDiameterExpr.hpp"
00035 #include "SundanceCurveNormExpr.hpp"
00036 #include "SundanceCellVectorExpr.hpp"
00037 #include "SundanceSpatialDerivSpecifier.hpp"
00038 #include "SundanceDiscreteFunction.hpp"
00039 #include "SundanceDiscreteFuncElement.hpp"
00040 #include "SundanceCellJacobianBatch.hpp"
00041 #include "SundanceOut.hpp"
00042 #include "PlayaTabs.hpp"
00043 #include "PlayaExceptions.hpp"
00044 #include "SundanceCurveIntegralCalc.hpp"
00045 
00046 #include "Teuchos_BLAS.hpp"
00047 
00048 
00049 using namespace Sundance;
00050 using namespace Teuchos;
00051 using namespace Playa;
00052 
00053 TEUCHOS_TIMER(coordEvalTimer, "Quad mediator: coord eval")
00054 
00055 CurveEvalMediator::CurveEvalMediator(
00056     const Mesh& mesh,   const ParametrizedCurve& paramcurve ,
00057     int cellDim, const QuadratureFamily& quad) : StdFwkEvalMediator(mesh, cellDim) ,
00058     numEvaluationCases_(-1),
00059     quad_(quad),
00060     numQuadPtsForMaxCell_(0) ,
00061     paramcurve_(paramcurve) {
00062 
00063   // set the cell types
00064   maxCellType_ = mesh.cellType(mesh.spatialDim() );
00065   curveCellType_ = mesh.cellType(mesh.spatialDim()-1);
00066 
00067   // detect the (constant) quadrature points (we must have the same nr quad points per cell!!!)
00068   Array<Point> quadPoints;
00069   Array<double> quadWeights;
00070   quad.getPoints(curveCellType_ , quadPoints , quadWeights );
00071   numQuadPtsForMaxCell_ = quadWeights.size();
00072 }
00073 
00074 Time& CurveEvalMediator::coordEvaluationTimer()
00075 {
00076   return coordEvalTimer();
00077 }
00078 
00079 void CurveEvalMediator::setCellType(const CellType& cellType,
00080   const CellType& maxCellType,
00081   bool isInternBdry) 
00082 {
00083   StdFwkEvalMediator::setCellType(cellType, maxCellType, isInternBdry);
00084 
00085   maxCellType_ = maxCellType;
00086 
00087   Tabs tab;
00088 
00089   SUNDANCE_MSG2(verb(), tab <<  "CurveEvalMediator::setCellType: cellType="
00090     << cellType << " cellDim=" << cellDim() << " maxCellType=" << maxCellType);
00091   SUNDANCE_MSG2(verb(), tab << "integration spec: =" 
00092     << integrationCellSpec());
00093   SUNDANCE_MSG2(verb(), tab << "forbid cofacet integrations =" 
00094     << forbidCofacetIntegrations());
00095   if (isInternalBdry()) 
00096   {
00097     SUNDANCE_MSG2(verb(), tab << "working on internal boundary");
00098   }
00099 
00100 
00101   if (cellType != maxCellType)
00102   {
00103     // this should not happen
00104   }
00105   else 
00106   {
00107     Tabs tab1;
00108     SUNDANCE_MSG2(verb(), tab1 << "no need for facet cases; work on original cell"); 
00109     numEvaluationCases_ = 1;
00110   }
00111 
00112   SUNDANCE_MSG2(verb(), tab << "num eval cases =" << numEvaluationCases_); 
00113 
00114 }
00115 
00116 int CurveEvalMediator::numQuadPts(const CellType& ct) const
00117 {
00118   if (ct != maxCellType_){
00119     // throw error
00120     TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "CurveEvalMediator::numQuadPts , cellType must be maxCellType_:" << ct);
00121   }
00122   return numQuadPtsForMaxCell_;
00123 }
00124 
00125 void CurveEvalMediator::evalCellDiameterExpr(const CellDiameterExpr* expr,
00126   RCP<EvalVector>& vec) const 
00127 {
00128   // this is the same for curve integral as for QuadratureEvalMediator
00129   Tabs tabs;
00130   SUNDANCE_MSG2(verb(),tabs 
00131     << "CurveEvalMediator evaluating cell diameter expr "
00132     << expr->toString());
00133 
00134   int nQuad = numQuadPts(cellType());
00135   int nCells = cellLID()->size();
00136 
00137   SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00138   Array<double> diameters;
00139   mesh().getCellDiameters(cellDim(), *cellLID(), diameters);
00140 
00141   vec->resize(nQuad*nCells);
00142   double * const xx = vec->start();
00143   int k=0;
00144   for (int c=0; c<nCells; c++)
00145   {
00146     double h = diameters[c];
00147     for (int q=0; q<nQuad; q++, k++) 
00148     {
00149       xx[k] = h;
00150     }
00151   }
00152 }
00153 
00154 void CurveEvalMediator::evalCellVectorExpr(const CellVectorExpr* expr,
00155   RCP<EvalVector>& vec) const 
00156 {
00157   Tabs tabs;
00158   SUNDANCE_MSG2(verb(),tabs 
00159     << "CurveEvalMediator evaluating cell vector expr "
00160     << expr->toString());
00161 
00162   int nQuad = numQuadPts(cellType());
00163   int nCells = cellLID()->size();
00164 
00165   vec->resize(nQuad*nCells);
00166 
00167   SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00168   int dir = expr->componentIndex();
00169 
00170   Array<Point> vectors;
00171   if (expr->isNormal())
00172   { 
00173     mesh().outwardNormals(*cellLID(), vectors);
00174   }
00175   else
00176   {
00177     TEUCHOS_TEST_FOR_EXCEPTION(cellDim() != 1, std::runtime_error,
00178       "unable to compute tangent vectors for cell dim = " << cellDim());
00179     mesh().tangentsToEdges(*cellLID(), vectors);
00180   }
00181     
00182   double * const xx = vec->start();
00183   int k=0;
00184   for (int c=0; c<nCells; c++)
00185   {
00186     double n = vectors[c][dir];
00187     for (int q=0; q<nQuad; q++, k++) 
00188     {
00189       xx[k] = n;
00190     }
00191   }
00192 }
00193 
00194 void CurveEvalMediator::evalCoordExpr(const CoordExpr* expr,
00195   RCP<EvalVector>& vec) const 
00196 {
00197   Tabs tabs;
00198   SUNDANCE_MSG2(verb(),tabs
00199     << "CurveEvalMediator evaluating coord expr "
00200     << expr->toString());
00201 
00202   TimeMonitor timer(coordEvalTimer());
00203   
00204   int nQuad = numQuadPts(cellType());
00205   int nCells = cellLID()->size();
00206   int d = expr->dir();
00207   
00208   SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad << ", nCells:" << nCells);
00209 
00210   // the size must be correct!!
00211   vec->resize(nQuad*nCells);
00212   double * const xx = vec->start();
00213   int k=0;
00214   int maxDim = mesh().spatialDim() ;
00215   Array<int> cellLIDs_tmp(1);
00216   Array<Point> phyPoints;
00217   Array<Point> refPoints;
00218   Array<Point> refDevs;
00219   Array<Point> refNormal;
00220 
00221   // iterate over each cell
00222   for (int c=0; c<nCells; c++)
00223   {
00224   int maxCellLID = (*cellLID())[c];
00225   cellLIDs_tmp[0] = (*cellLID())[c];
00226   // see if the mesh already contains the intersection points
00227   if ( mesh().hasCurvePoints( maxCellLID , paramcurve_.myID() ))
00228   {
00229     mesh().getCurvePoints( maxCellLID , paramcurve_.myID() , refPoints , refDevs , refNormal );
00230   }
00231   else // we have to calculate now the points
00232   {
00233     // calculate the intersection points
00234     CurveIntegralCalc::getCurveQuadPoints( maxCellType_ , maxCellLID , mesh() , paramcurve_ , quad_ ,
00235         refPoints, refDevs , refNormal);
00236     // store the intersection point in the mesh
00237     mesh().setCurvePoints( maxCellLID , paramcurve_.myID() , refPoints, refDevs , refNormal );
00238   }
00239 
00240   // calculate the physical coordinates
00241   mesh().pushForward( maxDim , cellLIDs_tmp , refPoints , phyPoints );
00242 
00243   // ---- return the desired component of the intersection point ----
00244     for (int q=0; q<nQuad; q++, k++)
00245     {
00246       xx[k] = phyPoints[q][d];
00247     }
00248   }
00249 }
00250 
00251 
00252 void CurveEvalMediator::evalCurveNormExpr(const CurveNormExpr* expr,
00253   RCP<EvalVector>& vec) const {
00254 
00255     Tabs tabs;
00256     SUNDANCE_MSG2(verb(),tabs
00257       << "CurveEvalMediator evaluating curve normal expr "
00258       << expr->toString());
00259 
00260     int nQuad = numQuadPts(cellType());
00261     int nCells = cellLID()->size();
00262     int d = expr->dir();
00263 
00264     SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00265 
00266     // the size of this vector must be correct !!!
00267     vec->resize(nCells*nQuad);
00268 
00269     double * const xx = vec->start();
00270     int k=0;
00271     Array<int> cellLIDs_tmp(1);
00272     Array<Point> phyPoints;
00273     Array<Point> refPoints;
00274     Array<Point> refDevs;
00275     Array<Point> refNormal;
00276 
00277     // iterate over each cell
00278     for (int c=0; c<nCells; c++)
00279     {
00280     int maxCellLID = (*cellLID())[c];
00281     cellLIDs_tmp[0] = (*cellLID())[c];
00282     // see if the mesh already contains the intersection points
00283     if ( mesh().hasCurvePoints( maxCellLID , paramcurve_.myID() ))
00284     {
00285       mesh().getCurvePoints( maxCellLID , paramcurve_.myID() , refPoints , refDevs , refNormal );
00286     }
00287     else // we have to calculate now the points
00288     {
00289       // calculate the intersection points
00290       CurveIntegralCalc::getCurveQuadPoints( maxCellType_ , maxCellLID , mesh() , paramcurve_ , quad_ ,
00291           refPoints, refDevs , refNormal);
00292       // store the intersection point in the mesh
00293       mesh().setCurvePoints( maxCellLID , paramcurve_.myID() , refPoints, refDevs , refNormal );
00294     }
00295 
00296     // ---- return the desired component of the unit normal vector ----
00297       for (int q=0; q<nQuad; q++, k++)
00298       {
00299         xx[k] = refNormal[q][d];
00300       }
00301     }
00302 }
00303 
00304 
00305 void CurveEvalMediator
00306 ::evalDiscreteFuncElement(const DiscreteFuncElement* expr,
00307   const Array<MultiIndex>& multiIndices,
00308   Array<RCP<EvalVector> >& vec) const
00309 {
00310 
00311     int verbo = dfVerb();
00312     Tabs tab1;
00313     SUNDANCE_MSG2(verbo , tab1
00314       << "CurveEvalMediator evaluating Discrete Function expr " << expr->toString());
00315 
00316     const DiscreteFunctionData* f = DiscreteFunctionData::getData(expr);
00317 
00318     TEUCHOS_TEST_FOR_EXCEPTION(f==0, std::logic_error,
00319       "QuadratureEvalMediator::evalDiscreteFuncElement() called "
00320       "with expr that is not a discrete function");
00321 
00322     SUNDANCE_MSG2(verbo , tab1 << "After casting DiscreteFunctionData" << expr->toString());
00323 
00324     RCP<Array<Array<double> > > localValues;
00325     Array<int> cellLIDs_tmp(1);
00326     Array<Point> phyPoints;
00327     Array<Point> refPoints;
00328     Array<Point> refDevs;
00329     Array<Point> refNormal;
00330     int nCells = cellLID()->size();
00331 
00332       RCP<const MapStructure> mapStruct;
00333     int myIndex = expr->myIndex();
00334     int nQuad = numQuadPtsForMaxCell_;
00335     Array<int> k(multiIndices.size(),0);
00336 
00337     Teuchos::BLAS<int,double> blas;
00338 
00339     SUNDANCE_MSG2(verbo , tab1 << "After declaring BLAS: " << expr->toString());
00340 
00341     // resize correctly the result vector
00342     for (int i=0; i<multiIndices.size(); i++)
00343     {
00344       vec[i]->resize(nCells*nQuad);
00345     }
00346 
00347     // loop over each cell
00348     for (int c=0; c<nCells; c++)
00349     {
00350       int maxCellLID = (*cellLID())[c];
00351         localValues = rcp(new Array<Array<double> >());
00352       SUNDANCE_MSG2(verbo , tab1 <<  "Cell:" << c <<  " of " << nCells << " , maxCellLID:" << maxCellLID );
00353 
00354       cellLIDs_tmp[0] = (*cellLID())[c];
00355       SUNDANCE_MSG2(verbo , tab1 <<  " Before calling f->getLocalValues:" << cellLIDs_tmp.size() <<
00356           " f==0 : " << (f==0) << " tmp:" << f->mesh().spatialDim());
00357       // - get local values from the DiscreteFunctionElementData
00358       mapStruct = f->getLocalValues(maxCellDim(), cellLIDs_tmp , *localValues);
00359       SUNDANCE_MSG2(verbo , tab1 <<  " After getting mapStruct:" << maxCellLID );
00360       SUNDANCE_MSG2(verbo , tab1 <<  " mapStruct->numBasisChunks():" << mapStruct->numBasisChunks() );
00361         int chunk = mapStruct->chunkForFuncID(myIndex);
00362         int funcIndex = mapStruct->indexForFuncID(myIndex);
00363         int nFuncs = mapStruct->numFuncs(chunk);
00364       SUNDANCE_MSG2(verbo , tab1 <<  " chunk:" << chunk );
00365       SUNDANCE_MSG2(verbo , tab1 <<  " funcIndex:" << funcIndex );
00366       SUNDANCE_MSG2(verbo , tab1 <<  " nFuncs:" << nFuncs );
00367 
00368       // the chunk of the function
00369         BasisFamily basis = rcp_dynamic_cast<BasisFamilyBase>(mapStruct->basis(chunk));
00370         int nNodesTotal = basis.nReferenceDOFsWithFacets(maxCellType(), maxCellType());
00371 
00372       // - get intersection (reference)points from the mesh (if not existent than compute them)
00373       if ( mesh().hasCurvePoints( maxCellLID , paramcurve_.myID() ))
00374       {
00375         mesh().getCurvePoints( maxCellLID , paramcurve_.myID() , refPoints , refDevs , refNormal );
00376       }
00377       else // we have to calculate now the points
00378       {
00379         // calculate the intersection points
00380         CurveIntegralCalc::getCurveQuadPoints( maxCellType_ , maxCellLID , mesh() , paramcurve_ , quad_ ,
00381             refPoints, refDevs , refNormal);
00382         // store the intersection point in the mesh
00383         mesh().setCurvePoints( maxCellLID , paramcurve_.myID() , refPoints, refDevs , refNormal );
00384       }
00385 
00386       // loop over each multi-index
00387       SUNDANCE_MSG2(verbo , tab1 << " multiIndices.size()" << multiIndices.size() );
00388       for (int i=0; i<multiIndices.size(); i++)
00389       {
00390 
00391         int nDerivResults = 1;
00392         if ( multiIndices[i].order() == 1 ) nDerivResults = maxCellDim();
00393 
00394           int pDir = 0;
00395           int derivNum = 1;
00396         MultiIndex mi;
00397         SUNDANCE_MSG2(verbo , tab1 << " before asking anything i = " << i);
00398         SUNDANCE_MSG2(verbo , tab1 << " multiindex order : " << multiIndices[i].order());
00399         SUNDANCE_MSG2(verbo , tab1 << " multiindex : " << multiIndices[i] );
00400         if (multiIndices[i].order() > 0){
00401           pDir = multiIndices[i].firstOrderDirection();
00402           mi[pDir] = 1;
00403           derivNum = mesh().spatialDim();
00404         }
00405 
00406         Array<Array<double> > result(nQuad*derivNum);
00407         Array<Array<Array<double> > > tmp;
00408 
00409         int offs = nNodesTotal;
00410 
00411         // resize the result vector
00412         for (int deriv = 0 ; deriv < derivNum ; deriv++)
00413         {
00414                     // test weather we have to compute derivative
00415           if (multiIndices[i].order() > 0){
00416             // in case of derivatives we set one dimension
00417             MultiIndex mi_tmp;
00418             mi_tmp[deriv] = 1;
00419             SpatialDerivSpecifier deriv(mi_tmp);
00420             SUNDANCE_MSG2(verbo , tab1 << "computing derivatives : " << deriv << " on reference cell ");
00421             basis.refEval( maxCellType_ , refPoints , deriv, tmp , verbo );
00422           } else {
00423             SpatialDerivSpecifier deriv(mi);
00424              // --- end eval basis functions
00425             SUNDANCE_MSG2(verbo , tab1 << "computing values reference cell ");
00426             basis.refEval( maxCellType_ , refPoints , deriv, tmp , verbo );
00427           }
00428 
00429           SUNDANCE_MSG2(verbo , tab1 << "resize result vector , offs:" << offs);
00430           for (int q=0; q<nQuad; q++){
00431             result[nQuad*deriv + q].resize(offs);
00432           }
00433 
00434             // copy the result in an other format
00435           SUNDANCE_MSG2(verbo , tab1 << "copy results ");
00436           int  offs1 = 0;
00437             for (int q=0; q<nQuad; q++)
00438             {
00439               offs1 = 0;
00440             for (int d=0; d<basis.dim(); d++)
00441             {
00442                int nNodes = tmp[d][q].size();
00443                  for (int n=0; n<nNodes; n++ , offs1++ )
00444                  {
00445                       result[nQuad*deriv + q][offs1] = tmp[d][q][n];
00446                  }
00447               }
00448           }
00449         }// loop over all dimensional derivative
00450 
00451                 // multiply the local results with the coefficients, (matrix vector OP)
00452         SUNDANCE_MSG2(verbo , tab1 << "summing up values , funcIndex:" << funcIndex << " offs:" << offs);
00453         for (int deriv = 0 ; deriv < derivNum ; deriv++)
00454         {
00455           for (int q=0; q<nQuad; q++)
00456           {
00457             double sum = 0.0;
00458             // sum over nodes
00459             for (int n = 0 ; n < offs ; n++){
00460               sum = sum + result[nQuad*deriv + q][n] * (*localValues)[chunk][funcIndex*offs + n];
00461             }
00462             // sum up the result in the 0th element
00463             result[nQuad*deriv + q][0] = sum;
00464           }
00465         }
00466 
00467           // multiply the result if necesary with the inverse of the Jacobian
00468         const CellJacobianBatch& J = JTrans();
00469         if (mi.order()==1)
00470         {
00471           Tabs tab1;
00472           Tabs tab2;
00473             SUNDANCE_MSG2(verbo, tab2 << "Jacobian batch nCells=" << J.numCells());
00474               SUNDANCE_MSG2(verbo, tab2 << "Jacobian batch cell dim=" << J.cellDim());
00475               SUNDANCE_MSG2(verbo, tab2 << "Jacobian batch spatial dim=" << J.spatialDim());
00476                     // we just multiply the derivative direction component
00477               Array<double> invJ;
00478               J.getInvJ(c, invJ);
00479           for (int q=0; q<nQuad; q++)
00480           {
00481             double sum = 0.0;
00482             for (int deriv = 0 ; deriv < derivNum ; deriv++)
00483             {
00484               // multiply one row from the J^{-T} matrix with the gradient vector
00485               sum = sum + result[nQuad*deriv + q][0] * invJ[derivNum*pDir + deriv];
00486             }
00487             // the resulting derivative on the physical cell in the "pDir" direction
00488             result[q][0] = sum;
00489           }
00490           }
00491 
00492         // --- just copy the result to the "vec" back, the result should be in the  "result[q][0]" place----
00493         //SUNDANCE_MSG2(verbo , tab1 << "copy results back ");
00494         double* vecPtr = vec[i]->start();
00495         for (int q=0; q<nQuad; q++, k[i]++)
00496         {
00497           vecPtr[k[i]] = result[q][0];
00498         }
00499         SUNDANCE_MSG2(verbo , tab1 << " END copy results back ");
00500       } // --- end loop multiindex
00501       SUNDANCE_MSG2(verbo , tab1 << " END loop over multiindex ");
00502   }// --- end loop over cells
00503   SUNDANCE_MSG2(verbo , tab1 << " END loop over cells ");
00504 }
00505 
00506 
00507 void CurveEvalMediator::print(std::ostream& os) const
00508 {
00509   // todo: implement this
00510 }
00511 
00512
SundanceCurveEvalMediator.cpp
