src/foam/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C

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