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 #ifndef SUNDANCE_MESH_H
00032 #define SUNDANCE_MESH_H
00033 
00034 #include "SundanceDefs.hpp"
00035 #include "SundanceMeshBase.hpp"
00036 #include "SundanceMaximalCofacetBatch.hpp"
00037 #include "SundanceIncrementallyCreatableMesh.hpp"
00038 #include "SundanceIdentityReorderer.hpp"
00039 #include "SundanceCellReorderer.hpp"
00040 #include "PlayaHandle.hpp"
00041 
00042 namespace Sundance
00043 {
00044 using namespace Teuchos;
00045 
00046 
00047 /**
00048  * Mesh is the user-level object representing discrete geometry. 
00049  * The Mesh class is a handle to a MeshBase, which is an abstract interface
00050  * for meshes. 
00051  */
00052 class Mesh : public Playa::Handle<MeshBase>
00053 {
00054 public:
00055 
00056   /* */
00057   HANDLE_CTORS(Mesh, MeshBase);
00058 
00059   /** */
00060   int id() const {return ptr()->id();}
00061 
00062   /** \brief Get the ordering convention used by this mesh */
00063   const MeshEntityOrder& meshOrder() const {return ptr()->meshOrder();}
00064     
00065   /** 
00066    * Get the spatial dimension of the mesh
00067    */
00068   int spatialDim() const {return ptr()->spatialDim();}
00069 
00070   /** 
00071    * Get the number of cells having dimension dim
00072    */
00073   int numCells(int dim) const {return ptr()->numCells(dim);}
00074 
00075   /** 
00076    * Return the position of the i-th node
00077    */
00078   Point nodePosition(int i) const {return ptr()->nodePosition(i);}
00079 
00080   /** 
00081    * Return a view of the i-th node's position
00082    */
00083   const double* nodePositionView(int i) const {return ptr()->nodePositionView(i);}
00084 
00085   /** 
00086    * Return the centroid position of the cellLID-th cell of dimension
00087    * cellDim.
00088    */
00089   Point centroid(int cellDim, int cellLID) const 
00090     {return ptr()->centroid(cellDim, cellLID);}
00091 
00092 
00093   /** 
00094    * Get the outward normals for the batch of cells of dimension
00095    * spatialDim()-1. If any cell in the batch is not on the boundary,
00096    * an exception is thrown. 
00097    *
00098    * \param cellLIDs [in] LIDs for the cells whose normals are to be
00099    * computed. 
00100    * \param outwardNormals [out] Outward normal unit vectors for each
00101    * cell in the batch.
00102    */
00103   void outwardNormals(
00104     const Array<int>& cellLIDs,
00105     Array<Point>& outwardNormals
00106     ) const 
00107     {ptr()->outwardNormals(cellLIDs, outwardNormals);}
00108 
00109   /** 
00110    * Get tangent vectors for a batch of edges
00111    *
00112    * \param cellLIDs [in] LIDs for the cells whose tangents are to be
00113    * computed. 
00114    * \param tangentVectors [out] Unit tangents for each cell
00115    */
00116   void tangentsToEdges(
00117     const Array<int>& cellLIDs,
00118     Array<Point>& tangentVectors
00119     ) const 
00120     {ptr()->tangentsToEdges(cellLIDs, tangentVectors);}
00121       
00122       
00123       
00124 
00125   /** 
00126    * Compute the jacobians of a batch of cells, returning the 
00127    * result via reference argument
00128    *
00129    * @param cellDim dimension of the cells whose Jacobians are to
00130    * be computed
00131    * @param cellLID local indices of the cells for which Jacobians
00132    * are to be computed
00133    * @param jBatch reference to the resulting Jacobian batch
00134    */
00135   void getJacobians(int cellDim, const Array<int>& cellLID,
00136     CellJacobianBatch& jBatch) const 
00137     {ptr()->getJacobians(cellDim, cellLID, jBatch);}
00138 
00139   /** 
00140    * Compute the diameters of a batch of cells,
00141    * result via reference argument
00142    *
00143    * @param cellDim dimension of the cells whose diameters are to
00144    * be computed
00145    * @param cellLID local indices of the cells for which diameters
00146    * are to be computed
00147    * @param diameters reference to the array of cell diameters
00148    */
00149   virtual void getCellDiameters(int cellDim, const Array<int>& cellLID,
00150     Array<double>& diameters) const 
00151     {ptr()->getCellDiameters(cellDim, cellLID, diameters);}
00152 
00153 
00154   /**
00155    * Map reference quadrature points to physical points on the
00156    * given cells. 
00157    */
00158   void pushForward(int cellDim, const Array<int>& cellLID,
00159     const Array<Point>& refQuadPts,
00160     Array<Point>& physQuadPts) const 
00161     {ptr()->pushForward(cellDim, cellLID, refQuadPts, physQuadPts);}
00162 
00163       
00164 
00165   /** 
00166    * Return the rank of the processor that owns the given cell
00167    */
00168   int ownerProcID(int cellDim, int cellLID) const 
00169     {return ptr()->ownerProcID(cellDim, cellLID);}
00170     
00171 
00172   /** 
00173    * Return the number of facets of the given cell
00174    */
00175   int numFacets(int cellDim, int cellLID, 
00176     int facetDim) const 
00177     {return ptr()->numFacets(cellDim, cellLID, facetDim);}
00178 
00179   /** 
00180    * Return the local ID of a facet cell
00181    * @param cellDim dimension of the cell whose facets are being obtained
00182    * @param cellLID local index of the cell whose
00183    * facets are being obtained
00184    * @param facetDim dimension of the desired facet
00185    * @param facetIndex index into the list of the cell's facets
00186    */
00187   int facetLID(int cellDim, int cellLID,
00188     int facetDim, int facetIndex,
00189     int& facetOrientation) const 
00190     {return ptr()->facetLID(cellDim, cellLID, 
00191         facetDim, facetIndex,
00192         facetOrientation);}
00193 
00194   /**
00195    * Return by reference argument an
00196    *  array containing the LIDs of the facetDim-dimensional 
00197    * facets of the given cell
00198    */
00199   void getFacetArray(int cellDim, int cellLID, int facetDim, 
00200     Array<int>& facetLIDs,
00201     Array<int>& facetOrientations) const 
00202     {ptr()->getFacetArray(cellDim, cellLID, 
00203         facetDim, facetLIDs,
00204         facetOrientations);}
00205 
00206   /** 
00207    * Return a view of an element's zero-dimensional facets
00208    */
00209   const int* elemZeroFacetView(int cellLID) const 
00210     {return ptr()->elemZeroFacetView(cellLID);}
00211 
00212   /** 
00213    * Return by reference argument an array containing
00214    * the LIDs of the facetDim-dimensional facets of the
00215    * given batch of cells 
00216    */
00217   void getFacetLIDs(int cellDim, 
00218     const Array<int>& cellLID,
00219     int facetDim,
00220     Array<int>& facetLID,
00221     Array<int>& facetOrientations) const 
00222     {ptr()->getFacetLIDs(cellDim, cellLID, 
00223         facetDim, facetLID, facetOrientations);}
00224 
00225 
00226   /** 
00227    * Return the number of maximal cofacets of the given cell
00228    */
00229   int numMaxCofacets(int cellDim, int cellLID) const 
00230     {return ptr()->numMaxCofacets(cellDim, cellLID);}
00231 
00232   /** 
00233    * Return the local ID of a maximal cofacet of a cell
00234    * @param cellDim dimension of the cell whose cofacets are being obtained
00235    * @param cellLID local index of the cell whose
00236    * cofacets are being obtained
00237    * @param cofacetIndex [in] index of the maximal cell 
00238    * into the list of the cell's cofacets
00239    * @param facetIndex [out] index of the calling cell
00240    * into the list of the maximal cell's facets
00241    */
00242   int maxCofacetLID(int cellDim, int cellLID,
00243     int cofacetIndex,
00244     int& facetIndex) const 
00245     {return ptr()->maxCofacetLID(cellDim, cellLID, cofacetIndex, facetIndex);}
00246 
00247   /** 
00248    * Get the LIDs of the maximal cofacets for a batch of 
00249    * codimension-one cells. 
00250    *
00251    * \param cellLIDs [in] array of LIDs of the cells whose cofacets are 
00252    * being obtained
00253    * \param cofacets [out] the batch of cofacets
00254    */
00255   void getMaxCofacetLIDs(const Array<int>& cellLIDs,
00256     MaximalCofacetBatch& cofacets) const 
00257     {ptr()->getMaxCofacetLIDs(cellLIDs, cofacets);}
00258 
00259 
00260 
00261   /** 
00262    * Find the cofacets of the given cell
00263    * @param cellDim dimension of the cell whose cofacets are being obtained
00264    * @param cellLID local index of the cell whose
00265    * cofacets are being obtained
00266    * @param cofacetDim dimension of the cofacets to get
00267    * @param cofacetLIDs LIDs of the cofacet
00268    */
00269   void getCofacets(int cellDim, int cellLID,
00270     int cofacetDim, Array<int>& cofacetLIDs) const 
00271     {ptr()->getCofacets(cellDim, cellLID, cofacetDim, cofacetLIDs);}
00272 
00273   /** 
00274    * Find the local ID of a cell given its global index
00275    */
00276   int mapGIDToLID(int cellDim, int globalIndex) const 
00277     {return ptr()->mapGIDToLID(cellDim, globalIndex);}
00278 
00279   /** 
00280    * Determine whether a given cell GID exists on this processor
00281    */
00282   bool hasGID(int cellDim, int globalIndex) const 
00283     {return ptr()->hasGID(cellDim, globalIndex);}
00284 
00285     
00286 
00287   /** 
00288    * Find the global ID of a cell given its local index
00289    */
00290   int mapLIDToGID(int cellDim, int localIndex) const 
00291     {return ptr()->mapLIDToGID(cellDim, localIndex);}
00292 
00293   /**
00294    * Get the type of the given cell
00295    */
00296   CellType cellType(int cellDim) const 
00297     {return ptr()->cellType(cellDim);}
00298 
00299   /** Get the label of the given cell */
00300   int label(int cellDim, int cellLID) const 
00301     {return ptr()->label(cellDim, cellLID);}
00302 
00303   /** Get the labels for a batch of cells */
00304   void getLabels(int cellDim, const Array<int>& cellLID, Array<int>& labels) const 
00305     {ptr()->getLabels(cellDim, cellLID, labels);}
00306 
00307   /** Set the label for the given cell */
00308   void setLabel(int cellDim, int cellLID, int label)
00309     {ptr()->setLabel(cellDim, cellLID, label);}
00310 
00311   /** Get the list of all labels defined for cells of the given dimension */
00312   Set<int> getAllLabelsForDimension(int cellDim) const 
00313     {return ptr()->getAllLabelsForDimension(cellDim);}
00314 
00315   /** 
00316    * Return the number of labels associated with the given dimension.
00317    */
00318   virtual int numLabels(int cellDim) const 
00319     {return getAllLabelsForDimension(cellDim).size();}
00320 
00321   /** 
00322    * Get the cells associated with a specified label. The array 
00323    * cellLID will be filled with those cells of dimension cellDim
00324    * having the given label.
00325    */
00326   void getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const 
00327     {ptr()->getLIDsForLabel(cellDim, label, cellLIDs);}
00328 
00329   /** Get the MPI communicator over which the mesh is distributed */
00330   const MPIComm& comm() const {return ptr()->comm();}
00331 
00332   /** \name Mesh creation methods */
00333   //@{
00334   /** Allocate space for an estimated number of vertices */
00335   void estimateNumVertices(int nPts) 
00336     {creatableMesh()->estimateNumVertices(nPts);}
00337     
00338   /** Allocate space for an estimated number of elements */
00339   void estimateNumElements(int nElems) 
00340     {creatableMesh()->estimateNumElements(nElems);}
00341 
00342     
00343   /** Add a new vertex to the mesh */
00344   int addVertex(int globalIndex, const Point& x,
00345     int ownerProcID, int label)
00346     {return creatableMesh()->addVertex(globalIndex, x, ownerProcID, label);}
00347 
00348   /** Add a new vertex to the mesh */
00349   int addElement(int globalIndex, const Array<int>& vertLID,
00350     int ownerProcID, int label)
00351     {return creatableMesh()->addElement(globalIndex, vertLID, 
00352         ownerProcID, label);}
00353 
00354   /** */
00355   void freezeTopology() 
00356     {creatableMesh()->freezeTopology();}
00357 
00358     
00359   /** \name Reordering */
00360   //@{
00361   /** Set the reordering method to be used with this mesh */
00362   void setReorderer(const CellReorderer& reorderer) 
00363     {ptr()->setReorderer(reorderer);}
00364 
00365   /** Set the reordering method to be used with this mesh */
00366   const CellReordererImplemBase* reorderer() const 
00367     {return ptr()->reorderer();}
00368   //@}
00369     
00370   /** */
00371   static CellReorderer& defaultReorderer()
00372     {
00373       static CellReorderer rtn = new IdentityReorderer();
00374       return rtn;
00375     }
00376 
00377   /** Work out global numberings for the cells of dimension cellDim */
00378   void assignIntermediateCellGIDs(int cellDim) 
00379     {
00380       if (!hasIntermediateGIDs(cellDim))
00381         ptr()->assignIntermediateCellGIDs(cellDim);
00382     }
00383 
00384   /** */
00385   bool hasIntermediateGIDs(int cellDim) const 
00386     {
00387       return ptr()->hasIntermediateGIDs(cellDim);
00388     }
00389 
00390   /** */
00391   void dump(const std::string& filename) const ;
00392 
00393   /** Test the consistency of the mesh numbering scheme 
00394    * across processors. This is meant as a check on Sundance's internal
00395    * logic rather than as a check on the validity of a user's mesh. */
00396   bool checkConsistency(const std::string& filename) const ;
00397 
00398   /** Test the consistency of the mesh numbering scheme 
00399    * across processors. This is meant as a check on Sundance's internal
00400    * logic rather than as a check on the validity of a user's mesh. */
00401   bool checkConsistency(std::ostream& os) const ;
00402 
00403   /** \name Functions for Mesh with hanging nodes */
00404     //@{
00405     /** Function returns true if the mesh allows hanging nodes (by refinement), false otherwise */
00406     bool allowsHangingHodes() const { return ptr()->allowsHangingHodes(); }
00407 
00408     /** Function returns true if the specified element is a "hanging" element
00409      * false otherwise */
00410     bool isElementHangingNode(int cellDim , int cellLID) const
00411         { return ptr()->isElementHangingNode(cellDim , cellLID); }
00412 
00413    /** Returns the index in the parent maxdim Cell of the refinement tree
00414     * @param maxCellLID [in] the LID of the cell */
00415     int indexInParent(int maxCellLID) const
00416         { return ptr()->indexInParent(maxCellLID); }
00417 
00418    /** How many children has a refined element. <br>
00419     * This function provides information of either we have bi or trisection */
00420    int maxChildren() const { return ptr()->maxChildren();}
00421 
00422    /** Function returns the facets of the maxdim parent cell (needed for HN treatment) */
00423    void returnParentFacets( int childCellLID , int dimFacets ,
00424                                  Array<int> &facetsLIDs , int &parentCellLIDs ) const {
00425       ptr()->returnParentFacets( childCellLID , dimFacets , facetsLIDs , parentCellLIDs );
00426    }
00427   //@}
00428 
00429 
00430   /** \name Special Weights Storage for Adaptive Cell Integration */
00431    //@{
00432   /** returns the status of the special weights if they are valid <br>
00433     *  These weights are usually computed for one setting of the curve (Adaptive Cell Integration)*/
00434   bool IsSpecialWeightValid() const {return ptr()->IsSpecialWeightValid();}
00435 
00436   /** specifies if the special weights are valid <br>
00437    *  if this is false then usually the special weights have to be recomputed */
00438   void setSpecialWeightValid(bool val) const { ptr()->setSpecialWeightValid(val);}
00439 
00440   /** deletes all the special weights*/
00441   void flushSpecialWeights() const { ptr()->flushSpecialWeights(); }
00442 
00443   /** verifies if the specified cell with the given dimension has special weights */
00444   bool hasSpecialWeight(int dim, int cellLID) const {return ptr()->hasSpecialWeight( dim, cellLID); }
00445 
00446   /** Sets the special weights */
00447   void setSpecialWeight(int dim, int cellLID, Array<double>& w) const {ptr()->setSpecialWeight(dim, cellLID, w);}
00448 
00449   /** Returns the special weights */
00450   void getSpecialWeight(int dim, int cellLID, Array<double>& w) const {ptr()->getSpecialWeight(dim, cellLID, w);}
00451    //@}
00452 
00453 
00454 
00455   /** \name Store the intersection/quadrature points for the curve/surf integral <br>
00456    *  for a curve or surf integral we need some quadrature points along the curve in one curve <br>
00457    *  These */
00458     //@{
00459 
00460   /** */
00461   bool IsCurvePointsValid() const {return ptr()->IsCurvePointsValid();}
00462 
00463   /**  */
00464   void setCurvePointsValid(bool val)  const {ptr()->setCurvePointsValid(val); }
00465 
00466   /** deletes all the curve points */
00467   void flushCurvePoints() const { ptr()->flushCurvePoints(); }
00468 
00469   /** verifies if the specified maxCell has already precalculated quadrature point for one curve */
00470   bool hasCurvePoints(int maxCellLID , int curveID) const { return ptr()->hasCurvePoints( maxCellLID , curveID); }
00471 
00472   /** Sets the points, curve derivatives and curve normals for one maxCell needed for curve/surf integral*/
00473   void setCurvePoints(int maxCellLID, int curveID , Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const
00474     {ptr()->setCurvePoints( maxCellLID, curveID , points , derivs , normals); }
00475 
00476   /** Gets the points, curve derivatives and curve normals for one maxCell needed for curve/surf integral*/
00477   void getCurvePoints(int maxCellLID,  int curveID , Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const
00478     {ptr()->getCurvePoints( maxCellLID,  curveID ,  points , derivs , normals); }
00479 
00480   //@}
00481 
00482 private:
00483   /** */
00484   IncrementallyCreatableMesh* creatableMesh();
00485 
00486 
00487   /** */
00488   bool checkVertexConsistency(std::ostream& os) const ;
00489   /** */
00490   bool checkCellConsistency(std::ostream& os, int dim) const ;
00491 
00492   /** */
00493   bool checkRemoteEntity(std::ostream& os, int p, int dim, int gid, 
00494     int owner, bool mustExist, int& lid) const ;
00495 
00496   /** */
00497   bool testIdentity(std::ostream& os, int a, int b, const std::string& msg) const ;
00498 
00499   /** */
00500   bool testIdentity(std::ostream& os, 
00501     const Array<int>& a,
00502     const Array<int>& b, const std::string& msg) const ;
00503 };
00504 }
00505 
00506 #endif
SundanceMesh.hpp
