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 * SundanceHNMesh2D.h 00032 * 00033 * Created on: April 30, 2010 00034 * Author: benk 00035 */ 00036 00037 #ifndef SUNDANCEHNMESH2D_H_ 00038 #define SUNDANCEHNMESH2D_H_ 00039 00040 #include "SundanceDefs.hpp" 00041 #include "SundanceMeshBase.hpp" 00042 #include "SundanceSet.hpp" 00043 #include "SundancePoint.hpp" 00044 #include "SundanceCellType.hpp" 00045 00046 #include "Teuchos_Array.hpp" 00047 #include "Teuchos_Hashtable.hpp" 00048 00049 #include "SundanceRefinementBase.hpp" 00050 #include "SundanceRefinementClass.hpp" 00051 00052 #include "SundanceDomainDefinition.hpp" 00053 00054 00055 namespace Sundance 00056 { 00057 using namespace Sundance; 00058 using namespace Teuchos; 00059 00060 /** 00061 * Class for 2D hierarchical structured quad Mesh <br> 00062 * The main internal idea of the mesh ist that there are different numberings of the elements 00063 * ID -> each element has an ID (also those which are completely outside the mesh domain, and those 00064 * which are not leaf) <br> 00065 * In the code sometimes internaly LID is used instead of ID due to legacy !! :-P <br> 00066 * GID -> all leaf elements which are visible to Sundance (only the leaf elements should be visible) <br> 00067 * LID -> all the elements which have GID and belong and either belong to the local processor or are in the ghost cells <br>*/ 00068 00069 class HNMesh2D : public MeshBase{ 00070 00071 public: 00072 /** 00073 * The Ctor for the dummy grid with hanging nodes */ 00074 HNMesh2D(int dim, 00075 const MPIComm& comm, 00076 const MeshEntityOrder& order); 00077 00078 /** 00079 * The Ctor for the HNMesh2D grid in 2D*/ 00080 void createMesh( 00081 double position_x, 00082 double position_y, 00083 double offset_x, 00084 double offset_y, 00085 int resolution_x, 00086 int resolution_y, 00087 const RefinementClass& refineClass, 00088 const MeshDomainDef& meshDomain ); 00089 00090 /** Dtor */ 00091 virtual ~HNMesh2D(){;} 00092 00093 /** 00094 * Get the number of cells having dimension dim 00095 */ 00096 virtual int numCells(int dim) const; 00097 00098 /** 00099 * Return the position of the i-th node 00100 */ 00101 virtual Point nodePosition(int i) const; 00102 00103 /** 00104 * Return a view of the i-th node's position 00105 */ 00106 const double* nodePositionView(int i) const; 00107 00108 00109 00110 /** 00111 * Compute the jacobians of a batch of cells, returning the 00112 * result via reference argument 00113 * 00114 * @param cellDim dimension of the cells whose Jacobians are to 00115 * be computed 00116 * @param cellLID local indices of the cells for which Jacobians 00117 * are to be computed 00118 * @param jBatch reference to the resulting Jacobian batch 00119 */ 00120 virtual void getJacobians(int cellDim, const Array<int>& cellLID, 00121 CellJacobianBatch& jBatch) const; 00122 00123 /** 00124 * Compute the diameters of a batch of cells, 00125 * result via reference argument 00126 * 00127 * @param cellDim dimension of the cells whose diameters are to 00128 * be computed 00129 * @param cellLID local indices of the cells for which diameters 00130 * are to be computed 00131 * @param diameters reference to the array of cell diameters 00132 */ 00133 virtual void getCellDiameters(int cellDim, const Array<int>& cellLID, 00134 Array<double>& cellDiameters) const; 00135 00136 00137 /** 00138 * Map reference quadrature points to physical points on the 00139 * given cells. 00140 */ 00141 virtual void pushForward(int cellDim, const Array<int>& cellLID, 00142 const Array<Point>& refQuadPts, 00143 Array<Point>& physQuadPts) const; 00144 00145 /** 00146 * Return the rank of the processor that owns the given cell 00147 */ 00148 virtual int ownerProcID(int cellDim, int cellLID) const; 00149 00150 /** 00151 * Return the number of facets of the given cell 00152 */ 00153 virtual int numFacets(int cellDim, int cellLID, 00154 int facetDim) const; 00155 00156 /** 00157 * Return the local ID of a facet cell 00158 * @param cellDim dimension of the cell whose facets are being obtained 00159 * @param cellLID local index of the cell whose 00160 * facets are being obtained 00161 * @param facetDim dimension of the desired facet 00162 * @param facetIndex index into the list of the cell's facets 00163 * @param facetOrientation orientation of the facet as seen from 00164 * the given cell (returned via reference) 00165 */ 00166 virtual int facetLID(int cellDim, int cellLID, 00167 int facetDim, int facetIndex, 00168 int& facetOrientation) const; 00169 00170 /** 00171 * Return by reference argument an array containing&(elemVerts_.value(cellLID, 0)) 00172 * the LIDs of the facetDim-dimensional facets of the 00173 * given batch of cells 00174 */ 00175 virtual void getFacetLIDs(int cellDim, 00176 const Array<int>& cellLID, 00177 int facetDim, 00178 Array<int>& facetLID, 00179 Array<int>& facetSign) const; 00180 00181 /** 00182 * Return a view of an element's zero-dimensional facets, 00183 * @return an array of integers with the indexes of the points which for it 00184 */ 00185 const int* elemZeroFacetView(int cellLID) const; 00186 00187 /** 00188 * Return the number of maximal cofacets of the given cell 00189 */ 00190 virtual int numMaxCofacets(int cellDim, int cellLID) const; 00191 00192 /** 00193 * Return the local ID of a maximal cofacet cell 00194 * @param cellDim dimension of the cell whose cofacets are being obtained 00195 * @param cellLID local index of the cell whose 00196 * cofacets are being obtained 00197 * @param cofacetIndex which maximal cofacet to get 00198 * @param facetIndex index of the cell cellLID into the list of the 00199 * maximal cell's facets 00200 */ 00201 virtual int maxCofacetLID(int cellDim, int cellLID, 00202 int cofacetIndex, 00203 int& facetIndex) const; 00204 /** 00205 * Get the LIDs of the maximal cofacets for a batch of codimension-one 00206 * cells. 00207 * 00208 * \param cellLIDs [in] array of LIDs of the cells whose cofacets are 00209 * being obtained 00210 * \param cofacets [out] the batch of cofacets 00211 */ 00212 virtual void getMaxCofacetLIDs(const Array<int>& cellLIDs, 00213 MaximalCofacetBatch& cofacets) const; 00214 00215 00216 /** 00217 * Find the cofacets of the given cell 00218 * @param cellDim dimension of the cell whose cofacets are being obtained 00219 * @param cellLID local index of the cell whose 00220 * cofacets are being obtained 00221 * @param cofacetDim dimension of the cofacets to get 00222 * @param cofacetLIDs LIDs of the cofacet 00223 */ 00224 void getCofacets(int cellDim, int cellLID, 00225 int cofacetDim, Array<int>& cofacetLIDs) const; 00226 00227 /** 00228 * Find the local ID of a cell given its global index 00229 */ 00230 virtual int mapGIDToLID(int cellDim, int globalIndex) const; 00231 00232 /** 00233 * Determine whether a given cell GID exists on this processor 00234 */ 00235 virtual bool hasGID(int cellDim, int globalIndex) const; 00236 00237 00238 /** 00239 * Find the global ID of a cell given its local index 00240 */ 00241 virtual int mapLIDToGID(int cellDim, int localIndex) const; 00242 00243 /** 00244 * Get the type of the given cell 00245 */ 00246 virtual CellType cellType(int cellDim) const; 00247 00248 00249 /** Get the label of the given cell */ 00250 virtual int label(int cellDim, int cellLID) const; 00251 00252 /** Get the labels for a batch of cells */ 00253 virtual void getLabels(int cellDim, const Array<int>& cellLID, 00254 Array<int>& labels) const; 00255 00256 00257 00258 /** Get the list of all labels defined for cells of the given dimension */ 00259 virtual Set<int> getAllLabelsForDimension(int cellDim) const; 00260 00261 /** 00262 * Get the cells associated with a specified label. The array 00263 * cellLID will be filled with those cells of dimension cellDim 00264 * having the given label. 00265 */ 00266 virtual void getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const; 00267 00268 00269 /** Set the label of the given cell */ 00270 virtual void setLabel(int cellDim, int cellLID, int label); 00271 00272 /** Coordinate intermediate cell definitions across processors */ 00273 virtual void assignIntermediateCellGIDs(int cellDim); 00274 00275 /** */ 00276 virtual bool hasIntermediateGIDs(int dim) const; 00277 00278 00279 /** Function returns true if the mesh allows hanging nodes (by refinement), 00280 * false otherwise */ 00281 virtual bool allowsHangingHodes() const { return true; } 00282 00283 /** Function returns true if the specified element is a "hanging" element 00284 * false otherwise <br> 00285 * @param dim [in] should be between 0 , D-1 00286 * @param cellLID [in] the local ID of the element */ 00287 virtual bool isElementHangingNode(int cellDim , int cellLID) const ; 00288 00289 /** Returns the index in the parent maxdim Cell of the refinement tree 00290 * @param maxCellLID [in] the LID of the cell */ 00291 virtual int indexInParent(int maxCellLID) const ; 00292 00293 /** How many children has a refined element. <br> 00294 * This function provides information of either we have bi or trisection */ 00295 virtual int maxChildren() const {return 9;} 00296 00297 /** Function returns the facets of the maxdim parent cell (needed for HN treatment) <br> 00298 * @param childCellLID [in] the LID of the maxdim cell, whos parents facets we want 00299 * @param dimFacets [in] the dimension of the facets which we want to have 00300 * @param facetsLIDs [out] the LID of the parents facets (all) in the defined order 00301 * @param parentCellLIDs [out] the maxdim parent cell LID */ 00302 virtual void returnParentFacets( int childCellLID , int dimFacets , 00303 Array<int> &facetsLIDs , int &parentCellLIDs ) const; 00304 00305 private: 00306 00307 /** For HN , returns parent facets, if the facet is not leaf, then return -1 at that place */ 00308 int facetLID_tree(int cellDim, int cellLID, 00309 int facetDim, int facetIndex) const; 00310 00311 /** adds one vertex to the mesh */ 00312 void addVertex(int vertexLID , int ownerProc , bool isHanging , 00313 double coordx , double coordy , const Array<int> &maxCoF); 00314 00315 /** adds one edge to the mesh */ 00316 void addEdge(int edgeLID , int ownerProc , bool isHanging , int edgeVertex , 00317 bool isProcBoundary , bool isMeshBoundary , 00318 const Array<int> &vertexLIDs , const Array<int> &maxCoF); 00319 00320 /** adds one cell(2D) to the mesh <br> 00321 * cell must be always leaf*/ 00322 void addCell(int cellLID , int ownerProc , 00323 int indexInParent , int parentCellLID , int level , 00324 const Array<int> &edgeLIDs , const Array<int> &vertexLIDs); 00325 00326 /** creates the mesh on the coarsest level as it is specified */ 00327 void createCoarseMesh(); 00328 00329 /** Does one single refinement iteration. <br> 00330 * Iterates trough all cells which are owned by the processor and refines if necessary */ 00331 bool oneRefinementIteration(); 00332 00333 /** refine the given cell by cellID, (method assumes that this cell can be refined)*/ 00334 void refineCell(int cellLID); 00335 00336 /** Create Leaf numbering */ 00337 void createLeafNumbering(); 00338 00339 /** Create Leaf numbering, but with a better load balancing for parallel case */ 00340 void createLeafNumbering_sophisticated(); 00341 00342 /** estimate the load of one cell*/ 00343 int estimateCellLoad(int ID); 00344 00345 /** marks the cells recursivly in the tree (and facets) owned by one processor*/ 00346 void markCellsAndFacets(int cellID , int procID); 00347 00348 /** The dimension of the grid*/ 00349 int _dimension; 00350 /** Number of processors */ 00351 int nrProc_; 00352 int myRank_; 00353 /** The communicator */ 00354 const MPIComm& _comm; 00355 /** */ 00356 double _pos_x; 00357 /** */ 00358 double _pos_y; 00359 /** */ 00360 double _ofs_x; 00361 /** */ 00362 double _ofs_y; 00363 /** */ 00364 int _res_x; 00365 /** */ 00366 int _res_y; 00367 /** */ 00368 mutable RefinementClass refineClass_; 00369 /** */ 00370 mutable MeshDomainDef meshDomain_; 00371 00372 //------ Point storage ---- 00373 /** all the global points index is [ID] */ 00374 Array<Point> points_; 00375 /** [3] the nr of ID per dim*/ 00376 Array<int> nrElem_; 00377 /** [3] the nr of owned elements per dim*/ 00378 Array<int> nrElemOwned_; 00379 00380 //----- Facets ----- ; -1 -> MeshBoundary , -2 -> ProcBoundary 00381 00382 /** [cellID][4]*/ 00383 Array< Array<int> > cellsPoints_; 00384 /** [cellID][4]*/ 00385 Array< Array<int> > cellsEdges_; 00386 /** [edgeID][2]*/ 00387 Array< Array<int> > edgePoints_; 00388 /** Information from the edge needs to be stored in one vertex <br> 00389 * We use the traditional mapping , information is put into the 0th Vertex of the Edge <br>*/ 00390 Array<int> edgeVertex_; 00391 00392 // ----- MaxCofacets ----; -1 -> MeshBoundary , -2 -> ProcBoundary 00393 00394 /** [edgeID][2] */ 00395 Array< Array<int> > edgeMaxCoF_; 00396 /** [pointID][4]*/ 00397 Array< Array<int> > pointMaxCoF_; 00398 00399 // ---- Different mesh boundaries ------ 00400 Array<bool> edgeIsProcBonudary_; 00401 Array<bool> edgeIsMeshDomainBonudary_; 00402 00403 //------ Element (processor) ownership ----- 00404 /** contains the ownership of the local elements [dim][ID] */ 00405 Array< Array< short signed int > > elementOwner_; 00406 00407 //---- hierarchical storage ----- 00408 00409 /** [cellID] , the child index in the parent */ 00410 Array<short int> indexInParent_; 00411 /** [cellID] , the LID of the parent cell */ 00412 Array<int> parentCellLID_; 00413 /** [cellID] , actual level of the cell*/ 00414 Array<short int> cellLevel_; 00415 /** [cellID] , if the element is leaf */ 00416 Array<bool> isCellLeaf_; 00417 /** [cellID] , if the cell is complete outside the user defined mesh domain*/ 00418 Array<bool> isCellOut_; 00419 /** [cellID] , children of the cell*/ 00420 Array< Array<int> > cellsChildren_; 00421 // ---- "hanging" info storage --- 00422 /** [pointID] , true if the node is hanging , false otherwise */ 00423 Array<bool> isPointHanging_; 00424 /** [edgeID] , true if the edge is hanging , false otherwise*/ 00425 Array<bool> isEdgeHanging_; 00426 00427 // ---- hanging element and refinement (temporary) storage --- 00428 00429 /** [vertexID] - > { h P1 LID , h P2 LID , h E1 LID , h E2 LID , h E3 LID } */ 00430 Array< Hashtable< int, Array<int> > > hangElmStore_; 00431 /** Neighbor Cell can mark the cell to provoke refinement */ 00432 Array<short int> refineCell_; 00433 00434 // ---- leaf mapping , GID and LID --- (points do not need this, all points are also leaf points) 00435 00436 /** [leaf_vertexLID] , the value must be a positive number */ 00437 Array<int> vertexLeafToLIDMapping_; 00438 /** [leaf_edgeLID] , the value must be a positive number */ 00439 Array<int> edgeLeafToLIDMapping_; 00440 /** [leaf_cellLID] , the value must be a positive number */ 00441 Array<int> cellLeafToLIDMapping_; 00442 /** [vertexLID] if vertex is inside the domain then > 0 , -1 otherwise */ 00443 Array<int> vertexLIDToLeafMapping_; 00444 /** [edgeLID] if edge is leaf(or inside the domain) then > 0 , -1 otherwise */ 00445 Array<int> edgeLIDToLeafMapping_; 00446 /** [cellLID] if cell is leaf(or inside the domain) then > 0 , -1 otherwise */ 00447 Array<int> cellLIDToLeafMapping_; 00448 00449 /** leaf LID numbering*/ 00450 int nrVertexLeafLID_; 00451 int nrCellLeafLID_; 00452 int nrEdgeLeafLID_; 00453 00454 /** [leaf_vertexGID] , the value must be a positive number */ 00455 Array<int> vertexLeafToGIDMapping_; 00456 /** [leaf_edgeGID] , the value must be a positive number */ 00457 Array<int> edgeLeafToGIDMapping_; 00458 /** [leaf_cellGID] , the value must be a positive number */ 00459 Array<int> cellLeafToGIDMapping_; 00460 /** [vertexGID] if vertex is inside the domain then > 0 , -1 otherwise */ 00461 Array<int> vertexGIDToLeafMapping_; 00462 /** [edgeGID] if edge is leaf(or inside the domain) then > 0 , -1 otherwise */ 00463 Array<int> edgeGIDToLeafMapping_; 00464 /** [cellGID] if cell is leaf(or inside the domain) then > 0 , -1 otherwise */ 00465 Array<int> cellGIDToLeafMapping_; 00466 00467 /** leaf GID numbering*/ 00468 int nrVertexLeafGID_; 00469 int nrCellLeafGID_; 00470 int nrEdgeLeafGID_; 00471 00472 // ------------- static data ------------------- 00473 00474 /** the offset in the X coordinate on the reference cell*/ 00475 static int offs_Points_x_[4]; 00476 00477 /** the offset in the Y coordinate on the reference cell*/ 00478 static int offs_Points_y_[4]; 00479 00480 /** stores the facet information on the reference Cell*/ 00481 static int edge_Points_localIndex[4][2]; 00482 00483 }; 00484 } 00485 00486 00487 #endif /* SUNDANCEHNMESH2D_H_ */