src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 void Foam::primitiveMesh::calcCellCentresAndVols() const
  35 {
  36     if (debug)
  37     {
  38         Pout<< "primitiveMesh::calcCellCentresAndVols() : "
  39             << "Calculating cell centres and cell volumes"
  40             << endl;
  41     }
  42
  43     // It is an error to attempt to recalculate cellCentres
  44     // if the pointer is already set
  45     if (cellCentresPtr_ || cellVolumesPtr_)
  46     {
  47         FatalErrorIn("primitiveMesh::calcCellCentresAndVols() const")
  48             << "Cell centres or cell volumes already calculated"
  49             << abort(FatalError);
  50     }
  51
  52     // set the accumulated cell centre to zero vector
  53     cellCentresPtr_ = new vectorField(nCells());
  54     vectorField& cellCtrs = *cellCentresPtr_;
  55
  56     // Initialise cell volumes to 0
  57     cellVolumesPtr_ = new scalarField(nCells());
  58     scalarField& cellVols = *cellVolumesPtr_;
  59
  60     // Make centres and volumes
  61     makeCellCentresAndVols(faceCentres(), faceAreas(), cellCtrs, cellVols);
  62
  63     if (debug)
  64     {
  65         Pout<< "primitiveMesh::calcCellCentresAndVols() : "
  66             << "Finished calculating cell centres and cell volumes"
  67             << endl;
  68     }
  69 }
  70
  71
  72 void Foam::primitiveMesh::makeCellCentresAndVols
  73 (
  74     const vectorField& fCtrs,
  75     const vectorField& fAreas,
  76     vectorField& cellCtrs,
  77     scalarField& cellVols
  78 ) const
  79 {
  80     // Clear the fields for accumulation
  81     cellCtrs = vector::zero;
  82     cellVols = 0.0;
  83
  84     const labelList& own = faceOwner();
  85     const labelList& nei = faceNeighbour();
  86
  87     // first estimate the approximate cell centre as the average of
  88     // face centres
  89
  90     vectorField cEst(nCells(), vector::zero);
  91     labelField nCellFaces(nCells(), 0);
  92
  93     forAll(own, facei)
  94     {
  95         cEst[own[facei]] += fCtrs[facei];
  96         nCellFaces[own[facei]] += 1;
  97     }
  98
  99     forAll(nei, facei)
 100     {
 101         cEst[nei[facei]] += fCtrs[facei];
 102         nCellFaces[nei[facei]] += 1;
 103     }
 104
 105     forAll(cEst, celli)
 106     {
 107         cEst[celli] /= nCellFaces[celli];
 108     }
 109
 110     forAll(own, facei)
 111     {
 112         // Calculate 3*face-pyramid volume
 113         scalar pyr3Vol =
 114             max(fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]), VSMALL);
 115
 116         // Calculate face-pyramid centre
 117         vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[own[facei]];
 118
 119         // Accumulate volume-weighted face-pyramid centre
 120         cellCtrs[own[facei]] += pyr3Vol*pc;
 121
 122         // Accumulate face-pyramid volume
 123         cellVols[own[facei]] += pyr3Vol;
 124     }
 125
 126     forAll(nei, facei)
 127     {
 128         // Calculate 3*face-pyramid volume
 129         scalar pyr3Vol =
 130             max(fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]), VSMALL);
 131
 132         // Calculate face-pyramid centre
 133         vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[nei[facei]];
 134
 135         // Accumulate volume-weighted face-pyramid centre
 136         cellCtrs[nei[facei]] += pyr3Vol*pc;
 137
 138         // Accumulate face-pyramid volume
 139         cellVols[nei[facei]] += pyr3Vol;
 140     }
 141
 142     cellCtrs /= cellVols;
 143     cellVols *= (1.0/3.0);
 144 }
 145
 146
 147 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 148
 149 const Foam::vectorField& Foam::primitiveMesh::cellCentres() const
 150 {
 151     if (!cellCentresPtr_)
 152     {
 153         calcCellCentresAndVols();
 154     }
 155
 156     return *cellCentresPtr_;
 157 }
 158
 159
 160 const Foam::scalarField& Foam::primitiveMesh::cellVolumes() const
 161 {
 162     if (!cellVolumesPtr_)
 163     {
 164         calcCellCentresAndVols();
 165     }
 166
 167     return *cellVolumesPtr_;
 168 }
 169
 170
 171 // ************************************************************************* //