src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.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     Efficient cell-centre calculation using face-addressing, face-centres and
  27     face-areas.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "primitiveMesh.H"
  32
  33 #include "tetPointRef.H"
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 void Foam::primitiveMesh::calcCellCentresAndVols() const
  38 {
  39     if (debug)
  40     {
  41         Pout<< "primitiveMesh::calcCellCentresAndVols() : "
  42             << "Calculating cell centres and cell volumes"
  43             << endl;
  44     }
  45
  46     // It is an error to attempt to recalculate cellCentres
  47     // if the pointer is already set
  48     if (cellCentresPtr_ || cellVolumesPtr_)
  49     {
  50         FatalErrorIn("primitiveMesh::calcCellCentresAndVols() const")
  51             << "Cell centres or cell volumes already calculated"
  52             << abort(FatalError);
  53     }
  54
  55     // set the accumulated cell centre to zero vector
  56     cellCentresPtr_ = new vectorField(nCells());
  57     vectorField& cellCtrs = *cellCentresPtr_;
  58
  59     // Initialise cell volumes to 0
  60     cellVolumesPtr_ = new scalarField(nCells());
  61     scalarField& cellVols = *cellVolumesPtr_;
  62
  63     // Make centres and volumes
  64     makeCellCentresAndVols(faceCentres(), faceAreas(), cellCtrs, cellVols);
  65
  66     if (debug)
  67     {
  68         Pout<< "primitiveMesh::calcCellCentresAndVols() : "
  69             << "Finished calculating cell centres and cell volumes"
  70             << endl;
  71     }
  72 }
  73
  74
  75 void Foam::primitiveMesh::makeCellCentresAndVols
  76 (
  77     const vectorField& fCtrs,
  78     const vectorField& fAreas,
  79     vectorField& cellCtrs,
  80     scalarField& cellVols
  81 ) const
  82 {
  83     // Clear the fields for accumulation
  84     cellCtrs = vector::zero;
  85     cellVols = 0.0;
  86
  87     const labelList& own = faceOwner();
  88     const labelList& nei = faceNeighbour();
  89
  90     // first estimate the approximate cell centre as the average of
  91     // face centres
  92
  93     vectorField cEst(nCells(), vector::zero);
  94     labelField nCellFaces(nCells(), 0);
  95
  96     forAll (own, facei)
  97     {
  98         cEst[own[facei]] += fCtrs[facei];
  99         nCellFaces[own[facei]] += 1;
 100     }
 101
 102     forAll (nei, facei)
 103     {
 104         cEst[nei[facei]] += fCtrs[facei];
 105         nCellFaces[nei[facei]] += 1;
 106     }
 107
 108     forAll(cEst, celli)
 109     {
 110         cEst[celli] /= nCellFaces[celli];
 111     }
 112
 113     const faceList& allFaces = faces();
 114     const pointField& allPoints = points();
 115
 116     forAll(own, faceI)
 117     {
 118         const face& f = allFaces[faceI];
 119
 120         if (f.size() == 3)
 121         {
 122             tetPointRef tpr
 123             (
 124                 allPoints[f[2]],
 125                 allPoints[f[1]],
 126                 allPoints[f[0]],
 127                 cEst[own[faceI]]
 128             );
 129
 130             scalar tetVol = tpr.mag();
 131
 132             // Accumulate volume-weighted tet centre
 133             cellCtrs[own[faceI]] += tetVol*tpr.centre();
 134
 135             // Accumulate tet volume
 136             cellVols[own[faceI]] += tetVol;
 137         }
 138         else
 139         {
 140             forAll(f, pI)
 141             {
 142                 tetPointRef tpr
 143                 (
 144                     allPoints[f[pI]],
 145                     allPoints[f.prevLabel(pI)],
 146                     fCtrs[faceI],
 147                     cEst[own[faceI]]
 148                 );
 149
 150                 scalar tetVol = tpr.mag();
 151
 152                 // Accumulate volume-weighted tet centre
 153                 cellCtrs[own[faceI]] += tetVol*tpr.centre();
 154
 155                 // Accumulate tet volume
 156                 cellVols[own[faceI]] += tetVol;
 157             }
 158         }
 159     }
 160
 161     forAll(nei, faceI)
 162     {
 163         const face& f = allFaces[faceI];
 164
 165         if (f.size() == 3)
 166         {
 167             tetPointRef tpr
 168             (
 169                 allPoints[f[0]],
 170                 allPoints[f[1]],
 171                 allPoints[f[2]],
 172                 cEst[nei[faceI]]
 173             );
 174
 175             scalar tetVol = tpr.mag();
 176
 177             // Accumulate volume-weighted tet centre
 178             cellCtrs[nei[faceI]] += tetVol*tpr.centre();
 179
 180             // Accumulate tet volume
 181             cellVols[nei[faceI]] += tetVol;
 182         }
 183         else
 184         {
 185             forAll(f, pI)
 186             {
 187                 tetPointRef tpr
 188                 (
 189                     allPoints[f[pI]],
 190                     allPoints[f.nextLabel(pI)],
 191                     fCtrs[faceI],
 192                     cEst[nei[faceI]]
 193                 );
 194
 195                 scalar tetVol = tpr.mag();
 196
 197                 // Accumulate volume-weighted tet centre
 198                 cellCtrs[nei[faceI]] += tetVol*tpr.centre();
 199
 200                 // Accumulate tet volume
 201                 cellVols[nei[faceI]] += tetVol;
 202             }
 203         }
 204     }
 205
 206     cellCtrs /= cellVols + VSMALL;
 207 }
 208
 209
 210 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 211
 212 const Foam::vectorField& Foam::primitiveMesh::cellCentres() const
 213 {
 214     if (!cellCentresPtr_)
 215     {
 216         calcCellCentresAndVols();
 217     }
 218
 219     return *cellCentresPtr_;
 220 }
 221
 222
 223 const Foam::scalarField& Foam::primitiveMesh::cellVolumes() const
 224 {
 225     if (!cellVolumesPtr_)
 226     {
 227         calcCellCentresAndVols();
 228     }
 229
 230     return *cellVolumesPtr_;
 231 }
 232
 233
 234 // ************************************************************************* //