primitiveMeshFaceCentresAndAreas.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-2020 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  Calculate the face centres and areas.
26 
27  Calculate the centre by breaking the face into triangles using the face
28  centre and area-weighted averaging their centres. This method copes with
29  small face-concavity.
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "primitiveMesh.H"
34 
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
39 {
40  if (debug)
41  {
42  Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
43  << "Calculating face centres and face areas"
44  << endl;
45  }
46 
47  // It is an error to attempt to recalculate faceCentres
48  // if the pointer is already set
49  if (faceCentresPtr_ || faceAreasPtr_ || magFaceAreasPtr_)
50  {
52  << "Face centres or face areas already calculated"
53  << abort(FatalError);
54  }
55 
56  faceCentresPtr_ = new vectorField(nFaces());
57  vectorField& fCtrs = *faceCentresPtr_;
58 
59  faceAreasPtr_ = new vectorField(nFaces());
60  vectorField& fAreas = *faceAreasPtr_;
61 
62  magFaceAreasPtr_ = new scalarField(nFaces());
63  scalarField& magfAreas = *magFaceAreasPtr_;
64 
65  makeFaceCentresAndAreas(points(), fCtrs, fAreas, magfAreas);
66 
67  if (debug)
68  {
69  Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
70  << "Finished calculating face centres and face areas"
71  << endl;
72  }
73 }
74 
75 
77 (
78  const pointField& p,
79  vectorField& fCtrs,
80  vectorField& fAreas,
81  scalarField& magfAreas
82 ) const
83 {
84  const faceList& fs = faces();
85 
86  forAll(fs, facei)
87  {
88  const labelList& f = fs[facei];
89  label nPoints = f.size();
90 
91  // If the face is a triangle, do a direct calculation for efficiency
92  // and to avoid round-off error-related problems
93  if (nPoints == 3)
94  {
95  fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
96  fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
97  }
98 
99  // For more complex faces, decompose into triangles
100  else
101  {
102  // Compute an estimate of the centre as the average of the points
103  point pAvg = p[f[0]];
104  for (label pi = 1; pi < nPoints; pi++)
105  {
106  pAvg += p[f[pi]];
107  }
108  pAvg /= nPoints;
109 
110  // Compute the face area normal and unit normal by summing up the
111  // normals of the triangles formed by connecting each edge to the
112  // point average.
113  vector sumA = Zero;
114  forAll(f, i)
115  {
116  const vector a =
117  (p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
118 
119  sumA += a;
120  }
121  const vector sumAHat = normalised(sumA);
122 
123  // Compute the area-weighted sum of the triangle centres. Note use
124  // the triangle area projected in the direction of the face normal
125  // as the weight, *not* the triangle area magnitude. Only the
126  // former makes the calculation independent of the initial estimate.
127  scalar sumAn = 0.0;
128  vector sumAnc = Zero;
129  forAll(f, i)
130  {
131  const vector a =
132  (p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
133  const vector c = p[f[i]] + p[f[f.fcIndex(i)]] + pAvg;
134 
135  const scalar an = a & sumAHat;
136 
137  sumAn += an;
138  sumAnc += an*c;
139  }
140 
141  // Complete calculating centres and areas. If the face is too small
142  // for the sums to be reliably divided then just set the centre to
143  // the initial estimate.
144  if (sumAn > vSmall)
145  {
146  fCtrs[facei] = (1.0/3.0)*sumAnc/sumAn;
147  }
148  else
149  {
150  fCtrs[facei] = pAvg;
151  }
152  fAreas[facei] = 0.5*sumA;
153  }
154 
155  magfAreas[facei] = max(mag(fAreas[facei]), vSmall);
156  }
157 }
158 
159 
160 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
161 
163 {
164  if (!faceCentresPtr_)
165  {
167  }
168 
169  return *faceCentresPtr_;
170 }
171 
172 
174 {
175  if (!faceAreasPtr_)
176  {
178  }
179 
180  return *faceAreasPtr_;
181 }
182 
183 
185 {
186  if (!magFaceAreasPtr_)
187  {
189  }
190 
191  return *magFaceAreasPtr_;
192 }
193 
194 
195 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void makeFaceCentresAndAreas(const pointField &p, vectorField &fCtrs, vectorField &fAreas, scalarField &magfAreas) const
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
label nFaces() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const pointField & points() const =0
Return mesh points.
const dimensionedScalar c
Speed of light in a vacuum.
const scalarField & magFaceAreas() const
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
labelList f(nPoints)
const vectorField & faceCentres() const
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
void calcFaceCentresAndAreas() const
Calculate face centres and areas.
const vectorField & faceAreas() const
virtual const faceList & faces() const =0
Return faces.
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
Field< vector > vectorField
Specialisation of Field<T> for vector.