primitiveMeshCellCentresAndVols.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 
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  {
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 
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 = 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(), 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  fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]);
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  fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]);
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  forAll(cellCtrs, celli)
143  {
144  if (mag(cellVols[celli]) > vSmall)
145  {
146  cellCtrs[celli] /= cellVols[celli];
147  }
148  else
149  {
150  cellCtrs[celli] = cEst[celli];
151  }
152  }
153 
154  cellVols *= (1.0/3.0);
155 }
156 
157 
158 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
159 
161 {
162  if (!cellCentresPtr_)
163  {
165  }
166 
167  return *cellCentresPtr_;
168 }
169 
170 
172 {
173  if (!cellVolumesPtr_)
174  {
176  }
177 
178  return *cellVolumesPtr_;
179 }
180 
181 
182 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const labelList & faceNeighbour() const =0
Face face-neighbour addressing.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void calcCellCentresAndVols() const
Calculate cell centres and volumes.
static const zero Zero
Definition: zero.H:97
void makeCellCentresAndVols(const vectorField &fCtrs, const vectorField &fAreas, vectorField &cellCtrs, scalarField &cellVols) const
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const vectorField & faceCentres() const
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const vectorField & faceAreas() const
virtual const labelList & faceOwner() const =0
Face face-owner addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const scalarField & cellVolumes() const