FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 
38 void Foam::primitiveMesh::calcFaceCentresAndAreas() const
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  {
51  FatalErrorIn("primitiveMesh::calcFaceCentresAndAreas() const")
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 
73 void Foam::primitiveMesh::makeFaceCentresAndAreas
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 = vector::zero;
97  scalar sumA = 0.0;
98  vector sumAc = vector::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  fCtrs[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
122  fAreas[facei] = 0.5*sumN;
123  }
124  }
125 }
126 
127 
128 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
129 
131 {
132  if (!faceCentresPtr_)
133  {
134  calcFaceCentresAndAreas();
135  }
136 
137  return *faceCentresPtr_;
138 }
139 
140 
142 {
143  if (!faceAreasPtr_)
144  {
145  calcFaceCentresAndAreas();
146  }
147 
148  return *faceAreasPtr_;
149 }
150 
151 
152 // ************************ vim: set sw=4 sts=4 et: ************************ //