cell.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 \*---------------------------------------------------------------------------*/
25 
26 #include "cell.H"
27 #include "pyramidPointFaceRef.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 const char* const Foam::cell::typeName = "cell";
32 
33 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
34 
36 {
37  // return the unordered list of vertex labels supporting the cell
38 
39  // count the maximum size of all vertices
40  label maxVert = 0;
41 
42  const labelList& faces = *this;
43 
44  forAll(faces, facei)
45  {
46  maxVert += f[faces[facei]].size();
47  }
48 
49  // set the fill-in list
50  labelList p(maxVert);
51 
52  // in the first face there is no duplicates
53  const labelList& first = f[faces[0]];
54 
55  forAll(first, pointi)
56  {
57  p[pointi] = first[pointi];
58  }
59 
60  // re-use maxVert to count the real vertices
61  maxVert = first.size();
62 
63  // go through the rest of the faces. For each vertex, check if the point is
64  // already inserted (up to maxVert, which now carries the number of real
65  // points. If not, add it at the end of the list.
66  for (label facei = 1; facei < faces.size(); facei++)
67  {
68  const labelList& curFace = f[faces[facei]];
69 
70  forAll(curFace, pointi)
71  {
72  const label curPoint = curFace[pointi];
73 
74  bool found = false;
75 
76  for (label checkI = 0; checkI < maxVert; checkI++)
77  {
78  if (curPoint == p[checkI])
79  {
80  found = true;
81  break;
82  }
83  }
84 
85  if (!found)
86  {
87  p[maxVert] = curPoint;
88 
89  maxVert++;
90  }
91  }
92  }
93 
94  // reset the size of the list
95  p.setSize(maxVert);
96 
97  return p;
98 }
99 
100 
102 (
103  const faceUList& f,
104  const pointField& meshPoints
105 ) const
106 {
108 
109  pointField p(pointLabels.size());
110 
111  forAll(p, i)
112  {
113  p[i] = meshPoints[pointLabels[i]];
114  }
115 
116  return p;
117 }
118 
119 
121 {
122  // return the lisf of cell edges
123 
124  const labelList& curFaces = *this;
125 
126  // create a list of edges
127  label maxNoEdges = 0;
128 
129  forAll(curFaces, facei)
130  {
131  maxNoEdges += f[curFaces[facei]].nEdges();
132  }
133 
134  edgeList allEdges(maxNoEdges);
135  label nEdges = 0;
136 
137  forAll(curFaces, facei)
138  {
139  const edgeList curFaceEdges = f[curFaces[facei]].edges();
140 
141  forAll(curFaceEdges, faceEdgeI)
142  {
143  const edge& curEdge = curFaceEdges[faceEdgeI];
144 
145  bool edgeFound = false;
146 
147  for (label addedEdgeI = 0; addedEdgeI < nEdges; addedEdgeI++)
148  {
149  if (allEdges[addedEdgeI] == curEdge)
150  {
151  edgeFound = true;
152 
153  break;
154  }
155  }
156 
157  if (!edgeFound)
158  {
159  // Add the new edge onto the list
160  allEdges[nEdges] = curEdge;
161  nEdges++;
162  }
163  }
164  }
165 
166  allEdges.setSize(nEdges);
167 
168  return allEdges;
169 }
170 
171 
173 (
174  const pointField& p,
175  const faceUList& f
176 ) const
177 {
178  // When one wants to access the cell centre and magnitude, the
179  // functionality on the mesh level should be used in preference to the
180  // functions provided here. They do not rely to the functionality
181  // implemented here, provide additional checking and are more efficient.
182  // The cell::centre and cell::mag functions may be removed in the future.
183 
184  // WARNING!
185  // In the old version of the code, it was possible to establish whether any
186  // of the pyramids had a negative volume, caused either by the poor
187  // estimate of the cell centre or by the fact that the cell was defined
188  // inside out. In the new description of the cell, this can only be
189  // established with the reference to the owner-neighbour face-cell
190  // relationship, which is not available on this level. Thus, all the
191  // pyramids are believed to be positive with no checking.
192 
193  // first calculate the approximate cell centre as the average of all
194  // face centres
195  vector cEst = Zero;
196  scalar sumArea = 0;
197 
198  const labelList& faces = *this;
199 
200  forAll(faces, facei)
201  {
202  scalar a = f[faces[facei]].mag(p);
203  cEst += f[faces[facei]].centre(p)*a;
204  sumArea += a;
205  }
206 
207  cEst /= sumArea + vSmall;
208 
209  // Calculate the centre by breaking the cell into pyramids and
210  // volume-weighted averaging their centres
211  vector sumVc = Zero;
212 
213  scalar sumV = 0;
214 
215  forAll(faces, facei)
216  {
217  // calculate pyramid volume. If it is greater than zero, OK.
218  // If not, the pyramid is inside-out. Create a face with the opposite
219  // order and recalculate pyramid centre!
220  scalar pyrVol = pyramidPointFaceRef(f[faces[facei]], cEst).mag(p);
221  vector pyrCentre = pyramidPointFaceRef(f[faces[facei]], cEst).centre(p);
222 
223  // if pyramid inside-out because face points inwards invert
224  // N.B. pyramid remains unchanged
225  if (pyrVol < 0)
226  {
227  pyrVol = -pyrVol;
228  }
229 
230  sumVc += pyrVol*pyrCentre;
231  sumV += pyrVol;
232  }
233 
234  return sumVc/(sumV + vSmall);
235 }
236 
237 
238 Foam::scalar Foam::cell::mag
239 (
240  const pointField& p,
241  const faceUList& f
242 ) const
243 {
244  // When one wants to access the cell centre and magnitude, the
245  // functionality on the mesh level should be used in preference to the
246  // functions provided here. They do not rely to the functionality
247  // implemented here, provide additional checking and are more efficient.
248  // The cell::centre and cell::mag functions may be removed in the future.
249 
250  // WARNING! See cell::centre
251 
252  // first calculate the approximate cell centre as the average of all
253  // face centres
254  vector cEst = Zero;
255  scalar nCellFaces = 0;
256 
257  const labelList& faces = *this;
258 
259  forAll(faces, facei)
260  {
261  cEst += f[faces[facei]].centre(p);
262  nCellFaces += 1;
263  }
264 
265  cEst /= nCellFaces;
266 
267  // Calculate the magnitude by summing the mags of the pyramids
268  scalar v = 0;
269 
270  forAll(faces, facei)
271  {
272  v += ::Foam::mag(pyramidPointFaceRef(f[faces[facei]], cEst).mag(p));
273  }
274 
275  return v;
276 }
277 
278 
279 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
280 
281 bool Foam::operator==(const cell& a, const cell& b)
282 {
283  // Trivial reject: faces are different size
284  if (a.size() != b.size())
285  {
286  return false;
287  }
288 
289  List<bool> fnd(a.size(), false);
290 
291  forAll(b, bI)
292  {
293  label curLabel = b[bI];
294 
295  bool found = false;
296 
297  forAll(a, aI)
298  {
299  if (a[aI] == curLabel)
300  {
301  found = true;
302  fnd[aI] = true;
303  break;
304  }
305  }
306 
307  if (!found)
308  {
309  return false;
310  }
311  }
312 
313  // check if all faces on a were marked
314  bool result = true;
315 
316  forAll(fnd, aI)
317  {
318  result = (result && fnd[aI]);
319  }
320 
321  return result;
322 }
323 
324 
325 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
labelList pointLabels(nPoints, -1)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
T & first()
Return the first element of the list.
Definition: UListI.H:114
edgeList edges(const faceUList &) const
Return cell edges.
Definition: cell.C:120
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
static const zero Zero
Definition: zero.H:97
point centre(const pointField &, const faceUList &) const
Returns cell centre.
Definition: cell.C:173
pyramid< point, const point &, const face & > pyramidPointFaceRef
scalar mag(const pointField &, const faceUList &) const
Returns cell volume.
Definition: cell.C:239
void setSize(const label)
Reset size of List.
Definition: List.C:281
labelList labels(const faceUList &) const
Return labels of cell vertices.
Definition: cell.C:35
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
pointField points(const faceUList &, const pointField &) const
Return the cell vertices.
Definition: cell.C:102
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField & p
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
bool found
static const char *const typeName
Definition: cell.H:65