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 * SundanceHNMesh3D.h 00032 * 00033 * Created on: May 30, 2010 00034 * Author: benk 00035 */ 00036 00037 #ifndef SUNDANCEHNMESH3D_H_ 00038 #define SUNDANCEHNMESH3D_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 3D 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 HNMesh3D : public MeshBase{ 00070 00071 public: 00072 /** 00073 * The Ctor for the dummy grid with hanging nodes */ 00074 HNMesh3D(int dim, 00075 const MPIComm& comm, 00076 const MeshEntityOrder& order); 00077 00078 /** 00079 * The Ctor for the HNMesh3D grid in 3D*/ 00080 void createMesh( 00081 double position_x, 00082 double position_y, 00083 double position_z, 00084 double offset_x, 00085 double offset_y, 00086 double offset_z, 00087 int resolution_x, 00088 int resolution_y, 00089 int resolution_z, 00090 const RefinementClass& refineClass, 00091 const MeshDomainDef& meshDomain ); 00092 00093 /** Dtor */ 00094 virtual ~HNMesh3D(){;} 00095 00096 /** 00097 * Get the number of cells having dimension dim 00098 */ 00099 virtual int numCells(int dim) const; 00100 00101 /** 00102 * Return the position of the i-th node 00103 */ 00104 virtual Point nodePosition(int i) const; 00105 00106 /** 00107 * Return a view of the i-th node's position 00108 */ 00109 const double* nodePositionView(int i) const; 00110 00111 00112 00113 /** 00114 * Compute the jacobians of a batch of cells, returning the 00115 * result via reference argument 00116 * 00117 * @param cellDim dimension of the cells whose Jacobians are to 00118 * be computed 00119 * @param cellLID local indices of the cells for which Jacobians 00120 * are to be computed 00121 * @param jBatch reference to the resulting Jacobian batch 00122 */ 00123 virtual void getJacobians(int cellDim, const Array<int>& cellLID, 00124 CellJacobianBatch& jBatch) const; 00125 00126 /** 00127 * Compute the diameters of a batch of cells, 00128 * result via reference argument 00129 * 00130 * @param cellDim dimension of the cells whose diameters are to 00131 * be computed 00132 * @param cellLID local indices of the cells for which diameters 00133 * are to be computed 00134 * @param diameters reference to the array of cell diameters 00135 */ 00136 virtual void getCellDiameters(int cellDim, const Array<int>& cellLID, 00137 Array<double>& cellDiameters) const; 00138 00139 00140 /** 00141 * Map reference quadrature points to physical points on the 00142 * given cells. 00143 */ 00144 virtual void pushForward(int cellDim, const Array<int>& cellLID, 00145 const Array<Point>& refQuadPts, 00146 Array<Point>& physQuadPts) const; 00147 00148 /** 00149 * Return the rank of the processor that owns the given cell 00150 */ 00151 virtual int ownerProcID(int cellDim, int cellLID) const; 00152 00153 /** 00154 * Return the number of facets of the given cell 00155 */ 00156 virtual int numFacets(int cellDim, int cellLID, 00157 int facetDim) const; 00158 00159 /** 00160 * Return the local ID of a facet cell 00161 * @param cellDim dimension of the cell whose facets are being obtained 00162 * @param cellLID local index of the cell whose 00163 * facets are being obtained 00164 * @param facetDim dimension of the desired facet 00165 * @param facetIndex index into the list of the cell's facets 00166 * @param facetOrientation orientation of the facet as seen from 00167 * the given cell (returned via reference) 00168 */ 00169 virtual int facetLID(int cellDim, int cellLID, 00170 int facetDim, int facetIndex, 00171 int& facetOrientation) const; 00172 00173 /** 00174 * Return by reference argument an array containing&(elemVerts_.value(cellLID, 0)) 00175 * the LIDs of the facetDim-dimensional facets of the 00176 * given batch of cells 00177 */ 00178 virtual void getFacetLIDs(int cellDim, 00179 const Array<int>& cellLID, 00180 int facetDim, 00181 Array<int>& facetLID, 00182 Array<int>& facetSign) const; 00183 00184 /** 00185 * Return a view of an element's zero-dimensional facets, 00186 * @return an array of integers with the indexes of the points which for it 00187 */ 00188 const int* elemZeroFacetView(int cellLID) const; 00189 00190 /** 00191 * Return the number of maximal cofacets of the given cell 00192 */ 00193 virtual int numMaxCofacets(int cellDim, int cellLID) const; 00194 00195 /** 00196 * Return the local ID of a maximal cofacet cell 00197 * @param cellDim dimension of the cell whose cofacets are being obtained 00198 * @param cellLID local index of the cell whose 00199 * cofacets are being obtained 00200 * @param cofacetIndex which maximal cofacet to get 00201 * @param facetIndex index of the cell cellLID into the list of the 00202 * maximal cell's facets 00203 */ 00204 virtual int maxCofacetLID(int cellDim, int cellLID, 00205 int cofacetIndex, 00206 int& facetIndex) const; 00207 /** 00208 * Get the LIDs of the maximal cofacets for a batch of codimension-one 00209 * cells. 00210 * 00211 * \param cellLIDs [in] array of LIDs of the cells whose cofacets are 00212 * being obtained 00213 * \param cofacets [out] the batch of cofacets 00214 */ 00215 virtual void getMaxCofacetLIDs(const Array<int>& cellLIDs, 00216 MaximalCofacetBatch& cofacets) const; 00217 00218 00219 /** 00220 * Find the cofacets of the given cell 00221 * @param cellDim dimension of the cell whose cofacets are being obtained 00222 * @param cellLID local index of the cell whose 00223 * cofacets are being obtained 00224 * @param cofacetDim dimension of the cofacets to get 00225 * @param cofacetLIDs LIDs of the cofacet 00226 */ 00227 void getCofacets(int cellDim, int cellLID, 00228 int cofacetDim, Array<int>& cofacetLIDs) const; 00229 00230 /** 00231 * Find the local ID of a cell given its global index 00232 */ 00233 virtual int mapGIDToLID(int cellDim, int globalIndex) const; 00234 00235 /** 00236 * Determine whether a given cell GID exists on this processor 00237 */ 00238 virtual bool hasGID(int cellDim, int globalIndex) const; 00239 00240 00241 /** 00242 * Find the global ID of a cell given its local index 00243 */ 00244 virtual int mapLIDToGID(int cellDim, int localIndex) const; 00245 00246 /** 00247 * Get the type of the given cell 00248 */ 00249 virtual CellType cellType(int cellDim) const; 00250 00251 00252 /** Get the label of the given cell */ 00253 virtual int label(int cellDim, int cellLID) const; 00254 00255 /** Get the labels for a batch of cells */ 00256 virtual void getLabels(int cellDim, const Array<int>& cellLID, 00257 Array<int>& labels) const; 00258 00259 00260 00261 /** Get the list of all labels defined for cells of the given dimension */ 00262 virtual Set<int> getAllLabelsForDimension(int cellDim) const; 00263 00264 /** 00265 * Get the cells associated with a specified label. The array 00266 * cellLID will be filled with those cells of dimension cellDim 00267 * having the given label. 00268 */ 00269 virtual void getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const; 00270 00271 00272 /** Set the label of the given cell */ 00273 virtual void setLabel(int cellDim, int cellLID, int label); 00274 00275 /** Coordinate intermediate cell definitions across processors */ 00276 virtual void assignIntermediateCellGIDs(int cellDim); 00277 00278 /** */ 00279 virtual bool hasIntermediateGIDs(int dim) const; 00280 00281 00282 /** Function returns true if the mesh allows hanging nodes (by refinement), 00283 * false otherwise */ 00284 virtual bool allowsHangingHodes() const { return true; } 00285 00286 /** Function returns true if the specified element is a "hanging" element 00287 * false otherwise <br> 00288 * @param dim [in] should be between 0 , D-1 00289 * @param cellLID [in] the local ID of the element */ 00290 virtual bool isElementHangingNode(int cellDim , int cellLID) const ; 00291 00292 /** Returns the index in the parent maxdim Cell of the refinement tree 00293 * @param maxCellLID [in] the LID of the cell */ 00294 virtual int indexInParent(int maxCellLID) const ; 00295 00296 /** How many children has a refined element. <br> 00297 * This function provides information of either we have bi or trisection */ 00298 virtual int maxChildren() const {return 27;} 00299 00300 /** Function returns the facets of the maxdim parent cell (needed for HN treatment) <br> 00301 * @param childCellLID [in] the LID of the maxdim cell, whos parents facets we want 00302 * @param dimFacets [in] the dimension of the facets which we want to have 00303 * @param facetsLIDs [out] the LID of the parents facets (all) in the defined order 00304 * @param parentCellLIDs [out] the maxdim parent cell LID */ 00305 virtual void returnParentFacets( int childCellLID , int dimFacets , 00306 Array<int> &facetsLIDs , int &parentCellLIDs ) const; 00307 00308 private: 00309 00310 /** For HN , returns parent facets, if the facet is not leaf, then return -1 at that place */ 00311 int facetLID_tree(int cellDim, int cellLID, 00312 int facetDim, int facetIndex) const; 00313 00314 /** adds one vertex to the mesh */ 00315 void addVertex(int vertexLID , int ownerProc , bool isHanging , 00316 double coordx , double coordy , double coordz , const Array<int> &maxCoF); 00317 00318 /** adds one edge to the mesh */ 00319 void addEdge(int edgeLID , int ownerProc , bool isHanging , int edgeOrientation , 00320 const Array<int> &vertexLIDs , const Array<int> &maxCoF); 00321 00322 /** adds one edge to the mesh */ 00323 void addFace(int faceLID , int ownerProc , bool isHanging , int faceOrientation , 00324 const Array<int> &vertexLIDs , const Array<int> &edgeLIDs , 00325 const Array<int> &maxCoF); 00326 00327 /** adds one cell(3D) to the mesh <br> 00328 * cell must be always leaf*/ 00329 void addCell(int cellLID , int ownerProc , 00330 int indexInParent , int parentCellLID , int level , 00331 const Array<int> &faceLIDs , const Array<int> &edgeLIDs , 00332 const Array<int> &vertexLIDs) ; 00333 00334 /** creates the mesh on the coarsest level as it is specified */ 00335 void createCoarseMesh(); 00336 00337 /** Does one single refinement iteration. <br> 00338 * Iterates trough all cells which are owned by the processor and refines if necessary */ 00339 bool oneRefinementIteration(); 00340 00341 /** refine the given cell by cellID, (method assumes that this cell can be refined)*/ 00342 void refineCell(int cellLID); 00343 00344 /** Create Leaf numbering */ 00345 void createLeafNumbering(); 00346 00347 /** Create Leaf numbering, but with a better load balancing for parallel case */ 00348 void createLeafNumbering_sophisticated(); 00349 00350 /** estimate the load of one cell*/ 00351 int estimateCellLoad(int ID); 00352 00353 /** marks the cells recursivly in the tree (and facets) owned by one processor*/ 00354 void markCellsAndFacets(int cellID , int procID); 00355 00356 /** this updates the array with the local index of vertex, edge and face in the static array <br> 00357 * this method is only used in the coarse mesh creation */ 00358 void updateLocalCoarseNumbering(int ix , int iy , int iz , int Nx , int Ny); 00359 00360 /** Used for refinemement 00361 * @param cellDim the requested elements dimension (if exists) 00362 * @param usedge , whould we look for that element in edge or in face 00363 * @param parentID , ID of the parent edge or face 00364 * @param parentOffset , the offset in the parent */ 00365 int getHangingElement(int cellDim, bool useEdge, int parentID , int parentOffset); 00366 00367 /** Used for refinemement 00368 * @param cellDim the requested elements dimension (if exists) 00369 * @param cellID , ID of the new hanging cell 00370 * @param usedge , whould we look for that element in edge or in face 00371 * @param parentID , ID of the parent edge or face 00372 * @param parentOffset , the offset in the parent */ 00373 void addHangingElement(int cellDim, int cellID ,bool useEdge, int parentID , int parentOffset); 00374 00375 /** 00376 * @param cellDim 00377 * @param cellID */ 00378 int numMaxCofacets_ID(int cellDim, int cellID); 00379 00380 /** Number of processors */ 00381 int nrProc_; 00382 int myRank_; 00383 /** The communicator */ 00384 const MPIComm& _comm; 00385 /** */ 00386 double _pos_x; 00387 /** */ 00388 double _pos_y; 00389 /** */ 00390 double _pos_z; 00391 /** */ 00392 double _ofs_x; 00393 /** */ 00394 double _ofs_y; 00395 /** */ 00396 double _ofs_z; 00397 /** */ 00398 int _res_x; 00399 /** */ 00400 int _res_y; 00401 /** */ 00402 int _res_z; 00403 /** */ 00404 mutable RefinementClass refineClass_; 00405 /** */ 00406 mutable MeshDomainDef meshDomain_; 00407 00408 //------ Point storage ---- 00409 /** all the global points index is [ID] */ 00410 Array<Point> points_; 00411 /** [4] the nr of ID per dim*/ 00412 Array<int> nrElem_; 00413 /** [4] the nr of owned elements per dim*/ 00414 Array<int> nrElemOwned_; 00415 00416 //----- Facets ----- -1 -> MeshBoundary 00417 00418 /** [cellID][8]*/ 00419 Array< Array<int> > cellsPoints_; 00420 /** [cellID][12]*/ 00421 Array< Array<int> > cellsEdges_; 00422 /** [cellID][6]*/ 00423 Array< Array<int> > cellsFaces_; 00424 /** [faceID][4]*/ 00425 Array< Array<int> > faceEdges_; 00426 /** [faceID][4]*/ 00427 Array< Array<int> > facePoints_; 00428 /** [edgeID][2]*/ 00429 Array< Array<int> > edgePoints_; 00430 00431 /** stores the edge orientation {0,1,2} 00432 * [edgeID]*/ 00433 Array<short int> edgeOrientation_; 00434 /** stores the face orientation {0,1,2} 00435 * [faceID]*/ 00436 Array<short int> faceOrientation_; 00437 00438 // ----- MaxCofacets ----; -1 -> MeshBoundary , -2 -> ProcBoundary 00439 00440 /** [faceID][2] */ 00441 Array< Array<int> > faceMaxCoF_; 00442 /** [edgeID][4] */ 00443 Array< Array<int> > edgeMaxCoF_; 00444 /** [pointID][8]*/ 00445 Array< Array<int> > pointMaxCoF_; 00446 00447 //------ Element (processor) ownership ----- 00448 00449 /** contains the ownership of the local elements [dim][ID] */ 00450 Array< Array<short int> > elementOwner_; 00451 00452 //---- hierarchical storage ----- 00453 00454 /** [cellID] , the child index in the parent */ 00455 Array<short int> indexInParent_; 00456 /** [cellID] , the LID of the parent cell */ 00457 Array<int> parentCellLID_; 00458 /** [cellID] , actual level of the cell*/ 00459 Array<short int> cellLevel_; 00460 /** [cellID] , if the element is leaf */ 00461 Array<bool> isCellLeaf_; 00462 /** [cellID] , if the cell is complete outside the user defined mesh domain*/ 00463 Array<bool> isCellOut_; 00464 /** [cellID] , children of the cell*/ 00465 Array< Array<int> > cellsChildren_; 00466 00467 // ---- "hanging" info storage --- 00468 /** [pointID] , true if the node is hanging , false otherwise */ 00469 Array<bool> isPointHanging_; 00470 /** [edgeID] , true if the edge is hanging , false otherwise*/ 00471 Array<bool> isEdgeHanging_; 00472 /** [faceID] , true if the face is hanging , false otherwise*/ 00473 Array<bool> isFaceHanging_; 00474 00475 // ---- hanging element and refinement (temporary) storage --- 00476 00477 /** [edgeID] - > { h P0 LID , h P1 LID , h E0 LID , h E1 LID , h E2 LID } <br> 00478 * These elements will only be put on not hanging when we access it 3 times <br> 00479 * together 5 elements*/ 00480 Hashtable< int, Array<int> > edgeHangingElmStore_; 00481 00482 /** the counter for each edge which stores hanging node information, when the counter reaches 2 (3 read access)<br> 00483 * then the hanging elements can be marked as not hanging <br> 00484 * At the beginning this should be equal with the numMaxCofacets */ 00485 Array<short int> hangingAccessCount_; 00486 00487 /** [faceID] - > { h P0 LID , h P1 LID , h P2 LID , h P3 LID , <br> 00488 * h E0 LID , h E1 LID , h E2 LID , h E3 LID , h E4 LID , h E5 LID , (4+)<br> 00489 * h E6 LID , h E7 LID , h E8 LID , h E9 LID , h E10 LID, h E11 LID , (16+) <br> 00490 * h F0 LID , h F1 LID , h F2 LID , h F3 LID , h F4 LID , h F5 LID , <br> 00491 * h F6 LID , h F7 LID , h F8 LID } <br> 00492 * together 16+9 = 25 elements*/ 00493 Hashtable< int, Array<int> > faceHangingElmStore_; 00494 00495 /** Neighbor Cell can mark the cell to provoke refinement */ 00496 Array<short int> refineCell_; 00497 00498 // ---- leaf mapping , GID and LID --- (points need this because) 00499 00500 /** [leaf_vertexLID] , the value must be a positive number */ 00501 Array<int> vertexLeafToLIDMapping_; 00502 /** [leaf_edgeLID] , the value must be a positive number */ 00503 Array<int> edgeLeafToLIDMapping_; 00504 /** [leaf_faceLID] , the value must be a positive number */ 00505 Array<int> faceLeafToLIDMapping_; 00506 /** [leaf_cellLID] , the value must be a positive number */ 00507 Array<int> cellLeafToLIDMapping_; 00508 /** [vertexLID] if vertex is inside the domain then > 0 , -1 otherwise */ 00509 Array<int> vertexLIDToLeafMapping_; 00510 /** [edgeLID] if edge is leaf(or inside the domain) then > 0 , -1 otherwise */ 00511 Array<int> edgeLIDToLeafMapping_; 00512 /** [faceLID] if face is leaf(or inside the domain) then > 0 , -1 otherwise */ 00513 Array<int> faceLIDToLeafMapping_; 00514 /** [cellLID] if cell is leaf(or inside the domain) then > 0 , -1 otherwise */ 00515 Array<int> cellLIDToLeafMapping_; 00516 00517 /** leaf LID numbering */ 00518 int nrVertexLeafLID_; 00519 int nrEdgeLeafLID_; 00520 int nrFaceLeafLID_; 00521 int nrCellLeafLID_; 00522 00523 /** [leaf_vertexGID] , the value must be a positive number */ 00524 Array<int> vertexLeafToGIDMapping_; 00525 /** [leaf_edgeGID] , the value must be a positive number */ 00526 Array<int> edgeLeafToGIDMapping_; 00527 /** [leaf_faceGID] , the value must be a positive number */ 00528 Array<int> faceLeafToGIDMapping_; 00529 /** [leaf_cellGID] , the value must be a positive number */ 00530 Array<int> cellLeafToGIDMapping_; 00531 /** [vertexGID] if vertex is inside the domain then > 0 , -1 otherwise */ 00532 Array<int> vertexGIDToLeafMapping_; 00533 /** [edgeGID] if edge is leaf(or inside the domain) then > 0 , -1 otherwise */ 00534 Array<int> edgeGIDToLeafMapping_; 00535 /** [faceGID] if edge is leaf(or inside the domain) then > 0 , -1 otherwise */ 00536 Array<int> faceGIDToLeafMapping_; 00537 /** [cellGID] if cell is leaf(or inside the domain) then > 0 , -1 otherwise */ 00538 Array<int> cellGIDToLeafMapping_; 00539 00540 /** leaf GID numbering */ 00541 int nrVertexLeafGID_; 00542 int nrEdgeLeafGID_; 00543 int nrFaceLeafGID_; 00544 int nrCellLeafGID_; 00545 00546 // ------------- static data ------------------- 00547 00548 /** the offset in the X coordinate on the reference cell*/ 00549 static int offs_Points_x_[8]; 00550 00551 /** the offset in the Y coordinate on the reference cell*/ 00552 static int offs_Points_y_[8]; 00553 00554 /** the offset in the Z coordinate on the reference cell*/ 00555 static int offs_Points_z_[8]; 00556 00557 /** stores the facet information on the reference Cell*/ 00558 static int edge_Points_localIndex[12][2]; 00559 00560 /** the edge orientation (local orientation)*/ 00561 static int edge_Orientation[12]; 00562 static int edge_MaxCofacetIndex[3][4]; 00563 static int edge_MaxCof[12]; 00564 00565 /** face edge-facet information */ 00566 static int face_Edges_localIndex[6][4]; 00567 00568 /** face point-facet information */ 00569 static int face_Points_localIndex[6][4]; 00570 00571 /** face orientation (local orientation)*/ 00572 static int face_Orientation[6]; 00573 static int face_MaxCofacetIndex[3][2]; 00574 static int face_MaxCof[6]; 00575 00576 /** used for coarse mesh creation for the vertex indexing */ 00577 static int vInd[8]; 00578 00579 /** used for coarse mesh creation for the edge indexing */ 00580 static int eInd[12]; 00581 00582 /** used for coarse mesh creation for the face indexing */ 00583 static int fInd[6]; 00584 00585 // the X and the Y coordinates of the newly 00586 static double vertex_X[64]; 00587 00588 static double vertex_Y[64]; 00589 00590 static double vertex_Z[64]; 00591 00592 // face index is above 20 00593 static int vertexToParentEdge[64]; 00594 // 00595 static int vertexInParentIndex[64]; 00596 // 00597 static int edgeToParentEdge[144]; 00598 // 00599 static int edgeInParentIndex[144]; 00600 // 00601 static int faceToParentFace[108]; 00602 // 00603 static int faceInParentIndex[108]; 00604 00605 }; 00606 } 00607 00608 00609 #endif /* SUNDANCEHNMESH3D_H_ */