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_MESHBASE_H 00032 #define SUNDANCE_MESHBASE_H 00033 00034 #include "SundanceDefs.hpp" 00035 #include "SundancePoint.hpp" 00036 #include "SundanceSet.hpp" 00037 #include "SundanceMap.hpp" 00038 #include "SundanceCellType.hpp" 00039 #include "SundanceObjectWithVerbosity.hpp" 00040 #include "SundanceObjectWithInstanceID.hpp" 00041 #include "PlayaMPIComm.hpp" 00042 #include "Teuchos_Array.hpp" 00043 #include "SundanceCellReorderer.hpp" 00044 #include "SundanceCellReordererImplemBase.hpp" 00045 00046 namespace Sundance { 00047 00048 00049 /** Identifier for ordering convention */ 00050 enum MeshEntityOrder {UFCMeshOrder, ExodusMeshOrder}; 00051 00052 class MaximalCofacetBatch; 00053 00054 using namespace Teuchos; 00055 using Playa::MPIComm; 00056 00057 class CellJacobianBatch; 00058 00059 00060 /** \brief Abstract interface to a mesh. 00061 * 00062 * \section Sundance_MeshBase_Outline_sec Outline 00063 * 00064 * <ul> 00065 * <li>\ref Sundance_MeshBase_Introduction_sec 00066 * <li>\ref Sundance_MeshBase_Definitions_sec 00067 * <li>\ref Sundance_MeshBase_Common_Function_Arguments_sec 00068 * <li>\ref Sundance_MeshBase_Refactorings_sec 00069 * </ul> 00070 * 00071 * \section Sundance_MeshBase_Introduction_sec Introduction 00072 * 00073 * A mesh is a collection of connected cells (also called elements). This 00074 * mesh interface is designed to allow for the representation of 1D, 2D, and 3D (and 00075 * even higher dimensional) meshes. A mesh contains both topological 00076 * information and geometric information about the cells and their subparts 00077 * (i.e. facets). Topological information refers to information about how 00078 * cells of various types and dimensions are connected to each other within 00079 * the mesh. Geometric information refers to information about the position, 00080 * size, and shape of cells within the mesh. 00081 * 00082 * Currently, this interface only handles meshes where all cells of a given 00083 * cell dimension have the same cell type. For example, this only allows 00084 * meshes with all triangles or all quads in 2D or all tetrahedrals in 3D. 00085 * 00086 * \todo Add some diagrams to show different cell types to help define cells, 00087 * facets, co-facets, vertices, etc. 00088 * 00089 * \section Sundance_MeshBase_Definitions_sec Definitions 00090 * 00091 * The following definitions are used in the specification of this mesh 00092 * interface. 00093 * 00094 *<ul> 00095 * 00096 * <li><b>Topological information</b> 00097 * 00098 * <ul> 00099 * 00100 * <li><b>Spatial Dimension (spatialDim)</b>: The maximum cell dimension of 00101 * any cell type in the mesh. For example, a 3D mesh will have 00102 * spatialDim=3, a 2D mesh will have spatialDim=2 and a 1D mesh will have 00103 * spatialDim=1. Typically, when we use the short hand 1D, 2D or 3D, we are 00104 * referring to the spatial dimension of the mesh. Note that 4D and 00105 * higher-dimensional meshes are also representable by this interface in 00106 * principle with some small changes (e.g. like making the size of a point 00107 * more flexible) 00108 * 00109 * <li><b>Cell</b>: A mesh object of some dimension. For example, in a 3D 00110 * mesh, cells are things like elements, faces, edges, and nodes. 00111 * 00112 * <li><b>Cell Type (cellType)</b>: An enumeration of type <tt>CellType</tt> 00113 * that lists some common cell types. For example, common cell types are 00114 * things like vertexes (0D), lines or edges (1D), triangles (2D), and 00115 * quadrilaterals (2D). 00116 * 00117 * <li><b>Cell Dimension (cellDim)</b>: The dimension of a cell. More 00118 * formally, the cell dimension is equal to the number of lower level facet 00119 * types in the cell. For example, A 3D tetrahedral cell (or element) has 00120 * cellDim=3 and it has the three lower-level facet types faces (cellDim=2), 00121 * edges (cellDim=1), and vertexes (cellDim=0). 00122 * 00123 * <li><b>Relative Cell</b>: The term "relative cell" is used below to refer 00124 * to the cell object for which we consider other cells called "facets" and 00125 * "co-facets". 00126 * 00127 * <li><b>Facet</b>: A lower dimensional cell directly connected to and 00128 * contained within a relative <em>parent</em> cell. For example, a 3D cell 00129 * (element) has three different facet types: faces (facetDim=2), edges 00130 * (facetDim=1), and vertexes (facetDim=0). A face (cellDim=2) relative 00131 * cell (or a 2D element) has two different facet types: edges (facetDim=1), 00132 * and vertexes (facetDim=0). A particular relative cell 00133 * <tt>(cellDim,cellLID)</tt> has 00134 * <tt>nf=numFacets(cellDim,cellLID,facetDim)</tt> facets of a particular 00135 * dimension <tt>facetDim</tt>. 00136 * 00137 * <li><b>Facet Dim (facetDim)</b>: The cell dimension of a given facet. 00138 * For example, the dimension of a face facet of a 3D element is facetDim=2. 00139 * 00140 * <li><b>Facet Index (facetIndex)</b>: This is the local index of the facet 00141 * with respect to its relative parent cell. For example, a 3D brick 00142 * relative cell/element will have 6 faces which can be accessed using the 00143 * facet indexes 0, 1, 2, 3, 4, and 5. 00144 * 00145 * <li><b>Co-facet</b>: A higher dimensional cell of which the given 00146 * relative cell is a facet. For example, the co-facets of a relative 00147 * face cell in a 3D mesh are the one or two 3D element cells that share 00148 * this face. In the case of a boundary face in a 3D mesh, only one 00149 * co-facet is connected to the boundary cell and that is its boundary 00150 * element. Note that in 3D, edges can have many different face and element 00151 * co-facets. There are some important properties of co-facets. Relative 00152 * cells of dimension cellDim=spatialDim-1 will have either one or two 00153 * co-facets of dimension cofacetDim=spatialDim. Therefore, relative cells 00154 * of the dimension cellDim=spatialDim-1 that have only one co-facet must be 00155 * a boundary cell (i.e. a boundary face in 3D, a boundary edge in 2D, or a 00156 * boundary vertex in 1D). Note that relative cells with dimension 00157 * cellDim=spatialDim-n, where n >=2, can have more than two co-facets in 00158 * general. 00159 * 00160 * <li><b>Co-facet Dimension (cofacetDim)</b>: The cell dimension of a 00161 * co-facet cell. 00162 * 00163 * <li><b>Maximal co-facet</b>: A co-facet of a relative cell whose 00164 * dimension is equal to the spatial dimension of the mesh. 00165 * 00166 * <em>Assertion</em>: <tt>maxCofacetDim == spatialDim</tt>. 00167 * 00168 * <li><b>Vertex</b>: A zero dimensional cell (cellDim=0). 00169 * 00170 * <em>Assertion</em>: <tt>vertexCellDim == 0</tt>. 00171 * 00172 * <li><b>Element</b>: The term "element" is used to refer to cells whose 00173 * dimension is equal to the spatial dimension of the mesh. Warning, the 00174 * term "element" is an overloaded term and is also used to refer to an 00175 * entry in an array of objects. 00176 * 00177 * <li><b>Facet orientation (facetOrientation)</b>: The orientation of a 00178 * facet with respect to its parent cell. This is given as an integer value 00179 * that requires agreement of interpretation. In 2D, an edge for instance 00180 * will have just a positive and a negative orientation. In 3D, a face can 00181 * have more than two orientations and each orientation is a permutation of 00182 * the different nodal listings on the face. The facet orientation integer 00183 * value is typically interpreted as an index into an array of possible 00184 * orientations. The mesh reader/creator and the assembly code have to 00185 * agree on how these orientations are interpreted. 00186 * 00187 * <li><b>Co-facet index (cofacetIndex)</b>: The relative index of a 00188 * co-facet of some relative cell. For example, an internal face in a 3D 00189 * mesh will have two co-facet element objects which can be accessed using 00190 * the co-facet indexes 0 and 1. 00191 * 00192 * <li><b>Label</b>: Every cell can be given a single integer value that 00193 * can be used in a variety of ways. For example, boundary faces in 3D can 00194 * be given an integer label to specify different types of boundary 00195 * conditions. This integer could also be used with bit manipulation to 00196 * allow for multiple classifications for a single cell. 00197 * 00198 * <li><b>Intermediate Cell</b>: An intermediate cell is a cell that has a 00199 * dimension cellDim in the range 0 < cellDim < spatialDim. For example, in 00200 * a 3D mesh, faces and edges are intermediate cell types. 00201 * 00202 * </ul> 00203 * 00204 * <li><b>Geometric information</b> 00205 * 00206 * <ul> 00207 * 00208 * <li><b>Point</b>: A point is an array of coordinate positions. Currently 00209 * this is stored as a array with three doubles (i.e. x=0, y=1, and z=2). 00210 * For 2D meshes only the x=0 and y=1 positions are used. For 1D meshes, 00211 * only the x=0 position is used. 00212 * 00213 * <li><b>Cell Centroid</b>: A point that represents the geometric center of 00214 * a cell. 00215 * 00216 * <li> <b>Cell Diameter</b>: The cell diameter is defined as the maximum 00217 * line segment that can be drawn in the cell. This gives a measure of how 00218 * big the cell is and, for instance, can be used in some types of 00219 * finite-element stabilization methods. 00220 * 00221 * <li><b>Cell Curvature</b>: The curvature of cell refers to the shape of a 00222 * cell. This is typically used to represent the shapes of faces in 3D or 00223 * edges in 2D which is critical in the implementation of higher-order 00224 * higher-order discrimination methods. Cell curvature is currently not 00225 * supported by this mesh interface but could be in the future. 00226 * 00227 * <li><b>Reference Cell</b>: The term "reference cell" is used below to 00228 * refer that a cell type's reference cell which is a geometric description 00229 * of the cell in a simple coordinate system. By using a standard 00230 * definition of a "reference cell", we can decouple different types of 00231 * mathematical operations defined of a cell of a certain basic type (e.g. 00232 * triangles and quadrilaterals in 2D, and tetrahedral and bricks in 3D). 00233 * Typically, the nodes of a reference cell will have coordinates that 00234 * integer multiples of one (including zero). 00235 * 00236 * <li><b>Cell Jacobian</b>: The cell Jacobian is a transformation matrix 00237 * that maps cell vertex coordinate values (i.e. points) from the reference 00238 * cell to the physical cell. The inverse of the cell Jacobian maps points 00239 * from the physical cell to the reference cell. A cell with the cell 00240 * dimension cellDim has a cell Jacobian matrix of dimension cellDim x 00241 * cellDim. 00242 * 00243 * </ul> 00244 * 00245 * <li><b>Parallel Information</b> 00246 * 00247 * <ul> 00248 * 00249 * <li><b>Cell Ownership</b>: The ownership of a cell refers to which 00250 * process owns the cell. All other processes that need to access that cell 00251 * get a "ghost" copy of the cell locally. 00252 * 00253 * <li><b>Ghost Cell</b>: A ghost cell is a cell that is not owned by the 00254 * local process that is copied to a process in order to perform certain 00255 * computations. 00256 * 00257 * <li><b>Cell Local ID (cellLID)</b>: The LID of a cell is an index to 00258 * the cell within the local process and these IDs are typically ordered 00259 * from 0 to numLocalCells-1. Local cells can include both owned cells and 00260 * ghosted cells. 00261 * 00262 * <li><b>Cell Global ID (cellGID)</b>: The global ID of a cell is an index 00263 * that is unique across all processes and is used as a process-independent 00264 * identifier for the cell. 00265 * 00266 * </ul> 00267 * 00268 * </ul> 00269 * 00270 * \section Sundance_MeshBase_Common_Function_Arguments_sec Common Function Arguments and Pre/Post-Conditions 00271 * 00272 * Below, a common set of Mesh arguments are described along with relevant 00273 * pre- and post conditions. 00274 * 00275 * <ul> 00276 * 00277 * <li><b>cellDim</b> [in] The dimension of a cell 00278 * 00279 * <em>Precondition</em>: <tt>0 <= cellDim <= spatialDim</tt> 00280 * 00281 * <li><b>cellLID(s)</b> [in] Local ID (LID) of a cell(s). 00282 * 00283 * <em>Precondition</em>: <tt>0 <= cellLID < numCells(cellDim)</tt>, where 00284 * <tt>cellDim</tt> is the dimension of the given cell. 00285 * 00286 * <li><b>cellGID</b> [in] Global ID (GID) of a cell which is unique across 00287 * processes. 00288 * 00289 * <em>Precondition</em>: <tt>??? <= cellGID < ???</tt> (How do we write 00290 * this?). 00291 * 00292 * <li><b>facetDim</b> [in] The dimension of a given facet cell. 00293 * 00294 * <em>Precondition</em>: <tt>0 <= facetDim < cellDim</tt>, where 00295 * <tt>cellDim</tt> is the dimension of the relative cell 00296 * 00297 * <li><b>facetIndex</b> [in] Local index of the facet with respect to its 00298 * parent relative cell. 00299 * 00300 * <em>Precondition</em>: <tt>0 <= facetIndex < 00301 * numFacets(cellDim,cellLID,facetDim)</tt>, where <tt>facetDim</tt> is the 00302 * dimension of the facet of the relative parent cell 00303 * <tt>(cellDim,cellLID)</tt> 00304 * 00305 * <li><b>facetOrientation(s)</b> [out] The orientation of a facet with 00306 * respect to its parent cell. This is given as an integer value that 00307 * requires agreement of interpretation (see above). 00308 * 00309 * <li><b>cofacetDim</b> [in] Cell dimension of a co-facet cell. 00310 * 00311 * <em>Precondition</em>: <tt>0 <= cellDim < cofacetDim <= spatialDim</tt>, 00312 * where <tt>cellDim</tt> is the dimension of the relative cell 00313 * 00314 * <li><b>cofacetIndex</b> [in] Relative index of a co-facet of some relative 00315 * cell. 00316 * 00317 * <em>Precondition</em>: <tt>0 <= cofacetIndex < 00318 * numMaxCofacets(cellDim,cellLID)</tt>, where <tt>(cellDim,cellLID)</tt> is 00319 * the relative cell. 00320 * 00321 * </ul> 00322 * 00323 * \section Sundance_MeshBase_Refactorings_sec Possible Refactorings 00324 * 00325 * <ul> 00326 * 00327 * <li>Copy ObjectWithInstanceID base class into a special Nihilo tools 00328 * package or into Teuchos. 00329 * 00330 * <li>Inherit MeshBase from Teuchos::VerboseObject instead of 00331 * ObjectWithVerbosity. 00332 * 00333 * <li>Remove any concrete implementations and data from this base interface 00334 * to provide a clean abstract interface. This would include removing the 00335 * constructor. 00336 * 00337 * <li>Either remove the non-batch access functions all together, or make the 00338 * namings of batch and non-batch access functions more consistent. 00339 * 00340 * <li>Have all of the functions return just a single return value or when 00341 * more than one object is returned, return all values from the argument list. 00342 * Or, when the value that is being returned from the argument list is only 00343 * occasionally returned, change it to a pointer type and give it a default 00344 * value of NULL so that it can be safely ignored! This will actually speed 00345 * up performance in some cases. For example, <tt>facetOrientation</tt> is a 00346 * return argument that is often ignored and would benefit from making it an 00347 * optional argument. 00348 * 00349 * <li>Remove the staggered output function since the Teuchos::FancyOStream 00350 * should handle parallel output from every process well. 00351 * 00352 * <li>Add a new cell type called a "general" cell for very irregular cells 00353 * that defy fixed classification. This could be used to represent many of 00354 * the irregular cells that remain after non-conforming mesh refinement. 00355 * 00356 * <li>Make sure to tell the DOF map how to interpret the facet orientation 00357 * integers. Currently the interpretation is hard coded in Sundance. These 00358 * permutations need to be enumerated in a C++ enum and the DOF map must be 00359 * given the array of these enums to interpret the orientation integer values. 00360 * To be 100 percent safe, we need a way of allowing the mesh to directly 00361 * return the interpretation of these orientations. What you can do is to 00362 * require that the int value be an index into an an array of enumeration 00363 * values that corresponds to the cell type. This would require an enum type 00364 * and an array of enum values for every facet cell type. For the current 00365 * cell types for regular elements, this would be no problem. However, for 00366 * the facets of a "general" cell type, the facets may also be classified as 00367 * "general" cell types and therefore may not allow such an orientation array 00368 * to be defined. 00369 * 00370 * <li>Clearly define the coordinate system of the reference cell of each of 00371 * the standard cell types in order to clearly define interpretation of the 00372 * reference points passed into the <tt>pushForward()</tt> function. Without 00373 * a clear definition of the reference cell, these reference points are 00374 * meaningless. Actually, for maximal flexibility, each cell type should 00375 * define, at runtime, the coordinate system for the reference cell. Or, we 00376 * must just all agree on a single set of standard definitions of reference cells 00377 * for different cell types. 00378 * 00379 * </ul> 00380 */ 00381 class MeshBase 00382 : public ObjectWithClassVerbosity<MeshBase> 00383 , public ObjectWithInstanceID<MeshBase> 00384 { 00385 public: 00386 00387 /** \name Constructors/Destructor */ 00388 //@{ 00389 00390 /** \brief . */ 00391 MeshBase(int dim, const MPIComm& comm, 00392 const MeshEntityOrder& meshOrder); 00393 00394 /** \brief . */ 00395 virtual ~MeshBase(){;} 00396 00397 //@} 00398 00399 /** \name Topological Information */ 00400 //@{ 00401 00402 /** \brief Get the ordering convention used by this mesh */ 00403 const MeshEntityOrder& meshOrder() const {return order_;} 00404 00405 /** \brief Get the spatial dimension of the mesh 00406 * 00407 * <b>Postconditions:</b><ul> 00408 * <li><tt>0 < returnVal <= 3</tt> 00409 * </ul> 00410 */ 00411 int spatialDim() const {return dim_;} 00412 00413 /** \brief Get the number of local cells having dimension cellDim. 00414 * 00415 * \todo Change this to <tt>numLocalCells(cellDim)</tt>. 00416 */ 00417 virtual int numCells(int cellDim) const = 0 ; 00418 00419 /** \brief Return the number of facets of the given relative cell. 00420 * 00421 * \param cellDim 00422 * [in] The dimension of the relative cell 00423 * \param cellLID 00424 * [in] The LID of the relative cell 00425 * \param facetDim 00426 * [in] The dimension of the facets for the 00427 * relative cell 00428 * 00429 * <b>Postconditions:</b><ul> 00430 * <li><tt>returnVal >= 2</tt>: Every cell has two or more 00431 * facets of a given dimension. 00432 * </ul> 00433 */ 00434 virtual int numFacets(int cellDim, int cellLID, 00435 int facetDim) const = 0 ; 00436 00437 /** \brief Return the LID of a single facet cell (and optionally its 00438 * orientation) with respect to a relative parent cell. 00439 * 00440 * \param cellDim 00441 * [in] Dimension of the parent cell whose facets 00442 * are being obtained. 00443 * \param cellLID 00444 * [in] Local ID of the parent cell whose facets 00445 * are being obtained 00446 * \param facetDim 00447 * [in] Dimension of the desired facet. 00448 * \param facetIndex 00449 * [in] The relative index into the list of the 00450 * cell's facets. 00451 * \param facetOrientation 00452 * [out] The orientation of the facet w.r.t the 00453 * parent cell. 00454 * 00455 * \todo Change the facetOrientation argument to be a raw pointer that can 00456 * be NULL and therefore easily ignored. 00457 */ 00458 virtual int facetLID(int cellDim, int cellLID, 00459 int facetDim, int facetIndex, 00460 int& facetOrientation) const = 0 ; 00461 00462 /** \brief Return an array containing the LIDs of facets of dimension 00463 * facetDim for the given batch of relative parent cells. 00464 * 00465 * \param cellDim 00466 * [in] Dimension of the relative parent cells whose facets 00467 * are being obtained 00468 * \param cellLIDs 00469 * [in] Array of LIDs for the relative parent cells 00470 * \param facetDim 00471 * [in] Dimension of the desired facets 00472 * \param facetLIDs 00473 * [out] On exit this array gives the local facet IDs 00474 * for all of the cells given in cellLIDs in one flat array 00475 * with size = <tt>cellLIDs.size()*nf</tt>, where 00476 * <tt>nf=numFacets(cellDim,cellLIDs[j],facetDim)</tt>) where 00477 * <tt>j</tt> can be any index <tt>0 <= j < numCells(cellDim)</tt>. 00478 * Specifically, the local facet IDs for <tt>cellLIDs[k]</tt>, where 00479 * <tt>k=0...cellLIDs.size()-1</tt>, is given 00480 * by <tt>facetLIDs[k*nf+j]</tt>, where <tt>j=0...nf-1</tt>. 00481 * \param facetOrientations 00482 * [out] On exist, if <tt>facetDim > 0</tt>, this array gives 00483 * the integer orientation of the facet with respect to 00484 * its parent cell (see above definition of "Facet Orientation"). 00485 * Oon exist this array will be 00486 * resized to size = <tt>cellLIDs.size()*nf</tt>, where 00487 * <tt>nf=numFacets(cellDim,cellLID,facetDim)</tt>). 00488 * Specifically, the local facet orientation for the cell cellLIDs[k], where 00489 * <tt>k=0...cellLIDs.size()-1</tt>, is given 00490 * by <tt>facetOrientations[k*nf+j]</tt>, where <tt>j=0...nf-1</tt>. 00491 * If <tt>facetDim==0</tt> then this array is ignored. 00492 * 00493 * <b>Warning!</b> This function will only work for homogeneous cell types! 00494 * 00495 * \todo Change the facetOrientation argument to an <tt>Array<int>*</tt> 00496 * type so that it can be NULL and therefore ignored (which is a common use 00497 * case). 00498 */ 00499 virtual void getFacetLIDs(int cellDim, 00500 const Array<int>& cellLIDs, 00501 int facetDim, 00502 Array<int>& facetLIDs, 00503 Array<int>& facetOrientations) const = 0 ; 00504 00505 /** \brief Return an array containing the LIDs of the facets of 00506 * dimension facetDim of the single relative parent cell (optionally also 00507 * returning the facet orientations). 00508 * 00509 * \param cellDim 00510 * [in] Dimension of the parent cell whose facets 00511 * are being obtained. 00512 * \param cellLID 00513 * [in] Local ID of the parent cell whose facets 00514 * are being obtained 00515 * \param facetDim 00516 * [in] Dimension of the desired facets 00517 * \param facetLIDs 00518 * [out] On exit this array gives the local facet IDs 00519 * for the parent cell 00520 * with size = <tt>nf</tt>, where 00521 * <tt>nf=numFacets(cellDim,cellLID,facetDim)</tt>. 00522 * \param facetOrientations 00523 * [out] On exist, if <tt>facetDim > 0</tt>, this array gives 00524 * the integer orientation of the facet with respect to 00525 * its parent cell (see above definition of "Facet Orientation"). 00526 * On exist this array will be 00527 * resized to size = <tt>nf</tt>, where 00528 * <tt>nf=numFacets(cellDim,cellLID,facetDim)</tt>). 00529 * If <tt>facetDim==0</tt> this this array argument is ignored! 00530 * 00531 * The default implementation loops over calls to 00532 * <tt>facetLID()</tt>. Subclasses can provide a more efficient 00533 * implementation if desired. 00534 * 00535 * \todo Rename to something like getSingleCellFacetLIDs(...). 00536 * 00537 * \todo Change the facetOrientation argument to an <tt>Array<int>*</tt> 00538 * type so that it can be NULL and therefore ignored (which is a common use 00539 * case). 00540 */ 00541 void getFacetArray(int cellDim, int cellLID, int facetDim, 00542 Array<int>& facetLIDs, 00543 Array<int>& facetOrientations) const ; 00544 00545 /** \brief Return a view of an array of LIDs for a maximum-dimensional 00546 * cell's zero-dimensional facets (i.e. vertexes). 00547 * 00548 * \param maxCellLID 00549 * [in] Local ID of the maximum dimensional cell (i.e. element). 00550 * 00551 * <b>Preconditions:</b><ul> 00552 * <li>The cell <tt>maxCellLID</tt> must be a maximum-dimensional cell! 00553 * </ul> 00554 * 00555 * \returns A raw pointer into an array which stores the local process IDs 00556 * of the vertices. These IDs are the local process IDs, not the facet 00557 * indexes relative to the relative parent cell. The dimension of this 00558 * array is determined by the cell type given by 00559 * <tt>cellType(maxCellLID)</tt>. 00560 * 00561 * \todo Return this array as a Teuchos::ArrayView<int> object which will 00562 * have the dimension embedded in it and will have full range checking! 00563 * 00564 * \todo Rename to something like getMaxCellZeroFacetsLIDsView(...). 00565 */ 00566 virtual const int* elemZeroFacetView(int maxCellLID) const ; 00567 00568 /** \brief Return the number of maximal co-facets of the given cell. 00569 * 00570 * <b>Preconditions:</b><ul> 00571 * <li><tt>0 <= cellDim < spatialDim()</tt> 00572 * </ul> 00573 * 00574 * Note that if <tt>cellDim==spatialDim()-1</tt> and <tt>returnVal==1</tt> then 00575 * this cell must be a boundary cell (i.e. a boundary face in 3D, a boundary 00576 * edge in 1D, or a boundary node in 1D)! 00577 */ 00578 virtual int numMaxCofacets(int cellDim, int cellLID) const = 0; 00579 00580 /** \brief Return the LID of a maximal co-facet of a cell. 00581 * 00582 * \param cellDim 00583 * [in] Dimension of the cell whose co-facets are being obtained 00584 * \param cellLID 00585 * [in] Local ID of the cell whose co-facets are being obtained 00586 * \param cofacetIndex 00587 * [in] Index into the list of the cell's co-facets. 00588 * \param facetIndex 00589 * [out] Returns the local index into the facet array 00590 * for the relative cell (cellDim,cellLID) with respect 00591 * to the maximal co-facet. In other words 00592 * <tt>cellLID==facetLID(spatialDim(),returnVal,cellDim,facetIndex)</tt>. 00593 * 00594 * <b>Preconditions:</b><ul> 00595 * <li><tt>0 <= cellDim < spatialDim()</tt> 00596 * </ul> 00597 * 00598 * \returns The LID of the requested maximal co-facet. 00599 * 00600 * \todo Make the <tt>facetIndex</tt> an <tt>int*</tt> argument and give it a 00601 * default value of <tt>NUL</tt> so that it can be easily ignored! 00602 */ 00603 virtual int maxCofacetLID(int cellDim, int cellLID, 00604 int cofacetIndex, 00605 int& facetIndex) const = 0 ; 00606 /** 00607 * Get the LIDs of the maximal cofacets for a batch of codimension-one 00608 * cells. The default implementation simply loops over the cells in the 00609 * cellLID array, taking no advantage of any internal data structures. 00610 * 00611 * \param cellLIDs [in] array of LIDs of the cells whose cofacets are 00612 * being obtained 00613 * \param cofacets [out] 00614 * \param facetIndex [out] index of each calling cell 00615 * into the list of its maximal cofacet's facets 00616 */ 00617 virtual void getMaxCofacetLIDs(const Array<int>& cellLIDs, 00618 MaximalCofacetBatch& cofacets) const ; 00619 00620 /** \brief Return an array of the LIDs of all of the co-facets for a 00621 * given relative cell. 00622 * 00623 * \param cellDim 00624 * [in] Dimension of the relative cell whose co-facets are being obtained 00625 * \param cellLID 00626 * [in] Local index of the relative cell whose co-facets are being obtained 00627 * \param cofacetDim 00628 * [in] Dimension of the co-facets to get 00629 * \param cofacetLIDs 00630 * [out] Array containing the LIDs of the co-facets 00631 * for the given relative cell (cellDim,cellLID). 00632 * 00633 * \todo Change name to <tt>getCofacetArray()</tt> to be consistent with 00634 * <tt>getFacetArray()</tt>! 00635 */ 00636 virtual void getCofacets(int cellDim, int cellLID, 00637 int cofacetDim, Array<int>& cofacetLIDs) const = 0 ; 00638 00639 /** \brief Get the cell type of the given cell dimension. 00640 * 00641 * Note: This function by its very definition assumes that all cells of a 00642 * given dimension have the same cell type within a mesh! 00643 * 00644 * \todo This function must be changed in order to deal with mixed cell 00645 * types with the same cellDim! 00646 */ 00647 virtual CellType cellType(int cellDim) const = 0 ; 00648 00649 /** \brief Get the label of the given cell. 00650 * 00651 * \todo Change to getLabel(...)? 00652 */ 00653 virtual int label(int cellDim, int cellLID) const = 0 ; 00654 00655 /** \brief Get the labels for a batch of cells. 00656 * 00657 * \param cellDim 00658 * [in] Dimension of the parent cell whose facets 00659 * are being obtained 00660 * \param cellLIDs 00661 * [in] Array of cell LIDs 00662 * \param labels 00663 * [out] On exit, contains an array (<tt>size=cellLIDs.size()</tt>) 00664 * of the labels for each of the given cells. 00665 */ 00666 void getLabels(int cellDim, const Array<int>& cellLIDs, 00667 Array<int>& labels) const ; 00668 00669 /** \brief Set the label for the given cell. 00670 * 00671 * \todo Move this out of this base interface and into a mesh loading 00672 * interface? 00673 */ 00674 virtual void setLabel(int cellDim, int cellLID, int label) = 0 ; 00675 00676 /** Get the list of all labels defined for cells of the given dimension */ 00677 virtual Set<int> getAllLabelsForDimension(int cellDim) const ; 00678 00679 /** 00680 * Get the cells associated with a specified label. The array 00681 * cellLID will be filled with those cells of dimension cellDim 00682 * having the given label. 00683 */ 00684 virtual void getLIDsForLabel(int cellDim, int label, Array<int>& cellLIDs) const ; 00685 //@} 00686 00687 /** \name Geometric Information */ 00688 //@{ 00689 00690 /** \brief Return the position of a local vertex. 00691 * 00692 * \param vertexLID 00693 * [in] The LID of the vertex. 00694 * 00695 * <b>Preconditions:</b><ul> 00696 * <li><tt>0 <= vertexLID < this->numCells(0)</tt> 00697 * </ul> 00698 * 00699 * \todo Change the name of this function to 00700 * <tt>getVertexPosition(...)</tt>. 00701 */ 00702 virtual Point nodePosition(int vertexLID) const = 0 ; 00703 00704 /** \brief Return a const view into an raw array for the position of a local 00705 * vertex. 00706 * 00707 * \param vertexLID [in] The LID of the vertex 00708 * 00709 * <b>Preconditions:</b><ul> 00710 * <li><tt>0 <= vertexLID < this->numCells(0)</tt> 00711 * </ul> 00712 * 00713 * \returns a raw pointer into an array of doubles of length 00714 * <tt>spatialDim</tt>. 00715 * 00716 * \todo Return this array as a Teuchos::ArrayView<double> object which will 00717 * have the dimension embedded in it and will have full range checking! 00718 * 00719 * \todo Change this function name to <tt>getVertexPositionView()</tt>. 00720 */ 00721 virtual const double* nodePositionView(int vertexLID) const = 0 ; 00722 00723 /** \brief Return the centroid of a cell. 00724 * 00725 * The default implementation just averages the positions of the 00726 * zero-dimensional facets (i.e. vertexes). 00727 */ 00728 virtual Point centroid(int cellDim, int cellLID) const ; 00729 00730 /** \brief Compute (or get) the Jacobians for a batch of cells. 00731 * 00732 * \param cellDim 00733 * [in] Dimension of the cells whose Jacobians are to be computed 00734 * \param cellLID 00735 * [in] Local IDs of the cells for which Jacobians are to be computed 00736 * \param jBatch 00737 * [out] Batch of cell Jacobians. 00738 * 00739 * Warning! The default implementation returns an empty batch of cell 00740 * Jacobians! 00741 * 00742 * \todo Add a query function to tell if this feature is supported and then 00743 * add a precondition based on this query function! 00744 */ 00745 00746 virtual void getJacobians(int cellDim, const Array<int>& cellLID, 00747 CellJacobianBatch& jBatch) const { ; } 00748 00749 //bvbw virtual void getJacobians( 00750 // int cellDim, const Array<int>& cellLID, 00751 // CellJacobianBatch& jBatch 00752 // ) const 00753 // { jBatch.resize(0,0,0); } 00754 00755 /** \brief Compute the diameters of a batch of cells. 00756 * 00757 * \param cellDim 00758 * [in] Dimension of the cells whose diameters are to be computed 00759 * \param cellLIDs 00760 * [in] Local IDs of the cells for which diameters are to be computed 00761 * \param diameters 00762 * [out] Array (<tt>size = cellLIDs.size()</tt>) of cell diameters. 00763 * 00764 * Warning! The default implementation returns an empty array of cell 00765 * diameters! 00766 * ToDo: Change the default implementation to compute diameters based on 00767 * calls to the node position accessors. Going through the Mesh interface in 00768 * that way will be less efficient than a low-level implementation, but 00769 * would be a reasonable intermediate step for mesh developers. 00770 * - KL 5 Aug 2007. 00771 */ 00772 virtual void getCellDiameters( 00773 int cellDim, const Array<int>& cellLIDs, 00774 Array<double>& diameters 00775 ) const 00776 { diameters.resize(0); } 00777 00778 00779 /** 00780 * Get the outward normals for the batch of cells of dimension 00781 * spatialDim()-1. If any cell in the batch is not on the boundary, 00782 * an exception is thrown. 00783 * 00784 * \param cellLIDs [in] LIDs for the cells whose normals are to be 00785 * computed. 00786 * \param outwardNormals [out] Outward normal unit vectors for each 00787 * cell in the batch. 00788 */ 00789 virtual void outwardNormals( 00790 const Array<int>& cellLIDs, 00791 Array<Point>& outwardNormals 00792 ) const ; 00793 00794 /** 00795 * Get tangent vectors for a batch of edges 00796 * 00797 * \param cellLIDs [in] LIDs for the cells whose tangents are to be 00798 * computed. 00799 * \param tangentVectors [out] Unit tangents for each cell 00800 */ 00801 virtual void tangentsToEdges( 00802 const Array<int>& cellLIDs, 00803 Array<Point>& tangenVectors 00804 ) const ; 00805 00806 00807 /** \brief Map points from a reference cell to physical points for a batch 00808 * of cells. 00809 * 00810 * \param cellDim 00811 * [in] Dimension of the cells 00812 * \param cellLIDs 00813 * [in] Local IDs of a batch of cells 00814 * \param refPts 00815 * [in] Array of points on the single reference cell with respect 00816 * to the reference cell's coordinate system. Note that the 00817 * interpretation of these reference points is strictly determined 00818 * by the coordinate system of the cell type 00819 * <tt>cellType(cellDim)</tt> and must be clearly defined by this 00820 * interface. 00821 * \param physPts 00822 * [out] Array (<tt>size = cellLIDs.size()*refPts.size()</tt>) of 00823 * the physical points given in a flat array for the batch of 00824 * cells. Specifically, the physical points for each cell 00825 * <tt>cellLIDs[k]</tt>, for <tt>k=0...cellLIDs.size()-1</tt>, is 00826 * given by <tt>physPts[k*nrp+j]</tt>, for <tt>j=0...nrp-1</tt> 00827 * and <tt>nrp=refPts.size()</tt>. 00828 * 00829 * Warning! The default implementation returns an empty array of physical 00830 * points! 00831 * 00832 * \todo Add a query function to tell if this feature is supported and then 00833 * add a precondition based on this query function! 00834 */ 00835 virtual void pushForward( 00836 int cellDim, const Array<int>& cellLIDs, 00837 const Array<Point>& refPts, 00838 Array<Point>& physPts 00839 ) const 00840 { physPts.resize(0); } 00841 00842 //@} 00843 00844 /** \name Parallel Information */ 00845 //@{ 00846 00847 /** \brief Return the MPI communicator over which this mesh is distributed. */ 00848 const MPIComm& comm() const {return comm_;} 00849 00850 /** \brief Return the rank of the processor that owns the given cell. 00851 * 00852 * If <tt>returnVal==comm().getRank()</tt> then this cell is owned by this 00853 * process. 00854 */ 00855 virtual int ownerProcID(int cellDim, int cellLID) const = 0 ; 00856 00857 /** \brief Determine whether a given cell GID exists on this processor. */ 00858 virtual bool hasGID(int cellDim, int cellGID) const = 0 ; 00859 00860 /** \brief Find the LID of a cell given its GID. 00861 * 00862 * \param cellDim 00863 * [in] Dimension of the cell 00864 * \param cellGID 00865 * [in] Global ID of the cell 00866 * 00867 * <b>Preconditions:</b><ul> 00868 * <li>hasGID(cellDim,cellGID)==true</tt> 00869 * </ul> 00870 * 00871 * <b>Postconditions:</b><ul> 00872 * <li>0 <= returnVal < numCells(cellDim)</tt> 00873 * </ul> 00874 */ 00875 virtual int mapGIDToLID(int cellDim, int cellGID) const = 0 ; 00876 00877 /** \brief Find the global ID of a cell given its LID. 00878 * 00879 * \param cellDim 00880 * [in] Dimension of the cell 00881 * \param cellLID 00882 * [in] Local ID of the cell 00883 * 00884 * <b>Postconditions:</b><ul> 00885 * <li><tt>returnVal >= 0 </tt> 00886 * </ul> 00887 */ 00888 virtual int mapLIDToGID(int cellDim, int cellLID) const = 0 ; 00889 00890 /** \brief Return if cells of dimension cellDim have been assigned global 00891 * IDs or not. 00892 * 00893 * \param cellDim 00894 * [in] Dimension of the cell 00895 */ 00896 virtual bool hasIntermediateGIDs(int cellDim) const = 0 ; 00897 00898 /** \brief Assign global IDs to cells of dimension cellDim. 00899 * 00900 * \param cellDim 00901 * [in] Dimension of the cell 00902 * 00903 * <b>Postconditions:</b><ul> 00904 * <li><tt>hasIntermediateGIDs(cellDim)==true</tt> 00905 * </ul> 00906 */ 00907 virtual void assignIntermediateCellGIDs(int cellDim) = 0 ; 00908 00909 //@} 00910 00911 /** \name Reordering */ 00912 //@{ 00913 00914 /** \brief Set the reordering strategy to be used with this mesh. */ 00915 void setReorderer(const CellReorderer& reorderer) 00916 {reorderer_ = reorderer.createInstance(this);} 00917 00918 /** \brief Get the reordering strategy to be used with this mesh. */ 00919 const CellReordererImplemBase* reorderer() const 00920 {return reorderer_.get();} 00921 00922 //@} 00923 00924 00925 00926 /** \name Functions for Mesh with hanging nodes */ 00927 //@{ 00928 /** Function returns true if the mesh allows hanging nodes (by refinement), 00929 * false otherwise */ 00930 virtual bool allowsHangingHodes() const { return false; } 00931 00932 /** Function returns true if the specified element is a "hanging" element 00933 * false otherwise <br> 00934 * @param cellDim [in] should be between 0 , D-1 00935 * @param cellLID [in] the local ID of the element */ 00936 virtual bool isElementHangingNode(int cellDim , int cellLID) const { return false; } 00937 00938 /** Returns the index in the parent maxdim Cell of the refinement tree 00939 * @param maxCellLID [in] the LID of the cell */ 00940 virtual int indexInParent(int maxCellLID) const { return 0; } 00941 00942 /** How many children has a refined element. <br> 00943 * This function provides information of either we have bi or trisection */ 00944 virtual int maxChildren() const { return 0;} 00945 00946 /** Function returns the facets of the parent cell (needed for HN treatment) <br> 00947 * @param childCellLID [in] the LID of the maxdim cell, whos parents facets we want 00948 * @param dimFacets [in] the dimension of the facets which we want to have 00949 * @param facetsLIDs [out] the LID of the parents facets (all) in the defined order 00950 * @param parentCellLIDs [out] the maxdim parent cell LID */ 00951 virtual void returnParentFacets( int childCellLID , int dimFacets , 00952 Array<int> &facetsLIDs , int &parentCellLIDs ) const { } 00953 //@} 00954 00955 00956 00957 /** \name Store special weights in the mesh (for Adaptive Cell Integration) */ 00958 //@{ 00959 00960 /** returns the status of the special weights if they are valid <br> 00961 * These weights are usually computed for one setting of the curve (Adaptive Cell Integration)*/ 00962 virtual bool IsSpecialWeightValid() const {return validWeights_;} 00963 00964 /** specifies if the special weights are valid <br> 00965 * if this is false then usually the special weights have to be recomputed */ 00966 virtual void setSpecialWeightValid(bool val) const { validWeights_ = val;} 00967 00968 /** deletes all special weights so those have to be recreated*/ 00969 virtual void flushSpecialWeights() const; 00970 00971 /** verifies if the specified cell with the given dimension has special weights */ 00972 virtual bool hasSpecialWeight(int dim, int cellLID) const; 00973 00974 /** Sets the special weights */ 00975 virtual void setSpecialWeight(int dim, int cellLID, Array<double>& w) const; 00976 00977 /** Returns the special weights */ 00978 virtual void getSpecialWeight(int dim, int cellLID, Array<double>& w) const; 00979 //@} 00980 00981 00982 /** \name Store the intersection/quadrature points for the curve/surf integral <br> 00983 * for a curve or surf integral we need some quadrature points along the curve in one curve <br> 00984 * These */ 00985 //@{ 00986 00987 /** */ 00988 virtual bool IsCurvePointsValid() const {return curvePoints_Are_Valid_;} 00989 00990 /** */ 00991 virtual void setCurvePointsValid(bool val) const { curvePoints_Are_Valid_ = val;} 00992 00993 /** detletes all the points and its normals which have been stored */ 00994 virtual void flushCurvePoints() const; 00995 00996 /** verifies if the specified maxCell has already precalculated quadrature point for one curve */ 00997 virtual bool hasCurvePoints(int maxCellLID , int curveID) const; 00998 00999 /** Sets the points, curve derivatives and curve normals for one maxCell needed for curve/surf integral*/ 01000 virtual void setCurvePoints(int maxCellLID, int curveID , Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const; 01001 01002 /** Gets the points, curve derivatives and curve normals for one maxCell needed for curve/surf integral*/ 01003 virtual void getCurvePoints(int maxCellLID, int curveID , Array<Point>& points , Array<Point>& derivs , Array<Point>& normals) const; 01004 01005 private: 01006 virtual int mapCurveID_to_Index(int curveID) const; 01007 public: 01008 //@} 01009 01010 /** \name Output */ 01011 //@{ 01012 01013 /** \brief Set whether to stagger output in parallel. Set to true for readable 01014 * output in parallel debugging sessions. 01015 * 01016 * This should be normally, as it causes one synchronization point per 01017 * process. 01018 * 01019 * \todo Get rid of this once we update to use Teuchos::FancyOStream since 01020 * parallel outputting will be done in a readable way automatically! 01021 */ 01022 static bool& staggerOutput() {static bool rtn=false; return rtn;} 01023 01024 //@} 01025 01026 private: 01027 01028 int dim_; 01029 01030 MPIComm comm_; 01031 01032 MeshEntityOrder order_; 01033 01034 RCP<CellReordererImplemBase> reorderer_; 01035 01036 01037 01038 /** flag to indicate if the weights stored are valid */ 01039 mutable bool validWeights_; 01040 01041 /** Object to store the special weights for integration */ 01042 mutable Array < Sundance::Map< int , Array<double> > > specialWeights_; 01043 01044 01045 01046 /** true if the curve did not moved, false if those points are not reusable*/ 01047 mutable bool curvePoints_Are_Valid_; 01048 01049 /** how many curves are participating in curve integrals*/ 01050 mutable int nrCurvesForIntegral_; 01051 01052 /** store intersection informations to calculate the curve integral*/ 01053 mutable Array < Sundance::Map< int , Array<Point> > > curvePoints_; 01054 01055 /** store directional derivative informations to calculate the curve integral*/ 01056 mutable Array < Sundance::Map< int , Array<Point> > > curveDerivative_; 01057 01058 /** store normal directional used in the curve or in the surf integral <br> 01059 * in case of the surf integral it is the cross product from the integral */ 01060 mutable Array < Sundance::Map< int , Array<Point> > > curveNormal_; 01061 01062 /** map the curve ID to index in the previous arrays */ 01063 mutable Sundance::Map< int , int > curveID_to_ArrayIndex_; 01064 }; 01065 01066 01067 } // namespace Sundance 01068 01069 #endif