src/tetDecompositionFiniteElement/tetPolyMeshFaceDecomp/tetPolyMeshFaceDecomp.C

   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 Description
  26     The tetrahedral decomposition of the mesh based on the underlying polyMesh.
  27     Uses minimal storage - the tet decomposition is provided on a cell-by-cell
  28     basis
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "tetPolyMeshFaceDecomp.H"
  33 #include "demandDrivenData.H"
  34 #include "tetPolyPatchFields.H"
  35 #include "tetPointFields.H"
  36 #include "elementFields.H"
  37 #include "tetFemMatrices.H"
  38 #include "mapPolyMesh.H"
  39 #include "MapTetFemFields.H"
  40 #include "tetPolyMeshMapperFaceDecomp.H"
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 defineTypeNameAndDebug(Foam::tetPolyMeshFaceDecomp, 0);
  45
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49 namespace Foam
  50 {
  51
  52 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  53
  54 void tetPolyMeshFaceDecomp::clearOut() const
  55 {
  56     nPoints_ = -1;
  57     nEdges_ = -1;
  58     nTets_ = -1;
  59     deleteDemandDrivenData(lduPtr_);
  60     maxNPointsForCell_ = -1;
  61
  62     clearOutParPointData();
  63 }
  64
  65 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  66
  67 // Construct from components
  68 tetPolyMeshFaceDecomp::tetPolyMeshFaceDecomp(const polyMesh& pMesh)
  69 :
  70     GeoMesh<polyMesh>(pMesh),
  71     tetFemSolution(pMesh),
  72     boundary_(*this, pMesh.boundaryMesh()),
  73     faceOffset_(mesh_.nPoints()),
  74     cellOffset_(faceOffset_ + mesh_.nFaces()),
  75     nPoints_(-1),
  76     nEdges_(-1),
  77     nTets_(-1),
  78     lduPtr_(NULL),
  79     maxNPointsForCell_(-1),
  80     parPointsPtr_(NULL),
  81     parEdgesPtr_(NULL)
  82 {
  83     if (debug)
  84     {
  85         Info<< "tetPolyMesh::tetPolyMesh(const polyMesh&) : "
  86             << "Creating tetPolyMesh" << endl;
  87     }
  88
  89     addParallelPointPatch();
  90 }
  91
  92
  93 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  94
  95 tetPolyMeshFaceDecomp::~tetPolyMeshFaceDecomp()
  96 {
  97     if (debug)
  98     {
  99         Info<< "tetPolyMesh::~tetPolyMesh() : "
 100             << "Deleting tetPolyMesh" << endl;
 101     }
 102
 103     clearOut();
 104 }
 105
 106
 107 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 108
 109 // Return number of edges in decomposition for a face
 110 label tetPolyMeshFaceDecomp::nEdgesForFace(const label faceID) const
 111 {
 112     // It is all the edges around the circumference + the internal diags
 113
 114     // Note: compiler bug.  It should be possible to do this without wrapping
 115     const faceList& f = mesh_.faces();
 116     return 2*f[faceID].size();
 117 }
 118
 119
 120 // Return number of edges in decomposition connected to a given point
 121 label tetPolyMeshFaceDecomp::nEdgesForPoint(const label pointID) const
 122 {
 123     const label startFaceOwn = lduAddr().ownerStartAddr()[pointID];
 124     const label endFaceOwn = lduAddr().ownerStartAddr()[pointID + 1];
 125
 126     const label startFaceNbr = lduAddr().losortStartAddr()[pointID];
 127     const label endFaceNbr = lduAddr().losortStartAddr()[pointID + 1];
 128
 129     // pointID is the owner of the first lot and the neighbour of the second lot
 130     return (endFaceOwn - startFaceOwn) + (endFaceNbr - startFaceNbr);
 131 }
 132
 133
 134 // Return number of edges in decomposition connected to a given point
 135 labelList tetPolyMeshFaceDecomp::edgesForPoint(const label pointID) const
 136 {
 137     const label startFaceOwn = lduAddr().ownerStartAddr()[pointID];
 138     const label endFaceOwn = lduAddr().ownerStartAddr()[pointID + 1];
 139
 140     const label startFaceNbr = lduAddr().losortStartAddr()[pointID];
 141     const label endFaceNbr = lduAddr().losortStartAddr()[pointID + 1];
 142
 143     const unallocLabelList& losort = lduAddr().losortAddr();
 144
 145     // pointID is the owner of the first lot and the neighbour of the second lot
 146     labelList edgeIndices(nEdgesForPoint(pointID), -1);
 147
 148     label i = 0;
 149
 150     // owner side
 151     for (label edgeLabel = startFaceOwn; edgeLabel < endFaceOwn; edgeLabel++)
 152     {
 153         edgeIndices[i] = edgeLabel;
 154         i++;
 155     }
 156
 157     // neighbour side
 158     for (label edgeLabel = startFaceNbr; edgeLabel < endFaceNbr; edgeLabel++)
 159     {
 160         edgeIndices[i] = losort[edgeLabel];
 161         i++;
 162     }
 163
 164     return edgeIndices;
 165 }
 166
 167
 168 void tetPolyMeshFaceDecomp::updateMesh
 169 (
 170     const tetPolyMeshMapperFaceDecomp& mapper
 171 )
 172 {
 173     // It is assumed that polyMesh::morph() has already been called.
 174     if (debug)
 175     {
 176         Info<< "void tetPolyMeshFaceDecomp::updateMesh: "
 177             << "Mesh update on topological change" << endl;
 178     }
 179
 180     // Clear out all data
 181     clearOut();
 182
 183     // Reset face and cell offset (topology has changed)
 184     faceOffset_ = mesh_.nPoints();
 185     cellOffset_ = faceOffset_ + mesh_.nFaces();
 186
 187     boundary_.updateMesh();
 188
 189     // Map all the tetPointFields in the database
 190
 191     MapGeometricFields
 192     <
 193         scalar,
 194         tetPolyPatchField,
 195         tetPolyMeshMapperFaceDecomp,
 196         tetPointMesh
 197     >(mapper);
 198
 199     MapGeometricFields
 200     <
 201         vector,
 202         tetPolyPatchField,
 203         tetPolyMeshMapperFaceDecomp,
 204         tetPointMesh
 205     >(mapper);
 206
 207     MapGeometricFields
 208     <
 209         tensor,
 210         tetPolyPatchField,
 211         tetPolyMeshMapperFaceDecomp,
 212         tetPointMesh
 213     >(mapper);
 214
 215     // Map all the element fields in the database
 216
 217     MapGeometricFields
 218     <
 219         scalar,
 220         elementPatchField,
 221         tetPolyMeshMapperFaceDecomp,
 222         elementMesh
 223     >(mapper);
 224
 225     MapGeometricFields
 226     <
 227         vector,
 228         elementPatchField,
 229         tetPolyMeshMapperFaceDecomp,
 230         elementMesh
 231     >(mapper);
 232
 233     MapGeometricFields
 234     <
 235         tensor,
 236         elementPatchField,
 237         tetPolyMeshMapperFaceDecomp,
 238         elementMesh
 239     >(mapper);
 240 }
 241
 242
 243 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 244
 245 bool tetPolyMeshFaceDecomp::operator!=(const tetPolyMeshFaceDecomp& bm) const
 246 {
 247     return &bm != this;
 248 }
 249
 250
 251 bool tetPolyMeshFaceDecomp::operator==(const tetPolyMeshFaceDecomp& bm) const
 252 {
 253     return &bm == this;
 254 }
 255
 256
 257 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 258
 259 } // End namespace Foam
 260
 261 // ************************************************************************* //