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_BASICSIMPLICIALMESH_H 00032 #define SUNDANCE_BASICSIMPLICIALMESH_H 00033 00034 00035 00036 #include "SundanceDefs.hpp" 00037 #include "SundanceMeshBase.hpp" 00038 #include "SundanceSet.hpp" 00039 #include "SundanceBasicVertexView.hpp" 00040 #include "SundanceArrayOfTuples.hpp" 00041 #include "SundanceIncrementallyCreatableMesh.hpp" 00042 #include "Teuchos_Array.hpp" 00043 #include "Teuchos_Hashtable.hpp" 00044 00045 namespace Sundance 00046 { 00047 00048 /** 00049 * A no-frills parallel simplicial mesh 00050 */ 00051 class BasicSimplicialMesh : public IncrementallyCreatableMesh 00052 { 00053 public: 00054 /** */ 00055 BasicSimplicialMesh(int dim, const MPIComm& comm, 00056 const MeshEntityOrder& order); 00057 00058 /** */ 00059 virtual ~BasicSimplicialMesh(){;} 00060 00061 /** 00062 * Get the number of cells having dimension dim 00063 */ 00064 virtual int numCells(int dim) const ; 00065 00066 /** 00067 * Return the position of the i-th node 00068 */ 00069 virtual Point nodePosition(int i) const {return points_[i];} 00070 00071 /** 00072 * Return a view of the i-th node's position 00073 */ 00074 const double* nodePositionView(int i) const {return &(points_[i][0]);} 00075 00076 00077 00078 /** 00079 * Compute the jacobians of a batch of cells, returning the 00080 * result via reference argument 00081 * 00082 * @param cellDim dimension of the cells whose Jacobians are to 00083 * be computed 00084 * @param cellLID local indices of the cells for which Jacobians 00085 * are to be computed 00086 * @param jBatch reference to the resulting Jacobian batch 00087 */ 00088 virtual void getJacobians(int cellDim, const Array<int>& cellLID, 00089 CellJacobianBatch& jBatch) const ; 00090 00091 /** 00092 * Compute the diameters of a batch of cells, 00093 * result via reference argument 00094 * 00095 * @param cellDim dimension of the cells whose diameters are to 00096 * be computed 00097 * @param cellLID local indices of the cells for which diameters 00098 * are to be computed 00099 * @param diameters reference to the array of cell diameters 00100 */ 00101 virtual void getCellDiameters(int cellDim, const Array<int>& cellLID, 00102 Array<double>& diameters) const ; 00103 00104 00105 /** 00106 * Map reference quadrature points to physical points on the 00107 * given cells. 00108 */ 00109 virtual void pushForward(int cellDim, const Array<int>& cellLID, 00110 const Array<Point>& refQuadPts, 00111 Array<Point>& physQuadPts) const ; 00112 00113 00114 00115 /** 00116 * Return the rank of the processor that owns the given cell 00117 */ 00118 virtual int ownerProcID(int cellDim, int cellLID) const ; 00119 00120 00121 00122 /** 00123 * Return the number of facets of the given cell 00124 */ 00125 virtual int numFacets(int cellDim, int cellLID, 00126 int facetDim) const ; 00127 00128 /** 00129 * Return the local ID of a facet cell 00130 * @param cellDim dimension of the cell whose facets are being obtained 00131 * @param cellLID local index of the cell whose 00132 * facets are being obtained 00133 * @param facetDim dimension of the desired facet 00134 * @param facetIndex index into the list of the cell's facets 00135 * @param facetOrientation orientation of the facet as seen from 00136 * the given cell (returned via reference) 00137 */ 00138 virtual int facetLID(int cellDim, int cellLID, 00139 int facetDim, int facetIndex, 00140 int& facetOrientation) const ; 00141 00142 00143 /** 00144 * Return by reference argument an array containing 00145 * the LIDs of the facetDim-dimensional facets of the 00146 * given batch of cells 00147 */ 00148 virtual void getFacetLIDs(int cellDim, 00149 const Array<int>& cellLID, 00150 int facetDim, 00151 Array<int>& facetLID, 00152 Array<int>& facetOrientations) const ; 00153 00154 /** 00155 * Return a view of an element's zero-dimensional facets 00156 */ 00157 const int* elemZeroFacetView(int cellLID) const 00158 {return &(elemVerts_.value(cellLID, 0));} 00159 00160 /** 00161 * Return the number of maximal cofacets of the given cell 00162 */ 00163 virtual int numMaxCofacets(int cellDim, int cellLID) const ; 00164 00165 /** 00166 * Return the local ID of a maximal cofacet cell 00167 * @param cellDim dimension of the cell whose cofacets are being obtained 00168 * @param cellLID local index of the cell whose 00169 * cofacets are being obtained 00170 * @param cofacetIndex which maximal cofacet to get 00171 * @param cofacetIndex index of the cell cellLID into the list of the 00172 * maximal cell's facets 00173 */ 00174 virtual int maxCofacetLID(int cellDim, int cellLID, 00175 int cofacetIndex, 00176 int& facetIndex) const ; 00177 /** 00178 * Get the LIDs of the maximal cofacets for a batch of codimension-one 00179 * cells. 00180 * 00181 * \param cellLIDs [in] array of LIDs of the cells whose cofacets are 00182 * being obtained 00183 * \param cofacets [out] the batch of cofacets 00184 */ 00185 virtual void getMaxCofacetLIDs(const Array<int>& cellLIDs, 00186 MaximalCofacetBatch& cofacets) const ; 00187 00188 00189 /** 00190 * Find the cofacets of the given cell 00191 * @param cellDim dimension of the cell whose cofacets are being obtained 00192 * @param cellLID local index of the cell whose 00193 * cofacets are being obtained 00194 * @param cofacetDim dimension of the cofacets to get 00195 * @param cofacetLIDs LIDs of the cofacet 00196 */ 00197 void getCofacets(int cellDim, int cellLID, 00198 int cofacetDim, Array<int>& cofacetLIDs) const ; 00199 00200 /** 00201 * Find the local ID of a cell given its global index 00202 */ 00203 virtual int mapGIDToLID(int cellDim, int globalIndex) const ; 00204 00205 /** 00206 * Determine whether a given cell GID exists on this processor 00207 */ 00208 virtual bool hasGID(int cellDim, int globalIndex) const ; 00209 00210 00211 00212 /** 00213 * Find the global ID of a cell given its local index 00214 */ 00215 virtual int mapLIDToGID(int cellDim, int localIndex) const ; 00216 00217 /** 00218 * Get the type of the given cell 00219 */ 00220 virtual CellType cellType(int cellDim) const ; 00221 00222 00223 /** Get the label of the given cell */ 00224 virtual int label(int cellDim, int cellLID) const ; 00225 00226 /** Get the labels for a batch of cells */ 00227 virtual void getLabels(int cellDim, const Array<int>& cellLID, 00228 Array<int>& labels) const ; 00229 00230 00231 00232 /** Get the list of all labels defined for cells of the given dimension */ 00233 virtual Set<int> getAllLabelsForDimension(int cellDim) const ; 00234 00235 /** 00236 * Get the cells associated with a specified label. The array 00237 * cellLID will be filled with those cells of dimension cellDim 00238 * having the given label. 00239 */ 00240 virtual void getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const ; 00241 00242 /** \name Incremental creation methods */ 00243 //@{ 00244 /** 00245 * Add new new vertex to the mesh. 00246 * \param globalIndex the GID of the new vertex 00247 * \param x the spatial position of the new vertex 00248 * \param ownerProcID the processor that "owns" this vertex 00249 * \param label a label for this vertex (optionally used in setting 00250 * loads, boundary conditions, etc) 00251 * \return the LID of the vertex. 00252 */ 00253 virtual int addVertex(int globalIndex, const Point& x, 00254 int ownerProcID, int label); 00255 00256 /** 00257 * Add a new element to the mesh. 00258 * \param globalIndex the GID of the new element 00259 * \param vertexGIDs tuple of GIDs for the vertices defining this element. 00260 * \param ownerProcID the processor that "owns" this element 00261 * \param label a label for this element (optionally used in setting loads, 00262 * material properties, etc) 00263 * \return the LID of the element 00264 */ 00265 virtual int addElement(int globalIndex, const Array<int>& vertexGIDs, 00266 int ownerProcID, int label); 00267 00268 /** Set the label of the given cell */ 00269 virtual void setLabel(int cellDim, int cellLID, int label) 00270 { 00271 labels_[cellDim][cellLID] = label; 00272 } 00273 00274 /** Optional preallocation of space for an estimated number of vertices */ 00275 virtual void estimateNumVertices(int numVertices); 00276 00277 /** Optional preallocation of space for an estimated number of elements */ 00278 virtual void estimateNumElements(int numElements); 00279 00280 /** Coordinate intermediate cell definitions across processors */ 00281 virtual void assignIntermediateCellGIDs(int cellDim) ; 00282 00283 /** */ 00284 virtual bool hasIntermediateGIDs(int dim) const 00285 { 00286 if (dim==1) return hasEdgeGIDs_; 00287 return hasFaceGIDs_; 00288 } 00289 00290 //@} 00291 private: 00292 00293 /** 00294 * Add a new face, first checking to see if it already exists. 00295 * This function is called from within addElement(), 00296 * not by the user, and is therefore private. 00297 * 00298 * \param vertGID The sorted GIDs for the three vertices of the face 00299 * \param vertLID The LIDs for the three vertices of the face 00300 * \param edgeLID{1,2,3} The LIDs for the three edges of the face 00301 * \return LID of this face 00302 */ 00303 int addFace(const Array<int>& vertLID, 00304 const Array<int>& vertGID, 00305 const Array<int>& edgeGID, 00306 int elemLID, 00307 int elemGID); 00308 00309 /** 00310 * Add a new edge, first checking to see if it already exists. 00311 * This function is called from within addElement(), 00312 * not by the user, and is therefore private. 00313 * 00314 * \param vertLID{1,2} 00315 * \param elemLID LID of the element that is adding the edge 00316 * \param elemGID GID of the element that is adding the edge 00317 * \param myFacetNumber facet number of the edge within the element 00318 * \return LID of this edge 00319 */ 00320 int addEdge(int vertLID1, int vertLID2, int elemLID, int elemGID, 00321 int myFacetNumber); 00322 00323 /** 00324 * Check for the presence of the edge (vertLID1, vertLID2) in the mesh. 00325 * \return the edge LID if the edge exists, -1 otherwise 00326 */ 00327 int checkForExistingEdge(int vertLID1, int vertLID2); 00328 00329 /** 00330 * Check whether the face defined by the given vertices exists 00331 * in the mesh. Returns -1 if the face does not exist. 00332 * Called during the synchronization of intermediate cell GIDs. 00333 * \param vertGID{1,2,3} the global indices of the vertices defining the face 00334 * \return the LID of the face 00335 */ 00336 int lookupFace(const Array<int>& vertGID) ; 00337 00338 /** */ 00339 void synchronizeGIDNumbering(int dim, int localCount) ; 00340 00341 /** */ 00342 void resolveEdgeOwnership(int cellDim); 00343 00344 /** */ 00345 std::string cellStr(int dim, const int* verts) const ; 00346 00347 /** */ 00348 std::string cellToStr(int dim, int cellLID) const ; 00349 00350 /** */ 00351 std::string printCells(int dim) const ; 00352 00353 /** */ 00354 void synchronizeNeighborLists(); 00355 00356 00357 00358 /** Number of cells of each dimension. */ 00359 Array<int> numCells_; 00360 00361 /** coordinates of vertices. The index into the array is the vertex LID .*/ 00362 Array<Point> points_; 00363 00364 /** pairs of local vertex indices for the edges, each pair ordered 00365 * from lower to higher <i>global</i> vertex index in order to define 00366 * an absolute edge orientation. Because global vertex indices are used, all 00367 * processors will agree on this orientation, regardless of the orientation 00368 * of the edge as seen by the element first defining it. 00369 * The first index into this 2D array is the edge LID, the second the vertex 00370 * number within the edge. 00371 * */ 00372 ArrayOfTuples<int> edgeVerts_; 00373 00374 /** Tuples of local vertex indices for the faces, with each tuple ordered from 00375 * lowest to highest <i>global</i> index in order to define an absolute edge 00376 * orientation. Because global vertex indices are used, all 00377 * processors will agree on this orientation, regardless of the orientation 00378 * of the face as seen by the element first defining it. 00379 * The first index into this 2D array is the face LID, the second the 00380 * vertex number within the face. */ 00381 ArrayOfTuples<int> faceVertLIDs_; 00382 00383 /** Tuples of global vertex indices for the faces, with each tuple ordered from 00384 * lowest to highest <i>global</i> index in order to define an absolute edge 00385 * orientation. Because global vertex indices are used, all 00386 * processors will agree on this orientation, regardless of the orientation 00387 * of the face as seen by the element first defining it. 00388 * The first index into this 2D array is the face LID, the second the 00389 * vertex number within the face. 00390 * 00391 * Notice that we duplicate the face vertex storage, storing both the 00392 * vertex LIDs and vertex GIDs for each face. This lets us do quick comparison 00393 * with the sorted GID array in order to identify pre-existing faces, while 00394 * also making it possible to retrieve face vertex LID information without 00395 * doing hashtable lookups. 00396 * 00397 */ 00398 ArrayOfTuples<int> faceVertGIDs_; 00399 00400 /** Tuples of local indices for the edges of all faces. The first index 00401 * into this 2D array is the face LID, the second the edge number. */ 00402 ArrayOfTuples<int> faceEdges_; 00403 00404 /** Tuples of edge signs for the faces. The edge sign indicates 00405 * whether the orientation of the edge as given by moving around the face 00406 * is parallel or antiparallel to the absolute orientation of the edge. */ 00407 ArrayOfTuples<int> faceEdgeSigns_; 00408 00409 /** tuples of local vertex indices for the elements. The first index into this 00410 * 2D array is the element LID, the second is the vertex number. */ 00411 ArrayOfTuples<int> elemVerts_; 00412 00413 /** tuples of local edge indices for the elements. The first index into 00414 * this 2D array is the element LID, the second is the edge number. */ 00415 ArrayOfTuples<int> elemEdges_; 00416 00417 /** tuples of edge orientations for the elements, indicating whether 00418 * the orientation of the edge as given by moving around the element 00419 * is parallel or antiparallel to the absolute orientation of the edge. 00420 * The first index into this 2D array is the element LID, the second 00421 * the edge number. 00422 * */ 00423 ArrayOfTuples<int> elemEdgeSigns_; 00424 00425 /** tuples of face LIDs for the elements. The first index is the 00426 * element LID, the second the face number. */ 00427 ArrayOfTuples<int> elemFaces_; 00428 00429 /** tuples of face rotations for the elements, defined relative to the 00430 * absolute orientation of the face. */ 00431 ArrayOfTuples<int> elemFaceRotations_; 00432 00433 /** table for mapping vertex set -> face index */ 00434 Hashtable<VertexView, int> vertexSetToFaceIndexMap_; 00435 00436 /** array of face cofacets for the edges. The first index 00437 * is the edge LID, the second the cofacet number. */ 00438 Array<Array<int> > edgeFaces_; 00439 00440 /** array of element cofacets for the edges. The first index 00441 * is the edge LID, the second the cofacet number. */ 00442 Array<Array<int> > edgeCofacets_; 00443 00444 /** array of element cofacets for the faces. The first index is the 00445 * face LID, the second the cofacet number. */ 00446 Array<Array<int> > faceCofacets_; 00447 00448 /** array of edge cofacets for the vertices. The first index is the 00449 * vertex LID, the second the edge cofacet number. */ 00450 Array<Array<int> > vertEdges_; 00451 00452 /** array of face cofacet LIDs for the vertices. The first index is the 00453 * vertex LID, the second the cofacet number. */ 00454 Array<Array<int> > vertFaces_; 00455 00456 /** array of maximal cofacets for the vertices. The first index is the 00457 * vertex LID, the second the cafacet number. */ 00458 Array<Array<int> > vertCofacets_; 00459 00460 /** array of edge partners for the vertices. The partners are other 00461 * vertices sharing an edge with the specified vertex. */ 00462 Array<Array<int> > vertEdgePartners_; 00463 00464 /** map from local to global cell indices. The first index into this 00465 * 2D array is the cell dimension, the second the cell LID. */ 00466 Array<Array<int> > LIDToGIDMap_; 00467 00468 /** map from global to local cell indices. The array index is the 00469 * cell dimension. The hashtable key is the cell GID, the value the 00470 * cell LID. */ 00471 Array<Hashtable<int, int> > GIDToLIDMap_; 00472 00473 /** Array of labels for the cells */ 00474 Array<Array<int> > labels_; 00475 00476 /** Array of owning processor IDs for the cells. Each cell is owned by 00477 * a single processor that is responsible for assigning global indices, 00478 * DOF numbers, and so on. */ 00479 Array<Array<int> > ownerProcID_; 00480 00481 /** 00482 * Pointer to the pointer at the base of the face vertex GID array. This is 00483 * used to get small-array "views" of the face vertex GID array without 00484 * making copies, resulting in a significant performance improvement 00485 * in the vertex set hashtable lookups to identify pre-existing faces. 00486 * 00487 * We use double rather than single indirection here because as elements 00488 * are added, the face vertex GID array will often be resized, thus changing 00489 * the base pointer. Each vertex set "view" keeps a pointer to the base pointer, 00490 * so that it always remains synchronized with the array of face vertices. 00491 * 00492 * IMPORTANT: any time faceVertGIDs_ is resized, faceVertGIDBase_[0] must 00493 * be reset to the base of the faceVertGIDs_ array so that the vertex 00494 * sets are pointing to the right place. 00495 */ 00496 Array<int*> faceVertGIDBase_; 00497 00498 00499 /** flag indicating whether the edge GIDs have been synchronized */ 00500 bool hasEdgeGIDs_; 00501 00502 /** flag indicating whether the face GIDs have been synchronized */ 00503 bool hasFaceGIDs_; 00504 00505 /** Set of all neighbor processors sharing data with this one */ 00506 Set<int> neighbors_; 00507 00508 /** Whether the neighbor lists have been synchronized across 00509 * processors */ 00510 bool neighborsAreSynchronized_; 00511 00512 00513 }; 00514 00515 } 00516 00517 00518 00519 00520 #endif