face.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-2021 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 "face.H"
27 #include "triFace.H"
28 #include "triPointRef.H"
29 #include "mathematicalConstants.H"
30 #include "Swap.H"
31 #include "ConstCirculator.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 const char* const Foam::face::typeName = "face";
36 
37 
38 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
39 
41 :
42  labelList(f)
43 {}
44 
45 
46 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
47 
48 int Foam::face::compare(const face& a, const face& b)
49 {
50  // Basic rule: we assume that the sequence of labels in each list
51  // will be circular in the same order (but not necessarily in the
52  // same direction or from the same starting point).
53 
54  // Trivial reject: faces are different size
55  label sizeA = a.size();
56  label sizeB = b.size();
57 
58  if (sizeA != sizeB || sizeA == 0)
59  {
60  return 0;
61  }
62  else if (sizeA == 1)
63  {
64  if (a[0] == b[0])
65  {
66  return 1;
67  }
68  else
69  {
70  return 0;
71  }
72  }
73 
74  ConstCirculator<face> aCirc(a);
75  ConstCirculator<face> bCirc(b);
76 
77  // Rotate face b until its element matches the starting element of face a.
78  do
79  {
80  if (aCirc() == bCirc())
81  {
82  // Set bCirc fulcrum to its iterator and increment the iterators
83  bCirc.setFulcrumToIterator();
84  ++aCirc;
85  ++bCirc;
86 
87  break;
88  }
90 
91  // If the circulator has stopped then faces a and b do not share a matching
92  // point. Doesn't work on matching, single element face.
93  if (!bCirc.circulate())
94  {
95  return 0;
96  }
97 
98  // Look forwards around the faces for a match
99  do
100  {
101  if (aCirc() != bCirc())
102  {
103  break;
104  }
105  }
106  while
107  (
110  );
111 
112  // If the circulator has stopped then faces a and b matched.
113  if (!aCirc.circulate())
114  {
115  return 1;
116  }
117  else
118  {
119  // Reset the circulators back to their fulcrum
120  aCirc.setIteratorToFulcrum();
121  bCirc.setIteratorToFulcrum();
122  ++aCirc;
123  --bCirc;
124  }
125 
126  // Look backwards around the faces for a match
127  do
128  {
129  if (aCirc() != bCirc())
130  {
131  break;
132  }
133  }
134  while
135  (
138  );
139 
140  // If the circulator has stopped then faces a and b matched.
141  if (!aCirc.circulate())
142  {
143  return -1;
144  }
145 
146  return 0;
147 }
148 
149 
150 bool Foam::face::sameVertices(const face& a, const face& b)
151 {
152  label sizeA = a.size();
153  label sizeB = b.size();
154 
155  // Trivial reject: faces are different size
156  if (sizeA != sizeB)
157  {
158  return false;
159  }
160  // Check faces with a single vertex
161  else if (sizeA == 1)
162  {
163  if (a[0] == b[0])
164  {
165  return true;
166  }
167  else
168  {
169  return false;
170  }
171  }
172 
173  forAll(a, i)
174  {
175  // Count occurrences of a[i] in a
176  label aOcc = 0;
177  forAll(a, j)
178  {
179  if (a[i] == a[j]) aOcc++;
180  }
181 
182  // Count occurrences of a[i] in b
183  label bOcc = 0;
184  forAll(b, j)
185  {
186  if (a[i] == b[j]) bOcc++;
187  }
188 
189  // Check if occurrences of a[i] in a and b are the same
190  if (aOcc != bOcc) return false;
191  }
192 
193  return true;
194 }
195 
196 
197 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
198 
200 {
201  if (size() > 1)
202  {
203  label ci = 0;
204  for (label i=1; i<size(); i++)
205  {
206  if (operator[](i) != operator[](ci))
207  {
208  operator[](++ci) = operator[](i);
209  }
210  }
211 
212  if (operator[](ci) != operator[](0))
213  {
214  ci++;
215  }
216 
217  setSize(ci);
218  }
219 
220  return size();
221 }
222 
223 
225 {
226  const label n = size();
227 
228  if (n > 2)
229  {
230  for (label i=1; i < (n+1)/2; ++i)
231  {
232  Swap(operator[](i), operator[](n-i));
233  }
234  }
235 }
236 
237 
239 {
240  return centre(UIndirectList<point>(ps, *this));
241 }
242 
243 
245 {
246  return area(UIndirectList<point>(ps, *this));
247 }
248 
249 
251 {
252  return normalised(area(points));
253 }
254 
255 
257 {
258  // Reverse the label list and return
259  // The starting points of the original and reverse face are identical.
260 
261  const labelList& f = *this;
262  labelList newList(size());
263 
264  newList[0] = f[0];
265 
266  for (label pointi = 1; pointi < newList.size(); pointi++)
267  {
268  newList[pointi] = f[size() - pointi];
269  }
270 
271  return face(move(newList));
272 }
273 
274 
276 {
277  const labelList& f = *this;
278 
279  forAll(f, localIdx)
280  {
281  if (f[localIdx] == globalIndex)
282  {
283  return localIdx;
284  }
285  }
286 
287  return -1;
288 }
289 
290 
291 Foam::scalar Foam::face::sweptVol
292 (
293  const pointField& oldPoints,
294  const pointField& newPoints
295 ) const
296 {
297  // This Optimisation causes a small discrepancy between the swept-volume of
298  // opposite faces of complex cells with triangular faces opposing polygons.
299  // It could be used without problem for tetrahedral cells
300  // if (size() == 3)
301  // {
302  // return
303  // (
304  // triPointRef
305  // (
306  // oldPoints[operator[](0)],
307  // oldPoints[operator[](1)],
308  // oldPoints[operator[](2)]
309  // ).sweptVol
310  // (
311  // triPointRef
312  // (
313  // newPoints[operator[](0)],
314  // newPoints[operator[](1)],
315  // newPoints[operator[](2)]
316  // )
317  // )
318  // );
319  // }
320 
321  scalar sv = 0;
322 
323  // Calculate the swept volume by breaking the face into triangles and
324  // summing their swept volumes.
325  // Changed to deal with small concavity by using a central decomposition
326 
327  point centreOldPoint = centre(oldPoints);
328  point centreNewPoint = centre(newPoints);
329 
330  label nPoints = size();
331 
332  for (label pi=0; pi<nPoints-1; ++pi)
333  {
334  // Note: for best accuracy, centre point always comes last
335  sv += triPointRef
336  (
337  centreOldPoint,
338  oldPoints[operator[](pi)],
339  oldPoints[operator[](pi + 1)]
340  ).sweptVol
341  (
343  (
344  centreNewPoint,
345  newPoints[operator[](pi)],
346  newPoints[operator[](pi + 1)]
347  )
348  );
349  }
350 
351  sv += triPointRef
352  (
353  centreOldPoint,
354  oldPoints[operator[](nPoints-1)],
355  oldPoints[operator[](0)]
356  ).sweptVol
357  (
359  (
360  centreNewPoint,
361  newPoints[operator[](nPoints-1)],
362  newPoints[operator[](0)]
363  )
364  );
365 
366  return sv;
367 }
368 
369 
371 (
372  const pointField& p,
373  const point& refPt,
374  scalar density
375 ) const
376 {
377  // If the face is a triangle, do a direct calculation
378  if (size() == 3)
379  {
380  return triPointRef
381  (
382  p[operator[](0)],
383  p[operator[](1)],
384  p[operator[](2)]
385  ).inertia(refPt, density);
386  }
387 
388  const point ctr = centre(p);
389 
390  tensor J = Zero;
391 
392  forAll(*this, i)
393  {
394  J += triPointRef
395  (
396  p[operator[](i)],
397  p[operator[](fcIndex(i))],
398  ctr
399  ).inertia(refPt, density);
400  }
401 
402  return J;
403 }
404 
405 
407 {
408  const labelList& points = *this;
409 
410  edgeList e(points.size());
411 
412  for (label pointi = 0; pointi < points.size() - 1; ++pointi)
413  {
414  e[pointi] = edge(points[pointi], points[pointi + 1]);
415  }
416 
417  // Add last edge
418  e.last() = edge(points.last(), points[0]);
419 
420  return e;
421 }
422 
423 
425 {
426  forAll(*this, i)
427  {
428  if (operator[](i) == e.start())
429  {
430  if (operator[](rcIndex(i)) == e.end())
431  {
432  // Reverse direction
433  return -1;
434  }
435  else if (operator[](fcIndex(i)) == e.end())
436  {
437  // Forward direction
438  return 1;
439  }
440 
441  // No match
442  return 0;
443  }
444  else if (operator[](i) == e.end())
445  {
446  if (operator[](rcIndex(i)) == e.start())
447  {
448  // Forward direction
449  return 1;
450  }
451  else if (operator[](fcIndex(i)) == e.start())
452  {
453  // Reverse direction
454  return -1;
455  }
456 
457  // No match
458  return 0;
459  }
460  }
461 
462  // Not found
463  return 0;
464 }
465 
466 
468 {
469  const edgeList& eds = f.edges();
470 
471  label longestEdgeI = -1;
472  scalar longestEdgeLength = -small;
473 
474  forAll(eds, edI)
475  {
476  scalar edgeLength = eds[edI].mag(pts);
477 
478  if (edgeLength > longestEdgeLength)
479  {
480  longestEdgeI = edI;
481  longestEdgeLength = edgeLength;
482  }
483  }
484 
485  return longestEdgeI;
486 }
487 
488 
489 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
edgeList edges() const
Return edges in face point ordering,.
Definition: face.C:406
T & operator[](const label)
Return element of UList.
Definition: UListI.H:167
face()
Construct null.
Definition: faceI.H:28
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
tensor inertia(const pointField &, const point &refPt=vector::zero, scalar density=1.0) const
Return the inertia tensor, with optional reference.
Definition: face.C:371
static vector area(const PointField &ps)
Return vector area given face points.
label which(const label globalIndex) const
Navigation through face vertices.
Definition: face.C:275
void setIteratorToFulcrum()
Set the iterator to the current position of the fulcrum.
vector normal(const pointField &) const
Return unit normal.
Definition: face.C:250
label collapse()
Collapse face by removing duplicate point labels.
Definition: face.C:199
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
face reverseFace() const
Return face with reverse direction.
Definition: face.C:256
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
void Swap(T &a, T &b)
Definition: Swap.H:43
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
int edgeDirection(const edge &) const
Return the edge direction on the face.
Definition: face.C:424
label nPoints
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
static int compare(const face &, const face &)
Compare faces.
Definition: face.C:48
static const zero Zero
Definition: zero.H:97
pointField points(const pointField &) const
Return the points corresponding to this face.
Definition: faceI.H:64
static const char *const typeName
Definition: face.H:91
static vector centre(const PointField &ps)
Return centre point given face points.
labelList f(nPoints)
scalar sweptVol(const pointField &oldPoints, const pointField &newPoints) const
Return the volume swept out by the face when its points move.
Definition: face.C:292
label longestEdge(const face &f, const pointField &pts)
Find the longest edge on a face. Face point labels index into pts.
Definition: face.C:467
static bool sameVertices(const face &, const face &)
Return true if the faces have the same vertices.
Definition: face.C:150
void setSize(const label)
Reset size of List.
Definition: List.C:281
bool circulate(const CirculatorBase::direction dir=CirculatorBase::direction::none)
Circulate around the list in the given direction.
label end() const
Return end vertex label.
Definition: edgeI.H:92
A List with indirect addressing.
Definition: fvMatrix.H:106
Swap its arguments.
void flip()
Flip the face in-place.
Definition: face.C:224
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
label n
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Walks over a container as if it were circular. The container must have the following members defined:...
T & last()
Return the last element of the list.
Definition: UListI.H:128
void setFulcrumToIterator()
Set the fulcrum to the current position of the iterator.
label size() const
Return the number of elements in the UList.
Definition: ListI.H:171
label start() const
Return start vertex label.
Definition: edgeI.H:81