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 "SundanceMeshBase.hpp"
00032 #include "SundanceMesh.hpp"
00033 #include "PlayaExceptions.hpp"
00034 #include "Teuchos_Time.hpp"
00035 #include "Teuchos_TimeMonitor.hpp"
00036 #include "PlayaTabs.hpp"
00037 
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Teuchos;
00041 using namespace Sundance;
00042 
00043 
00044 
00045 MeshBase::MeshBase(int dim, const MPIComm& comm, 
00046   const MeshEntityOrder& order) 
00047   : dim_(dim), 
00048     comm_(comm),
00049     order_(order),
00050     reorderer_(Mesh::defaultReorderer().createInstance(this)),
00051     validWeights_(true),
00052     specialWeights_(),
00053     curvePoints_Are_Valid_(true),
00054     nrCurvesForIntegral_(0),
00055     curvePoints_() ,
00056     curveDerivative_(),
00057     curveNormal_(),
00058     curveID_to_ArrayIndex_()
00059 { specialWeights_.resize(dim_);
00060   curvePoints_.resize(0);
00061   curveDerivative_.resize(0);
00062   curveNormal_.resize(0);
00063 }
00064 
00065 const int* MeshBase::elemZeroFacetView(int maxCellLID) const
00066 {
00067   TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00068     "MeshBase::elemZeroFacetView() is deprecated");
00069   return (const int*) 0;
00070 }
00071 
00072 
00073 Point MeshBase::centroid(int cellDim, int cellLID) const
00074 {
00075   if (cellDim==0) return nodePosition(cellLID);
00076   int dummy;
00077   Point x = nodePosition(facetLID(cellDim, cellLID, 0, 0, dummy));
00078   int nf = numFacets(cellDim, cellLID, 0);
00079   for (int f=1; f<nf; f++) 
00080     x += nodePosition(facetLID(cellDim, cellLID, 0, f, dummy));
00081   return x / ((double) nf);
00082 }
00083 
00084 void MeshBase::outwardNormals(
00085   const Array<int>& cellLIDs,
00086   Array<Point>& outwardNormals
00087   ) const 
00088 {
00089   int D = spatialDim();
00090   outwardNormals.resize(cellLIDs.size());
00091   for (int c=0; c<cellLIDs.size(); c++)
00092   {
00093     int f=-1;
00094     TEUCHOS_TEST_FOR_EXCEPTION(numMaxCofacets(D-1, cellLIDs[c]) > 1, 
00095       std::runtime_error,
00096       "cell #" << cellLIDs[c] << " is not a boundary cell");
00097     int maxLID = maxCofacetLID(D-1, cellLIDs[c], 0, f);
00098     Point cInterior = centroid(D, maxLID);
00099     Point cBdry = centroid(D-1, cellLIDs[c]);
00100     Point q = cBdry - cInterior;
00101     Point s;
00102     if (D==1) 
00103     {
00104       s = Point(1.0);
00105     }
00106     else if (D==2)
00107     {
00108       Point A = nodePosition(facetLID(D-1, cellLIDs[c], 0, 0, f));
00109       Point B = nodePosition(facetLID(D-1, cellLIDs[c], 0, 1, f));
00110       Point t = B - A;
00111       s = Point(-t[1], t[0]);
00112     }
00113     else 
00114     {
00115       Point A = nodePosition(facetLID(D-1, cellLIDs[c], 0, 0, f));
00116       Point B = nodePosition(facetLID(D-1, cellLIDs[c], 0, 1, f));
00117       Point C = nodePosition(facetLID(D-1, cellLIDs[c], 0, 2, f));
00118       s = cross(B-A, C-A);
00119     }
00120     if (q * s > 0.0)
00121     {
00122       outwardNormals[c] = s/::sqrt(s*s);
00123     }
00124     else
00125     {
00126       outwardNormals[c] = -s/::sqrt(s*s);
00127     }
00128   }
00129 }
00130 
00131 
00132 void  MeshBase::tangentsToEdges(
00133   const Array<int>& cellLIDs,
00134   Array<Point>& tangentVectors
00135   ) const 
00136 {
00137   TEUCHOS_TEST_FOR_EXCEPT(spatialDim() <= 1);
00138 
00139   tangentVectors.resize(cellLIDs.size());
00140 
00141   for (int c=0; c<cellLIDs.size(); c++)
00142   {
00143     int fOrient=1;
00144     Point A = nodePosition(facetLID(1, cellLIDs[c], 0, 0, fOrient));
00145     Point B = nodePosition(facetLID(1, cellLIDs[c], 0, 1, fOrient));
00146     Point t = B - A;
00147     tangentVectors[c] = t/(sqrt(t*t));
00148   }
00149 }
00150 
00151 
00152 
00153 
00154 void MeshBase::getFacetArray(int cellDim, int cellLID, int facetDim, 
00155   Array<int>& facetLIDs,
00156   Array<int>& facetOrientations) const
00157 {
00158   int nf = numFacets(cellDim, cellLID, facetDim);
00159   facetLIDs.resize(nf);
00160   facetOrientations.resize(nf);
00161   for (int f=0; f<nf; f++) 
00162   {
00163     facetLIDs[f] = facetLID(cellDim, cellLID, facetDim, f, 
00164       facetOrientations[f]);
00165   }
00166 }
00167 
00168 
00169 void MeshBase::getLabels(int cellDim, const Array<int>& cellLID, 
00170   Array<int>& labels) const
00171 {
00172   labels.resize(cellLID.size());
00173   for (int i=0; i<cellLID.size(); i++) labels[i] = label(cellDim, cellLID[i]);
00174 }
00175 
00176 
00177 void MeshBase::getMaxCofacetLIDs(
00178   const Array<int>& cellLIDs,
00179   MaximalCofacetBatch& cofacets) const
00180 {
00181   TEUCHOS_TEST_FOR_EXCEPT(cellLIDs.size()==0);
00182   int d = spatialDim() - 1;
00183   int nc = numMaxCofacets(d, cellLIDs[0]);  
00184 
00185   cofacets.reset(cellLIDs.size(), nc);
00186 
00187   for (int i=0; i<cellLIDs.size(); i++) 
00188   {
00189     int f1;
00190     int cofacet1 = maxCofacetLID(d, cellLIDs[i], 0, f1);
00191     if (nc==1) cofacets.addSingleCofacet(i, cofacet1, f1);
00192     else
00193     {
00194       int f2;
00195       int cofacet2 = maxCofacetLID(d, cellLIDs[i], 1, f2);
00196       cofacets.addTwoCofacets(i, cofacet1, f1, cofacet2, f2);
00197     }
00198   }
00199 }
00200 
00201 void MeshBase::getLIDsForLabel(int cellDim, int labelToCheck, 
00202   Array<int>& cellLIDs) const
00203 {
00204   cellLIDs.resize(0);
00205   int nc = numCells(cellDim); 
00206 
00207   for (int c=0; c<nc; c++)
00208   {
00209     int lc = label(cellDim, c);
00210     if (lc==labelToCheck) cellLIDs.append(c);
00211   }
00212 }
00213 
00214 Set<int> MeshBase::getAllLabelsForDimension(int cellDim) const
00215 {
00216   int nc = numCells(cellDim); 
00217   Set<int> rtn;
00218 
00219   for (int c=0; c<nc; c++)
00220   {
00221     rtn.put(label(cellDim, c));
00222   }
00223   return rtn;
00224 }
00225 
00226 // ===================== storing special weights for specfic cells  ======================
00227 
00228 bool MeshBase::hasSpecialWeight(int dim, int cellLID) const {
00229   if (specialWeights_[dim-1].containsKey(cellLID))
00230         return ((specialWeights_[dim-1].get(cellLID).size() > 0));
00231   else
00232     return false;
00233 }
00234 
00235 void MeshBase::setSpecialWeight(int dim, int cellLID, Array<double>& w) const {
00236   specialWeights_[dim-1].put(cellLID,w);
00237 }
00238 
00239 void MeshBase::getSpecialWeight(int dim, int cellLID, Array<double>& w) const {
00240   w = specialWeights_[dim-1].get(cellLID);
00241 }
00242 
00243 /** deletes all special weights so those have to be recreated*/
00244 void MeshBase::flushSpecialWeights() const {
00245   // delete all the special weights for all cells
00246   for (int d = 0 ; d < dim_ ; d++){
00247     // reset the map to a new and empty one
00248     specialWeights_[d] = Sundance::Map< int , Array<double> >();
00249   }
00250   validWeights_ = true;
00251 }
00252 
00253 // ===================== storing curve intersection/quadrature points ======================
00254 
00255 void MeshBase::flushCurvePoints() const {
00256   // all the curve points should be deleted
00257   // first reset the Map between curve ID and index
00258   curveID_to_ArrayIndex_ = Sundance::Map< int , int >();
00259   nrCurvesForIntegral_ = 0;
00260   curvePoints_Are_Valid_ = true;
00261   for (int c = 0 ; c < curvePoints_.size() ; c++ ){
00262     // here delete each Map which existed before
00263     curvePoints_[c] = Sundance::Map< int , Array<Point> >();
00264     curveDerivative_[c] = Sundance::Map< int , Array<Point> >();
00265     curveNormal_[c] = Sundance::Map< int , Array<Point> >();
00266   }
00267 }
00268 
00269 bool MeshBase::hasCurvePoints(int maxCellLID , int curveID) const {
00270     if (curvePoints_[mapCurveID_to_Index(curveID)].containsKey(maxCellLID))
00271     return ( curvePoints_[mapCurveID_to_Index(curveID)].get(maxCellLID).size() > 0 );
00272     else
00273       return false;
00274 }
00275 
00276 void MeshBase::setCurvePoints(int maxCellLID, int curveID ,
00277     Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const {
00278 
00279     Tabs tabs;
00280     int verbo = 0;
00281     SUNDANCE_MSG3(verbo, tabs << "MeshBase::setCurvePoints , nr:" << nrCurvesForIntegral_);
00282       curvePoints_[mapCurveID_to_Index(curveID)].put( maxCellLID , points );
00283       curveDerivative_[mapCurveID_to_Index(curveID)].put( maxCellLID , derivs );
00284       curveNormal_[mapCurveID_to_Index(curveID)].put( maxCellLID , normals );
00285 
00286 }
00287 
00288 void MeshBase::getCurvePoints(int maxCellLID, int curveID ,
00289     Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const {
00290 
00291     Tabs tabs;
00292     int verbo = 0;
00293     SUNDANCE_MSG3(verbo, tabs << "MeshBase::getCurvePoints , nr:" << nrCurvesForIntegral_);
00294       points = curvePoints_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00295       derivs = curveDerivative_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00296       normals = curveNormal_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00297 
00298 }
00299 
00300 int MeshBase::mapCurveID_to_Index(int curveID) const {
00301 
00302    Tabs tabs;
00303    int verbo = 0;
00304 
00305    SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index curveID:" << curveID);
00306      if (curveID_to_ArrayIndex_.containsKey(curveID)){
00307        SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index value found for curveID:" << curveID << " ret:" << curveID_to_ArrayIndex_.get(curveID));
00308          return curveID_to_ArrayIndex_.get(curveID);
00309      } else {
00310        SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index create new :" << nrCurvesForIntegral_);
00311          curveID_to_ArrayIndex_.put( curveID , nrCurvesForIntegral_ );
00312          SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index , increment ");
00313          nrCurvesForIntegral_++;
00314        curvePoints_.resize(nrCurvesForIntegral_);
00315        curveDerivative_.resize(nrCurvesForIntegral_);
00316        curveNormal_.resize(nrCurvesForIntegral_);
00317        SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index create new :" << nrCurvesForIntegral_);
00318          return nrCurvesForIntegral_-1;
00319      }
00320 
00321 }
SundanceMeshBase.cpp
