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 MeshBase::MeshBase(int dim, const MPIComm& comm, 
00045   const MeshEntityOrder& order) 
00046   : dim_(dim), 
00047     comm_(comm),
00048     order_(order),
00049     reorderer_(Mesh::defaultReorderer().createInstance(this)),
00050     validWeights_(true),
00051     specialWeights_(),
00052     curvePoints_Are_Valid_(true),
00053     nrCurvesForIntegral_(0),
00054     curvePoints_() ,
00055     curveDerivative_(),
00056     curveNormal_(),
00057     curveID_to_ArrayIndex_()
00058 { specialWeights_.resize(dim_);
00059   curvePoints_.resize(0);
00060   curveDerivative_.resize(0);
00061   curveNormal_.resize(0);
00062 }
00063 
00064 
00065 
00066 Point MeshBase::centroid(int cellDim, int cellLID) const
00067 {
00068   if (cellDim==0) return nodePosition(cellLID);
00069   int dummy;
00070   Point x = nodePosition(facetLID(cellDim, cellLID, 0, 0, dummy));
00071   int nf = numFacets(cellDim, cellLID, 0);
00072   for (int f=1; f<nf; f++) 
00073     x += nodePosition(facetLID(cellDim, cellLID, 0, f, dummy));
00074   return x / ((double) nf);
00075 }
00076 
00077 void MeshBase::outwardNormals(
00078   const Array<int>& cellLIDs,
00079   Array<Point>& outwardNormals
00080   ) const 
00081 {
00082   int D = spatialDim();
00083   outwardNormals.resize(cellLIDs.size());
00084   for (int c=0; c<cellLIDs.size(); c++)
00085   {
00086     int f=-1;
00087     TEUCHOS_TEST_FOR_EXCEPTION(numMaxCofacets(D-1, cellLIDs[c]) > 1, 
00088       std::runtime_error,
00089       "cell #" << cellLIDs[c] << " is not a boundary cell");
00090     int maxLID = maxCofacetLID(D-1, cellLIDs[c], 0, f);
00091     Point cInterior = centroid(D, maxLID);
00092     Point cBdry = centroid(D-1, cellLIDs[c]);
00093     Point q = cBdry - cInterior;
00094     Point s;
00095     if (D==1) 
00096     {
00097       s = Point(1.0);
00098     }
00099     else if (D==2)
00100     {
00101       Point A = nodePosition(facetLID(D-1, cellLIDs[c], 0, 0, f));
00102       Point B = nodePosition(facetLID(D-1, cellLIDs[c], 0, 1, f));
00103       Point t = B - A;
00104       s = Point(-t[1], t[0]);
00105     }
00106     else 
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 C = nodePosition(facetLID(D-1, cellLIDs[c], 0, 2, f));
00111       s = cross(B-A, C-A);
00112     }
00113     if (q * s > 0.0)
00114     {
00115       outwardNormals[c] = s/::sqrt(s*s);
00116     }
00117     else
00118     {
00119       outwardNormals[c] = -s/::sqrt(s*s);
00120     }
00121   }
00122 }
00123 
00124 
00125 void  MeshBase::tangentsToEdges(
00126   const Array<int>& cellLIDs,
00127   Array<Point>& tangentVectors
00128   ) const 
00129 {
00130   TEUCHOS_TEST_FOR_EXCEPT(spatialDim() <= 1);
00131 
00132   tangentVectors.resize(cellLIDs.size());
00133 
00134   for (int c=0; c<cellLIDs.size(); c++)
00135   {
00136     int fOrient=1;
00137     Point A = nodePosition(facetLID(1, cellLIDs[c], 0, 0, fOrient));
00138     Point B = nodePosition(facetLID(1, cellLIDs[c], 0, 1, fOrient));
00139     Point t = B - A;
00140     tangentVectors[c] = t/(sqrt(t*t));
00141   }
00142 }
00143 
00144 
00145 
00146 
00147 void MeshBase::getFacetArray(int cellDim, int cellLID, int facetDim, 
00148   Array<int>& facetLIDs,
00149   Array<int>& facetOrientations) const
00150 {
00151   int nf = numFacets(cellDim, cellLID, facetDim);
00152   facetLIDs.resize(nf);
00153   facetOrientations.resize(nf);
00154   for (int f=0; f<nf; f++) 
00155   {
00156     facetLIDs[f] = facetLID(cellDim, cellLID, facetDim, f, 
00157       facetOrientations[f]);
00158   }
00159 }
00160 
00161 
00162 void MeshBase::getLabels(int cellDim, const Array<int>& cellLID, 
00163   Array<int>& labels) const
00164 {
00165   labels.resize(cellLID.size());
00166   for (int i=0; i<cellLID.size(); i++) labels[i] = label(cellDim, cellLID[i]);
00167 }
00168 
00169 
00170 // ===================== storing special weights for specfic cells  ======================
00171 
00172 bool MeshBase::hasSpecialWeight(int dim, int cellLID) const {
00173   if (specialWeights_[dim-1].containsKey(cellLID))
00174         return ((specialWeights_[dim-1].get(cellLID).size() > 0));
00175   else
00176     return false;
00177 }
00178 
00179 void MeshBase::setSpecialWeight(int dim, int cellLID, Array<double>& w) const {
00180   specialWeights_[dim-1].put(cellLID,w);
00181 }
00182 
00183 void MeshBase::getSpecialWeight(int dim, int cellLID, Array<double>& w) const {
00184   w = specialWeights_[dim-1].get(cellLID);
00185 }
00186 
00187 /** deletes all special weights so those have to be recreated*/
00188 void MeshBase::flushSpecialWeights() const {
00189   // delete all the special weights for all cells
00190   for (int d = 0 ; d < dim_ ; d++){
00191     // reset the map to a new and empty one
00192     specialWeights_[d] = Sundance::Map< int , Array<double> >();
00193   }
00194   validWeights_ = true;
00195 }
00196 
00197 // ===================== storing curve intersection/quadrature points ======================
00198 
00199 void MeshBase::flushCurvePoints() const {
00200   // all the curve points should be deleted
00201   // first reset the Map between curve ID and index
00202   curveID_to_ArrayIndex_ = Sundance::Map< int , int >();
00203   nrCurvesForIntegral_ = 0;
00204   curvePoints_Are_Valid_ = true;
00205   for (int c = 0 ; c < curvePoints_.size() ; c++ ){
00206     // here delete each Map which existed before
00207     curvePoints_[c] = Sundance::Map< int , Array<Point> >();
00208     curveDerivative_[c] = Sundance::Map< int , Array<Point> >();
00209     curveNormal_[c] = Sundance::Map< int , Array<Point> >();
00210   }
00211 }
00212 
00213 bool MeshBase::hasCurvePoints(int maxCellLID , int curveID) const {
00214     if (curvePoints_[mapCurveID_to_Index(curveID)].containsKey(maxCellLID))
00215     return ( curvePoints_[mapCurveID_to_Index(curveID)].get(maxCellLID).size() > 0 );
00216     else
00217       return false;
00218 }
00219 
00220 void MeshBase::setCurvePoints(int maxCellLID, int curveID ,
00221     Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const {
00222 
00223     Tabs tabs;
00224     int verbo = 0;
00225     SUNDANCE_MSG3(verbo, tabs << "MeshBase::setCurvePoints , nr:" << nrCurvesForIntegral_);
00226       curvePoints_[mapCurveID_to_Index(curveID)].put( maxCellLID , points );
00227       curveDerivative_[mapCurveID_to_Index(curveID)].put( maxCellLID , derivs );
00228       curveNormal_[mapCurveID_to_Index(curveID)].put( maxCellLID , normals );
00229 
00230 }
00231 
00232 void MeshBase::getCurvePoints(int maxCellLID, int curveID ,
00233     Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const {
00234 
00235     Tabs tabs;
00236     int verbo = 0;
00237     SUNDANCE_MSG3(verbo, tabs << "MeshBase::getCurvePoints , nr:" << nrCurvesForIntegral_);
00238       points = curvePoints_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00239       derivs = curveDerivative_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00240       normals = curveNormal_[mapCurveID_to_Index(curveID)].get( maxCellLID );
00241 
00242 }
00243 
00244 int MeshBase::mapCurveID_to_Index(int curveID) const {
00245 
00246    Tabs tabs;
00247    int verbo = 0;
00248 
00249    SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index curveID:" << curveID);
00250      if (curveID_to_ArrayIndex_.containsKey(curveID)){
00251        SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index value found for curveID:" << curveID << " ret:" << curveID_to_ArrayIndex_.get(curveID));
00252          return curveID_to_ArrayIndex_.get(curveID);
00253      } else {
00254        SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index create new :" << nrCurvesForIntegral_);
00255          curveID_to_ArrayIndex_.put( curveID , nrCurvesForIntegral_ );
00256          SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index , increment ");
00257          nrCurvesForIntegral_++;
00258        curvePoints_.resize(nrCurvesForIntegral_);
00259        curveDerivative_.resize(nrCurvesForIntegral_);
00260        curveNormal_.resize(nrCurvesForIntegral_);
00261        SUNDANCE_MSG3(verbo, tabs << "MeshBase::mapCurveID_to_Index create new :" << nrCurvesForIntegral_);
00262          return nrCurvesForIntegral_-1;
00263      }
00264 
00265 }
SundanceMeshBase.cpp
