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 #ifdef MOVED_TO_BASE_CLASS 00179 /** 00180 * Get the LIDs of the maximal cofacets for a batch of codimension-one 00181 * cells. 00182 * 00183 * \param cellLIDs [in] array of LIDs of the cells whose cofacets are 00184 * being obtained 00185 * \param cofacets [out] the batch of cofacets 00186 */ 00187 virtual void getMaxCofacetLIDs(const Array<int>& cellLIDs, 00188 MaximalCofacetBatch& cofacets) const ; 00189 #endif 00190 00191 /** 00192 * Find the cofacets of the given cell 00193 * @param cellDim dimension of the cell whose cofacets are being obtained 00194 * @param cellLID local index of the cell whose 00195 * cofacets are being obtained 00196 * @param cofacetDim dimension of the cofacets to get 00197 * @param cofacetLIDs LIDs of the cofacet 00198 */ 00199 void getCofacets(int cellDim, int cellLID, 00200 int cofacetDim, Array<int>& cofacetLIDs) const ; 00201 00202 /** 00203 * Find the local ID of a cell given its global index 00204 */ 00205 virtual int mapGIDToLID(int cellDim, int globalIndex) const ; 00206 00207 /** 00208 * Determine whether a given cell GID exists on this processor 00209 */ 00210 virtual bool hasGID(int cellDim, int globalIndex) const ; 00211 00212 00213 00214 /** 00215 * Find the global ID of a cell given its local index 00216 */ 00217 virtual int mapLIDToGID(int cellDim, int localIndex) const ; 00218 00219 /** 00220 * Get the type of the given cell 00221 */ 00222 virtual CellType cellType(int cellDim) const ; 00223 00224 00225 /** Get the label of the given cell */ 00226 virtual int label(int cellDim, int cellLID) const ; 00227 00228 /** Get the labels for a batch of cells */ 00229 virtual void getLabels(int cellDim, const Array<int>& cellLID, 00230 Array<int>& labels) const ; 00231 00232 00233 00234 /** Get the list of all labels defined for cells of the given dimension */ 00235 virtual Set<int> getAllLabelsForDimension(int cellDim) const ; 00236 00237 /** 00238 * Get the cells associated with a specified label. The array 00239 * cellLID will be filled with those cells of dimension cellDim 00240 * having the given label. 00241 */ 00242 virtual void getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const ; 00243 00244 /** \name Incremental creation methods */ 00245 //@{ 00246 /** 00247 * Add new new vertex to the mesh. 00248 * \param globalIndex the GID of the new vertex 00249 * \param x the spatial position of the new vertex 00250 * \param ownerProcID the processor that "owns" this vertex 00251 * \param label a label for this vertex (optionally used in setting 00252 * loads, boundary conditions, etc) 00253 * \return the LID of the vertex. 00254 */ 00255 virtual int addVertex(int globalIndex, const Point& x, 00256 int ownerProcID, int label); 00257 00258 /** 00259 * Add a new element to the mesh. 00260 * \param globalIndex the GID of the new element 00261 * \param vertexGIDs tuple of GIDs for the vertices defining this element. 00262 * \param ownerProcID the processor that "owns" this element 00263 * \param label a label for this element (optionally used in setting loads, 00264 * material properties, etc) 00265 * \return the LID of the element 00266 */ 00267 virtual int addElement(int globalIndex, const Array<int>& vertexGIDs, 00268 int ownerProcID, int label); 00269 00270 /** Set the label of the given cell */ 00271 virtual void setLabel(int cellDim, int cellLID, int label) 00272 { 00273 labels_[cellDim][cellLID] = label; 00274 } 00275 00276 /** Optional preallocation of space for an estimated number of vertices */ 00277 virtual void estimateNumVertices(int numVertices); 00278 00279 /** Optional preallocation of space for an estimated number of elements */ 00280 virtual void estimateNumElements(int numElements); 00281 00282 /** Coordinate intermediate cell definitions across processors */ 00283 virtual void assignIntermediateCellGIDs(int cellDim) ; 00284 00285 /** */ 00286 virtual bool hasIntermediateGIDs(int dim) const 00287 { 00288 if (dim==1) return hasEdgeGIDs_; 00289 return hasFaceGIDs_; 00290 } 00291 00292 //@} 00293 private: 00294 00295 /** 00296 * Add a new face, first checking to see if it already exists. 00297 * This function is called from within addElement(), 00298 * not by the user, and is therefore private. 00299 * 00300 * \param vertGID The sorted GIDs for the three vertices of the face 00301 * \param vertLID The LIDs for the three vertices of the face 00302 * \param edgeLID{1,2,3} The LIDs for the three edges of the face 00303 * \return LID of this face 00304 */ 00305 int addFace(const Array<int>& vertLID, 00306 const Array<int>& vertGID, 00307 const Array<int>& edgeGID, 00308 int elemLID, 00309 int elemGID); 00310 00311 /** 00312 * Add a new edge, first checking to see if it already exists. 00313 * This function is called from within addElement(), 00314 * not by the user, and is therefore private. 00315 * 00316 * \param vertLID{1,2} 00317 * \param elemLID LID of the element that is adding the edge 00318 * \param elemGID GID of the element that is adding the edge 00319 * \param myFacetNumber facet number of the edge within the element 00320 * \return LID of this edge 00321 */ 00322 int addEdge(int vertLID1, int vertLID2, int elemLID, int elemGID, 00323 int myFacetNumber); 00324 00325 /** 00326 * Check for the presence of the edge (vertLID1, vertLID2) in the mesh. 00327 * \return the edge LID if the edge exists, -1 otherwise 00328 */ 00329 int checkForExistingEdge(int vertLID1, int vertLID2); 00330 00331 /** 00332 * Check whether the face defined by the given vertices exists 00333 * in the mesh. Returns -1 if the face does not exist. 00334 * Called during the synchronization of intermediate cell GIDs. 00335 * \param vertGID{1,2,3} the global indices of the vertices defining the face 00336 * \return the LID of the face 00337 */ 00338 int lookupFace(const Array<int>& vertGID) ; 00339 00340 /** */ 00341 void synchronizeGIDNumbering(int dim, int localCount) ; 00342 00343 /** */ 00344 void resolveEdgeOwnership(int cellDim); 00345 00346 /** */ 00347 std::string cellStr(int dim, const int* verts) const ; 00348 00349 /** */ 00350 std::string cellToStr(int dim, int cellLID) const ; 00351 00352 /** */ 00353 std::string printCells(int dim) const ; 00354 00355 /** */ 00356 void synchronizeNeighborLists(); 00357 00358 00359 00360 /** Number of cells of each dimension. */ 00361 Array<int> numCells_; 00362 00363 /** coordinates of vertices. The index into the array is the vertex LID .*/ 00364 Array<Point> points_; 00365 00366 /** pairs of local vertex indices for the edges, each pair ordered 00367 * from lower to higher <i>global</i> vertex index in order to define 00368 * an absolute edge orientation. Because global vertex indices are used, all 00369 * processors will agree on this orientation, regardless of the orientation 00370 * of the edge as seen by the element first defining it. 00371 * The first index into this 2D array is the edge LID, the second the vertex 00372 * number within the edge. 00373 * */ 00374 ArrayOfTuples<int> edgeVerts_; 00375 00376 /** Tuples of local vertex indices for the faces, with each tuple ordered from 00377 * lowest to highest <i>global</i> index in order to define an absolute edge 00378 * orientation. Because global vertex indices are used, all 00379 * processors will agree on this orientation, regardless of the orientation 00380 * of the face as seen by the element first defining it. 00381 * The first index into this 2D array is the face LID, the second the 00382 * vertex number within the face. */ 00383 ArrayOfTuples<int> faceVertLIDs_; 00384 00385 /** Tuples of global vertex indices for the faces, with each tuple ordered from 00386 * lowest to highest <i>global</i> index in order to define an absolute edge 00387 * orientation. Because global vertex indices are used, all 00388 * processors will agree on this orientation, regardless of the orientation 00389 * of the face as seen by the element first defining it. 00390 * The first index into this 2D array is the face LID, the second the 00391 * vertex number within the face. 00392 * 00393 * Notice that we duplicate the face vertex storage, storing both the 00394 * vertex LIDs and vertex GIDs for each face. This lets us do quick comparison 00395 * with the sorted GID array in order to identify pre-existing faces, while 00396 * also making it possible to retrieve face vertex LID information without 00397 * doing hashtable lookups. 00398 * 00399 */ 00400 ArrayOfTuples<int> faceVertGIDs_; 00401 00402 /** Tuples of local indices for the edges of all faces. The first index 00403 * into this 2D array is the face LID, the second the edge number. */ 00404 ArrayOfTuples<int> faceEdges_; 00405 00406 /** Tuples of edge signs for the faces. The edge sign indicates 00407 * whether the orientation of the edge as given by moving around the face 00408 * is parallel or antiparallel to the absolute orientation of the edge. */ 00409 ArrayOfTuples<int> faceEdgeSigns_; 00410 00411 /** tuples of local vertex indices for the elements. The first index into this 00412 * 2D array is the element LID, the second is the vertex number. */ 00413 ArrayOfTuples<int> elemVerts_; 00414 00415 /** tuples of local edge indices for the elements. The first index into 00416 * this 2D array is the element LID, the second is the edge number. */ 00417 ArrayOfTuples<int> elemEdges_; 00418 00419 /** tuples of edge orientations for the elements, indicating whether 00420 * the orientation of the edge as given by moving around the element 00421 * is parallel or antiparallel to the absolute orientation of the edge. 00422 * The first index into this 2D array is the element LID, the second 00423 * the edge number. 00424 * */ 00425 ArrayOfTuples<int> elemEdgeSigns_; 00426 00427 /** tuples of face LIDs for the elements. The first index is the 00428 * element LID, the second the face number. */ 00429 ArrayOfTuples<int> elemFaces_; 00430 00431 /** tuples of face rotations for the elements, defined relative to the 00432 * absolute orientation of the face. */ 00433 ArrayOfTuples<int> elemFaceRotations_; 00434 00435 /** table for mapping vertex set -> face index */ 00436 Hashtable<VertexView, int> vertexSetToFaceIndexMap_; 00437 00438 /** array of face cofacets for the edges. The first index 00439 * is the edge LID, the second the cofacet number. */ 00440 Array<Array<int> > edgeFaces_; 00441 00442 /** array of element cofacets for the edges. The first index 00443 * is the edge LID, the second the cofacet number. */ 00444 Array<Array<int> > edgeCofacets_; 00445 00446 /** array of element cofacets for the faces. The first index is the 00447 * face LID, the second the cofacet number. */ 00448 Array<Array<int> > faceCofacets_; 00449 00450 /** array of edge cofacets for the vertices. The first index is the 00451 * vertex LID, the second the edge cofacet number. */ 00452 Array<Array<int> > vertEdges_; 00453 00454 /** array of face cofacet LIDs for the vertices. The first index is the 00455 * vertex LID, the second the cofacet number. */ 00456 Array<Array<int> > vertFaces_; 00457 00458 /** array of maximal cofacets for the vertices. The first index is the 00459 * vertex LID, the second the cafacet number. */ 00460 Array<Array<int> > vertCofacets_; 00461 00462 /** array of edge partners for the vertices. The partners are other 00463 * vertices sharing an edge with the specified vertex. */ 00464 Array<Array<int> > vertEdgePartners_; 00465 00466 /** map from local to global cell indices. The first index into this 00467 * 2D array is the cell dimension, the second the cell LID. */ 00468 Array<Array<int> > LIDToGIDMap_; 00469 00470 /** map from global to local cell indices. The array index is the 00471 * cell dimension. The hashtable key is the cell GID, the value the 00472 * cell LID. */ 00473 Array<Hashtable<int, int> > GIDToLIDMap_; 00474 00475 /** Array of labels for the cells */ 00476 Array<Array<int> > labels_; 00477 00478 /** Array of owning processor IDs for the cells. Each cell is owned by 00479 * a single processor that is responsible for assigning global indices, 00480 * DOF numbers, and so on. */ 00481 Array<Array<int> > ownerProcID_; 00482 00483 /** 00484 * Pointer to the pointer at the base of the face vertex GID array. This is 00485 * used to get small-array "views" of the face vertex GID array without 00486 * making copies, resulting in a significant performance improvement 00487 * in the vertex set hashtable lookups to identify pre-existing faces. 00488 * 00489 * We use double rather than single indirection here because as elements 00490 * are added, the face vertex GID array will often be resized, thus changing 00491 * the base pointer. Each vertex set "view" keeps a pointer to the base pointer, 00492 * so that it always remains synchronized with the array of face vertices. 00493 * 00494 * IMPORTANT: any time faceVertGIDs_ is resized, faceVertGIDBase_[0] must 00495 * be reset to the base of the faceVertGIDs_ array so that the vertex 00496 * sets are pointing to the right place. 00497 */ 00498 Array<int*> faceVertGIDBase_; 00499 00500 00501 /** flag indicating whether the edge GIDs have been synchronized */ 00502 bool hasEdgeGIDs_; 00503 00504 /** flag indicating whether the face GIDs have been synchronized */ 00505 bool hasFaceGIDs_; 00506 00507 /** Set of all neighbor processors sharing data with this one */ 00508 Set<int> neighbors_; 00509 00510 /** Whether the neighbor lists have been synchronized across 00511 * processors */ 00512 bool neighborsAreSynchronized_; 00513 00514 00515 }; 00516 00517 } 00518 00519 00520 00521 00522 #endif