src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/tetPolyMeshCellDecomp.H

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