src/tetFiniteElement/tetPolyMesh/tetPolyMesh.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Class
  25     tetPolyMesh
  26
  27 Description
  28     The tetrahedral decomposition of the mesh based on the underlying polyMesh.
  29     Uses minimal storage - the tet decomposition is provided on a cell-by-cell
  30     basis.
  31
  32     This class will be used as the underlying structure for FEM mesh motion
  33     and therefore it needs to be able to create and store tet point-edge
  34     addressing for the solver.
  35
  36     Note on addressing:
  37         The problem of addressing for Finite elements is two-fold. The
  38         first issue is addressing in the FEM matrix. The FEM coefficients are
  39         associated with edges in the tet decomposition. In order for the FOAM
  40         solver to work correctly, the coefficients need to be ordered in the
  41         upper triangular order. It is convinient to keep the order of edges
  42         and off-diagonal coefficients the same; therefore, the order of edges
  43         in the tet mesh is fixed by the ordering of points. The construction
  44         of the list from the edge list in polyMesh will be described in
  45         calcPolyMeshAddressing.C.
  46
  47         The second issue is the construction of the FEM matrix.
  48         A fast method (Wavre 19/Sep/2003) is currently used:
  49
  50         1) take a group of elements.  The discretisation will be done
  51         over a group.  A goup should be such that I can allocate a
  52         dense square matrix for all elements in the group.  Do NOT
  53         initialise arrays to zero.  As only a part of the group will
  54         be addressed, it needs to be zeroed out only for the entries
  55         it will have to use.  As the number of vertices per group will
  56         be variable, the local dense matrix changes for every group.
  57
  58         2).  Allocate group-to-global and global-to-group arrays.  The
  59         second one is full of -1-s for vertices that are not in the
  60         group.  It will be initialised to -1 and is kept outside of
  61         the group.  This is done by going through all the vertices in
  62         the group and looking up a global-to-local index. If the index
  63         is -1, this global index has not been seen before: allocate
  64         next free local index to the vertex.
  65
  66         3) For every element in the local-to-global, fidn the matrix
  67         row in ldu addressing.  For each element in a row, using
  68         global-to-local, find the local off-diagonal index and zero it
  69         out.  These are the only ones being used!
  70
  71         4) Fill in local matrix based on the local vertex indices into
  72         the local dense matrix.
  73
  74         5) Assemble global matrix by doing the following:
  75             - go through all group-to-global indices
  76             - for each group-to-global grab the row
  77             - for each non-zero element of the row, look up
  78             global-to-group.  This gives the x,y coordinates in the
  79             dense matrix and and insert it.
  80
  81         6) Clear the global-to-local addressing by looking up the
  82         local-to-global indices and zeroing out only those entries.
  83
  84 SourceFiles
  85     tetPolyMesh.C
  86     calcTetPolyMeshGeometry.C
  87     calcTetPolyMeshAddressing.C
  88     calcTetPolyMeshParPointData.C
  89     addParallelPointPatch.C
  90
  91 \*---------------------------------------------------------------------------*/
  92
  93 #ifndef tetPolyMesh_H
  94 #define tetPolyMesh_H
  95
  96 #include "GeoMesh.H"
  97 #include "polyMesh.H"
  98 #include "lduMesh.H"
  99 #include "tetFemSolution.H"
 100 #include "primitiveMesh.H"
 101 #include "tetPolyBoundaryMesh.H"
 102 #include "tetCellList.H"
 103 #include "tetPolyMeshLduAddressing.H"
 104 #include "cellShapeList.H"
 105 #include "FieldFields.H"
 106
 107 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 108
 109 namespace Foam
 110 {
 111
 112 // Class forward declarations
 113 template<class Type> class SquareMatrix;
 114
 115 class tetPolyMeshMapper;
 116 class tetPolyMeshLduAddressing;
 117
 118 /*---------------------------------------------------------------------------*\
 119                      Class tetPolyMesh Declaration
 120 \*---------------------------------------------------------------------------*/
 121
 122 class tetPolyMesh
 123 :
 124     public GeoMesh<polyMesh>,
 125     public lduMesh,
 126     public tetFemSolution
 127 {
 128     // Permanent data
 129
 130         // Minimum mesh data
 131
 132         //- Boundary mesh
 133         tetPolyBoundaryMesh boundary_;
 134
 135         //- Offset in numbering to first face centre
 136         label faceOffset_;
 137
 138         //- Offset in numbering to first cell centre
 139         label cellOffset_;
 140
 141
 142     // Demand-driven data
 143
 144         //- Number of points
 145         mutable label nPoints_;
 146
 147         //- Number of edges
 148         mutable label nEdges_;
 149
 150         //- Number of tetrahedra
 151         mutable label nTets_;
 152
 153         //- LDU addressing
 154         mutable tetPolyMeshLduAddressing* lduPtr_;
 155
 156         //- Max number of points per cell
 157         mutable label maxNPointsForCell_;
 158
 159
 160         // Parallel mesh analysis tools
 161
 162             //- List of points shared between processor patches
 163             mutable labelList* parPointsPtr_;
 164
 165             //- List of edges shared between processor patches
 166             mutable edgeList* parEdgesPtr_;
 167
 168
 169     // Private Member Functions
 170
 171         //- Disallow default bitwise copy construct
 172         tetPolyMesh(const tetPolyMesh&);
 173
 174         //- Disallow default bitwise assignment
 175         void operator=(const tetPolyMesh&);
 176
 177
 178     // Private member functions to calculate demand driven data
 179
 180         //- Calculate addressing
 181         void calcAddressing() const;
 182
 183         //- Check for and create a parallel point patch
 184         void addParallelPointPatch();
 185
 186         //- Calculate points shared between parallel patches
 187         void calcParPointData() const;
 188         void clearOutParPointData() const;
 189
 190         //- Clear out all data
 191         void clearOut() const;
 192
 193
 194 public:
 195
 196     // Declare name of the class and its debug switch
 197     ClassName("tetPolyMesh");
 198
 199     typedef tetPolyMesh Mesh;
 200     typedef tetPolyBoundaryMesh BoundaryMesh;
 201
 202
 203     // Constructors
 204
 205         //- Construct from components
 206         explicit tetPolyMesh(const polyMesh& pMesh);
 207
 208     // Destructor
 209
 210         ~tetPolyMesh();
 211
 212
 213     // Member Functions
 214
 215         // Access
 216
 217             //- Return the top-level database
 218             const Time& time() const
 219             {
 220                 return mesh_.time();
 221             }
 222
 223             //- Return the object registry
 224             virtual const objectRegistry& thisDb() const
 225             {
 226                 return GeoMesh<polyMesh>::thisDb();
 227             }
 228
 229             //- Return number of points in decomposition
 230             label nPoints() const;
 231
 232             //- Return number of edges in decomposition
 233             label nEdges() const;
 234
 235             //- Return number of tetrahedra in decomposition
 236             label nTets() const;
 237
 238             //- Return number of cells in polyhedral mesh
 239             label nCells() const
 240             {
 241                 return mesh_.nCells();
 242             }
 243
 244             //- Return number of tetrahedra in decomposition for cell
 245             label nTetsForCell(const label cellID) const;
 246
 247             //- Return number of edges in decomposition for a face
 248             label nEdgesForFace(const label faceID) const;
 249
 250             //- Return number of edges in decomposition connected to a
 251             //  given point
 252             label nEdgesForPoint(const label pointID) const;
 253
 254             //- Return face offset
 255             label faceOffset() const
 256             {
 257                 return faceOffset_;
 258             }
 259
 260             //- Return cell offset
 261             label cellOffset() const
 262             {
 263                 return cellOffset_;
 264             }
 265
 266             //- Return list of edge labels coming out of a point
 267             labelList edgesForPoint(const label pointID) const;
 268
 269             //- Return points
 270             tmp<pointField> points() const;
 271
 272             //- Return complete list of cell shapes. All are tetrahedra
 273             cellShapeList tetCells() const;
 274
 275             //- Return reference to boundary mesh
 276             const tetPolyBoundaryMesh& boundary() const
 277             {
 278                 return boundary_;
 279             }
 280
 281             //- Return ldu addressing
 282             const lduAddressing& lduAddr() const;
 283
 284             //- Return a list of pointers for each patch
 285             //  with only those pointing to interfaces being set
 286             virtual lduInterfacePtrsList interfaces() const
 287             {
 288                 return boundary().interfaces();
 289             }
 290
 291             //- Owner
 292             const unallocLabelList& owner() const
 293             {
 294                 return lduAddr().lowerAddr();
 295             }
 296
 297             //- Neighbour
 298             const unallocLabelList& neighbour() const
 299             {
 300                 return lduAddr().upperAddr();
 301             }
 302
 303             //- Return tetrahedral decomposition for cell
 304             tetCellList tets(const label cellID) const;
 305
 306             //- Return max number of tets in a cell
 307             label maxNPointsForCell() const;
 308
 309             //- Fill buffer with tet decomposition addressing
 310             //  Used for FEM matrix assembly.
 311             //  localToGlobalBuffer - sized to max number of vertices per cell
 312             //                        in the mesh
 313             //  globalToLocalBuffer - sized to total number of points in
 314             //                        the mesh and initialised to -1
 315             label addressing
 316             (
 317                 const label cellID,
 318                 labelList& localToGlobalBuffer,
 319                 labelList& globalToLocalBuffer
 320             ) const;
 321
 322             //- Clear global to local addressing
 323             void clearAddressing
 324             (
 325                 const label cellID,
 326                 const label nCellPoints,
 327                 labelList& localToGlobalBuffer,
 328                 labelList& globalToLocalBuffer
 329             ) const;
 330
 331             //- Fill buffer with dot-products of shape functions
 332             // Used for FEM matrix assembly
 333             void gradNiDotGradNj
 334             (
 335                 const label cellID,
 336                 SquareMatrix<scalar>& denseMatrix,
 337                 const labelList& globalToLocalBuffer
 338             ) const;
 339
 340             //- Fill buffer with tensor products of shape functions
 341             // Used for FEM matrix assembly
 342             void gradNiGradNj
 343             (
 344                 const label cellID,
 345                 SquareMatrix<tensor>& denseMatrix,
 346                 const labelList& globalToLocalBuffer
 347             ) const;
 348
 349             //- Fill buffer with the volume integral distributed into vertices
 350             void volIntegral
 351             (
 352                 const label cellID,
 353                 scalarField& buffer,
 354                 const labelList& globalToLocalBuffer
 355             ) const;
 356
 357
 358             // Parallel mesh analysis data
 359
 360                 //- Return parallel info
 361                 const globalMeshData& globalData() const
 362                 {
 363                     return mesh_.globalData();
 364                 }
 365
 366                 //- Shared parallel points
 367                 const labelList& parallelPoints() const;
 368
 369                 //- Shared parallel edges
 370                 const edgeList& parallelEdges() const;
 371
 372
 373         // Edit
 374
 375             //- Update mesh topology using the morph engine
 376             void updateMesh(const tetPolyMeshMapper& mapper);
 377
 378
 379     // Member Operators
 380
 381         bool operator!=(const tetPolyMesh&) const;
 382         bool operator==(const tetPolyMesh&) const;
 383 };
 384
 385
 386 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 387
 388 } // End namespace Foam
 389
 390 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 391
 392 #endif
 393
 394 // ************************************************************************* //