primitiveMeshFaceCentresAndAreas.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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  Calulate 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_)
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  makeFaceCentresAndAreas(points(), fCtrs, fAreas);
63 
64  if (debug)
65  {
66  Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
67  << "Finished calculating face centres and face areas"
68  << endl;
69  }
70 }
71 
72 
74 (
75  const pointField& p,
76  vectorField& fCtrs,
77  vectorField& fAreas
78 ) const
79 {
80  const faceList& fs = faces();
81 
82  forAll(fs, facei)
83  {
84  const labelList& f = fs[facei];
85  label nPoints = f.size();
86 
87  // If the face is a triangle, do a direct calculation for efficiency
88  // and to avoid round-off error-related problems
89  if (nPoints == 3)
90  {
91  fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
92  fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
93  }
94  else
95  {
96  vector sumN = Zero;
97  scalar sumA = 0.0;
98  vector sumAc = Zero;
99 
100  point fCentre = p[f[0]];
101  for (label pi = 1; pi < nPoints; pi++)
102  {
103  fCentre += p[f[pi]];
104  }
105 
106  fCentre /= nPoints;
107 
108  for (label pi = 0; pi < nPoints; pi++)
109  {
110  const point& nextPoint = p[f[(pi + 1) % nPoints]];
111 
112  vector c = p[f[pi]] + nextPoint + fCentre;
113  vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
114  scalar a = mag(n);
115 
116  sumN += n;
117  sumA += a;
118  sumAc += a*c;
119  }
120 
121  // This is to deal with zero-area faces. Mark very small faces
122  // to be detected in e.g., processorPolyPatch.
123  if (sumA < ROOTVSMALL)
124  {
125  fCtrs[facei] = fCentre;
126  fAreas[facei] = Zero;
127  }
128  else
129  {
130  fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
131  fAreas[facei] = 0.5*sumN;
132  }
133  }
134  }
135 }
136 
137 
138 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
139 
141 {
142  if (!faceCentresPtr_)
143  {
145  }
146 
147  return *faceCentresPtr_;
148 }
149 
150 
152 {
153  if (!faceAreasPtr_)
154  {
156  }
157 
158  return *faceAreasPtr_;
159 }
160 
161 
162 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label nFaces() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual const pointField & points() const =0
Return mesh points.
void makeFaceCentresAndAreas(const pointField &p, vectorField &fCtrs, vectorField &fAreas) const
static const zero Zero
Definition: zero.H:91
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.
const dimensionedScalar c
Speed of light in a vacuum.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
label nPoints() const
Field< vector > vectorField
Specialisation of Field<T> for vector.