FreeFOAM The Cross-Platform CFD Toolkit
face.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 \*---------------------------------------------------------------------------*/
25 
26 #include "face.H"
27 #include <OpenFOAM/triFace.H>
28 #include <OpenFOAM/triPointRef.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 const char* const Foam::face::typeName = "face";
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
38 Foam::face::calcEdges(const pointField& points) const
39 {
40  tmp<vectorField> tedges(new vectorField(size()));
41  vectorField& edges = tedges();
42 
43  forAll(*this, i)
44  {
45  label ni = fcIndex(i);
46 
47  point thisPt = points[operator[](i)];
48  point nextPt = points[operator[](ni)];
49 
50  vector vec(nextPt - thisPt);
51  vec /= Foam::mag(vec) + VSMALL;
52 
53  edges[i] = vec;
54  }
55 
56  return tedges;
57 }
58 
59 
60 Foam::scalar Foam::face::edgeCos
61 (
62  const vectorField& edges,
63  const label index
64 ) const
65 {
66  label leftEdgeI = left(index);
67  label rightEdgeI = right(index);
68 
69  // note negate on left edge to get correct left-pointing edge.
70  return -(edges[leftEdgeI] & edges[rightEdgeI]);
71 }
72 
73 
74 Foam::label Foam::face::mostConcaveAngle
75 (
76  const pointField& points,
77  const vectorField& edges,
78  scalar& maxAngle
79 ) const
80 {
81  vector n(normal(points));
82 
83  label index = 0;
84  maxAngle = -GREAT;
85 
86  forAll(edges, i)
87  {
88  label leftEdgeI = left(i);
89  label rightEdgeI = right(i);
90 
91  vector edgeNormal = edges[rightEdgeI] ^ edges[leftEdgeI];
92 
93  scalar edgeCos = edges[leftEdgeI] & edges[rightEdgeI];
94  scalar edgeAngle = acos(max(-1.0, min(1.0, edgeCos)));
95 
96  scalar angle;
97 
98  if ((edgeNormal & n) > 0)
99  {
100  // Concave angle.
101  angle = mathematicalConstant::pi + edgeAngle;
102  }
103  else
104  {
105  // Convex angle. Note '-' to take into account that rightEdge
106  // and leftEdge are head-to-tail connected.
107  angle = mathematicalConstant::pi - edgeAngle;
108  }
109 
110  if (angle > maxAngle)
111  {
112  maxAngle = angle;
113  index = i;
114  }
115  }
116 
117  return index;
118 }
119 
120 
121 Foam::label Foam::face::split
122 (
123  const face::splitMode mode,
124  const pointField& points,
125  label& triI,
126  label& quadI,
127  faceList& triFaces,
128  faceList& quadFaces
129 ) const
130 {
131  label oldIndices = (triI + quadI);
132 
133  if (size() <= 2)
134  {
136  (
137  "face::split"
138  "(const face::splitMode, const pointField&, label&, label&"
139  ", faceList&, faceList&)"
140  )
141  << "Serious problem: asked to split a face with < 3 vertices"
142  << abort(FatalError);
143  }
144  if (size() == 3)
145  {
146  // Triangle. Just copy.
147  if (mode == COUNTTRIANGLE || mode == COUNTQUAD)
148  {
149  triI++;
150  }
151  else
152  {
153  triFaces[triI++] = *this;
154  }
155  }
156  else if (size() == 4)
157  {
158  if (mode == COUNTTRIANGLE)
159  {
160  triI += 2;
161  }
162  if (mode == COUNTQUAD)
163  {
164  quadI++;
165  }
166  else if (mode == SPLITTRIANGLE)
167  {
168  // Start at point with largest internal angle.
169  const vectorField edges(calcEdges(points));
170 
171  scalar minAngle;
172  label startIndex = mostConcaveAngle(points, edges, minAngle);
173 
174  label nextIndex = fcIndex(startIndex);
175  label splitIndex = fcIndex(nextIndex);
176 
177  // Create triangles
178  face triFace(3);
179  triFace[0] = operator[](startIndex);
180  triFace[1] = operator[](nextIndex);
181  triFace[2] = operator[](splitIndex);
182 
183  triFaces[triI++] = triFace;
184 
185  triFace[0] = operator[](splitIndex);
186  triFace[1] = operator[](fcIndex(splitIndex));
187  triFace[2] = operator[](startIndex);
188 
189  triFaces[triI++] = triFace;
190  }
191  else
192  {
193  quadFaces[quadI++] = *this;
194  }
195  }
196  else
197  {
198  // General case. Like quad: search for largest internal angle.
199 
200  const vectorField edges(calcEdges(points));
201 
202  scalar minAngle = 1;
203  label startIndex = mostConcaveAngle(points, edges, minAngle);
204 
205  scalar bisectAngle = minAngle/2;
206  vector rightEdge = edges[right(startIndex)];
207 
208  //
209  // Look for opposite point which as close as possible bisects angle
210  //
211 
212  // split candidate starts two points away.
213  label index = fcIndex(fcIndex(startIndex));
214 
215  label minIndex = index;
216  scalar minDiff = Foam::mathematicalConstant::pi;
217 
218  for(label i = 0; i < size() - 3; i++)
219  {
220  vector splitEdge
221  (
222  points[operator[](index)]
223  - points[operator[](startIndex)]
224  );
225  splitEdge /= Foam::mag(splitEdge) + VSMALL;
226 
227  const scalar splitCos = splitEdge & rightEdge;
228  const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
229  const scalar angleDiff = fabs(splitAngle - bisectAngle);
230 
231  if (angleDiff < minDiff)
232  {
233  minDiff = angleDiff;
234  minIndex = index;
235  }
236 
237  // go to next candidate
238  index = fcIndex(index);
239  }
240 
241 
242  // Split into two subshapes.
243  // face1: startIndex to minIndex
244  // face2: minIndex to startIndex
245 
246  // Get sizes of the two subshapes
247  label diff = 0;
248  if (minIndex > startIndex)
249  {
250  diff = minIndex - startIndex;
251  }
252  else
253  {
254  // folded around
255  diff = minIndex + size() - startIndex;
256  }
257 
258  label nPoints1 = diff + 1;
259  label nPoints2 = size() - diff + 1;
260 
261  // Collect face1 points
262  face face1(nPoints1);
263 
264  index = startIndex;
265  for (label i = 0; i < nPoints1; i++)
266  {
267  face1[i] = operator[](index);
268  index = fcIndex(index);
269  }
270 
271  // Collect face2 points
272  face face2(nPoints2);
273 
274  index = minIndex;
275  for (label i = 0; i < nPoints2; i++)
276  {
277  face2[i] = operator[](index);
278  index = fcIndex(index);
279  }
280 
281  // Split faces
282  face1.split(mode, points, triI, quadI, triFaces, quadFaces);
283  face2.split(mode, points, triI, quadI, triFaces, quadFaces);
284  }
285 
286  return (triI + quadI - oldIndices);
287 }
288 
289 
290 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
291 
293 :
294  labelList(f)
295 {}
296 
297 
298 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
299 
300 
301 // return
302 // 0: no match
303 // +1: identical
304 // -1: same face, but different orientation
305 int Foam::face::compare(const face& a, const face& b)
306 {
307  // Basic rule: we assume that the sequence of labels in each list
308  // will be circular in the same order (but not necessarily in the
309  // same direction or from the same starting point).
310 
311  // Trivial reject: faces are different size
312  label sizeA = a.size();
313  label sizeB = b.size();
314 
315  if (sizeA != sizeB)
316  {
317  return 0;
318  }
319 
320 
321  // Full list comparison
322  const label firstA = a[0];
323  label Bptr = -1;
324 
325  forAll (b, i)
326  {
327  if (b[i] == firstA)
328  {
329  Bptr = i; // 'found match' at element 'i'
330  break;
331  }
332  }
333 
334  // If no match was found, return 0
335  if (Bptr < 0)
336  {
337  return 0;
338  }
339 
340  // Now we must look for the direction, if any
341  label secondA = a[1];
342 
343  if (sizeA > 1 && (secondA == firstA || firstA == a[sizeA - 1]))
344  {
345  face ca = a;
346  ca.collapse();
347 
348  face cb = b;
349  cb.collapse();
350 
351  return face::compare(ca, cb);
352  }
353 
354  int dir = 0;
355 
356  // Check whether at top of list
357  Bptr++;
358  if (Bptr == b.size())
359  {
360  Bptr = 0;
361  }
362 
363  // Test whether upward label matches second A label
364  if (b[Bptr] == secondA)
365  {
366  // Yes - direction is 'up'
367  dir = 1;
368  }
369  else
370  {
371  // No - so look downwards, checking whether at bottom of list
372  Bptr -= 2;
373 
374  if (Bptr < 0)
375  {
376  // wraparound
377  Bptr += b.size();
378  }
379 
380  // Test whether downward label matches second A label
381  if (b[Bptr] == secondA)
382  {
383  // Yes - direction is 'down'
384  dir = -1;
385  }
386  }
387 
388  // Check whether a match was made at all, and exit 0 if not
389  if (dir == 0)
390  {
391  return 0;
392  }
393 
394  // Decrement size by 2 to account for first searches
395  sizeA -= 2;
396 
397  // We now have both direction of search and next element
398  // to search, so we can continue search until no more points.
399  label Aptr = 1;
400  if (dir > 0)
401  {
402  while (sizeA--)
403  {
404  Aptr++;
405  if (Aptr >= a.size())
406  {
407  Aptr = 0;
408  }
409 
410  Bptr++;
411  if (Bptr >= b.size())
412  {
413  Bptr = 0;
414  }
415 
416  if (a[Aptr] != b[Bptr])
417  {
418  return 0;
419  }
420  }
421  }
422  else
423  {
424  while (sizeA--)
425  {
426  Aptr++;
427  if (Aptr >= a.size())
428  {
429  Aptr = 0;
430  }
431 
432  Bptr--;
433  if (Bptr < 0)
434  {
435  Bptr = b.size() - 1;
436  }
437 
438  if (a[Aptr] != b[Bptr])
439  {
440  return 0;
441  }
442  }
443  }
444 
445  // They must be equal - return direction
446  return dir;
447 }
448 
449 
450 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
451 
452 
453 Foam::label Foam::face::collapse()
454 {
455  if (size() > 1)
456  {
457  label ci = 0;
458  for (label i=1; i<size(); i++)
459  {
460  if (operator[](i) != operator[](ci))
461  {
462  operator[](++ci) = operator[](i);
463  }
464  }
465 
466  if (operator[](ci) != operator[](0))
467  {
468  ci++;
469  }
470 
471  setSize(ci);
472  }
473 
474  return size();
475 }
476 
477 
478 Foam::point Foam::face::centre(const pointField& meshPoints) const
479 {
480  // Calculate the centre by breaking the face into triangles and
481  // area-weighted averaging their centres
482 
483  // If the face is a triangle, do a direct calculation
484  if (size() == 3)
485  {
486  return
487  (1.0/3.0)
488  *(
489  meshPoints[operator[](0)]
490  + meshPoints[operator[](1)]
491  + meshPoints[operator[](2)]
492  );
493  }
494 
495  label nPoints = size();
496 
497  point centrePoint = point::zero;
498  for (register label pI=0; pI<nPoints; pI++)
499  {
500  centrePoint += meshPoints[operator[](pI)];
501  }
502  centrePoint /= nPoints;
503 
504  scalar sumA = 0;
505  vector sumAc = vector::zero;
506 
507  for (register label pI=0; pI<nPoints; pI++)
508  {
509  const point& nextPoint = meshPoints[operator[]((pI + 1) % nPoints)];
510 
511  // Calculate 3*triangle centre
512  vector ttc
513  (
514  meshPoints[operator[](pI)]
515  + nextPoint
516  + centrePoint
517  );
518 
519  // Calculate 2*triangle area
520  scalar ta = Foam::mag
521  (
522  (meshPoints[operator[](pI)] - centrePoint)
523  ^ (nextPoint - centrePoint)
524  );
525 
526  sumA += ta;
527  sumAc += ta*ttc;
528  }
529 
530  if (sumA > VSMALL)
531  {
532  return sumAc/(3*sumA);
533  }
534  else
535  {
536  return centrePoint;
537  }
538 }
539 
540 
542 {
543  // Calculate the normal by summing the face triangle normals.
544  // Changed to deal with small concavity by using a central decomposition
545  //
546 
547  // If the face is a triangle, do a direct calculation to avoid round-off
548  // error-related problems
549  //
550  if (size() == 3)
551  {
552  return triPointRef
553  (
554  p[operator[](0)],
555  p[operator[](1)],
556  p[operator[](2)]
557  ).normal();
558  }
559 
560  label nPoints = size();
561 
562  register label pI;
563 
564  point centrePoint = vector::zero;
565  for (pI = 0; pI < nPoints; pI++)
566  {
567  centrePoint += p[operator[](pI)];
568  }
569  centrePoint /= nPoints;
570 
571  vector n = vector::zero;
572 
573  point nextPoint = centrePoint;
574 
575  for (pI = 0; pI < nPoints; pI++)
576  {
577  if (pI < nPoints - 1)
578  {
579  nextPoint = p[operator[](pI + 1)];
580  }
581  else
582  {
583  nextPoint = p[operator[](0)];
584  }
585 
586  // Note: for best accuracy, centre point always comes last
587  //
588  n += triPointRef
589  (
590  p[operator[](pI)],
591  nextPoint,
592  centrePoint
593  ).normal();
594  }
595 
596  return n;
597 }
598 
599 
601 {
602  // reverse the label list and return
603  // The starting points of the original and reverse face are identical.
604 
605  const labelList& f = *this;
606  labelList newList(size());
607 
608  newList[0] = f[0];
609 
610  for (label pointI = 1; pointI < newList.size(); pointI++)
611  {
612  newList[pointI] = f[size() - pointI];
613  }
614 
615  return face(xferMove(newList));
616 }
617 
618 
619 Foam::label Foam::face::which(const label globalIndex) const
620 {
621  label pointInFace = -1;
622  const labelList& f = *this;
623 
624  forAll (f, i)
625  {
626  if (f[i] == globalIndex)
627  {
628  pointInFace = i;
629  break;
630  }
631  }
632 
633  return pointInFace;
634 }
635 
636 
637 Foam::scalar Foam::face::sweptVol
638 (
639  const pointField& oldPoints,
640  const pointField& newPoints
641 ) const
642 {
643  scalar sv = 0;
644 
645  // Calculate the swept volume by breaking the face into triangles and
646  // summing their swept volumes.
647  // Changed to deal with small concavity by using a central decomposition
648 
649  point centreOldPoint = centre(oldPoints);
650  point centreNewPoint = centre(newPoints);
651 
652  label nPoints = size();
653 
654  point nextOldPoint = centreOldPoint;
655  point nextNewPoint = centreNewPoint;
656 
657  register label pI;
658 
659  for (pI = 0; pI < nPoints; pI++)
660  {
661  if (pI < nPoints - 1)
662  {
663  nextOldPoint = oldPoints[operator[](pI + 1)];
664  nextNewPoint = newPoints[operator[](pI + 1)];
665  }
666  else
667  {
668  nextOldPoint = oldPoints[operator[](0)];
669  nextNewPoint = newPoints[operator[](0)];
670  }
671 
672  // Note: for best accuracy, centre point always comes last
673  sv += triPointRef
674  (
675  centreOldPoint,
676  oldPoints[operator[](pI)],
677  nextOldPoint
678  ).sweptVol
679  (
681  (
682  centreNewPoint,
683  newPoints[operator[](pI)],
684  nextNewPoint
685  )
686  );
687  }
688 
689  return sv;
690 }
691 
692 
694 {
695  const labelList& points = *this;
696 
697  edgeList e(points.size());
698 
699  label pointI;
700 
701  for (pointI = 0; pointI < points.size() - 1; pointI++)
702  {
703  e[pointI] = edge(points[pointI], points[pointI + 1]);
704  }
705 
706  // add last edge
707  e[points.size() - 1] = edge(points[points.size() - 1], points[0]);
708 
709  return e;
710 }
711 
712 
714 {
715  forAll (*this, i)
716  {
717  if (operator[](i) == e.start())
718  {
719  if (operator[](rcIndex(i)) == e.end())
720  {
721  // reverse direction
722  return -1;
723  }
724  else if (operator[](fcIndex(i)) == e.end())
725  {
726  // forward direction
727  return 1;
728  }
729 
730  // no match
731  return 0;
732  }
733  else if (operator[](i) == e.end())
734  {
735  if (operator[](rcIndex(i)) == e.start())
736  {
737  // forward direction
738  return 1;
739  }
740  else if (operator[](fcIndex(i)) == e.start())
741  {
742  // reverse direction
743  return -1;
744  }
745 
746  // no match
747  return 0;
748  }
749  }
750 
751  // not found
752  return 0;
753 }
754 
755 
756 // Number of triangles directly known from number of vertices
757 Foam::label Foam::face::nTriangles(const pointField&) const
758 {
759  return nTriangles();
760 }
761 
762 
763 Foam::label Foam::face::triangles
764 (
765  const pointField& points,
766  label& triI,
767  faceList& triFaces
768 ) const
769 {
770  label quadI = 0;
771  faceList quadFaces;
772 
773  return split(SPLITTRIANGLE, points, triI, quadI, triFaces, quadFaces);
774 }
775 
776 
777 Foam::label Foam::face::nTrianglesQuads
778 (
779  const pointField& points,
780  label& triI,
781  label& quadI
782 ) const
783 {
784  faceList triFaces;
785  faceList quadFaces;
786 
787  return split(COUNTQUAD, points, triI, quadI, triFaces, quadFaces);
788 }
789 
790 
791 Foam::label Foam::face::trianglesQuads
792 (
793  const pointField& points,
794  label& triI,
795  label& quadI,
796  faceList& triFaces,
797  faceList& quadFaces
798 ) const
799 {
800  return split(SPLITQUAD, points, triI, quadI, triFaces, quadFaces);
801 }
802 
803 
804 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
805 
806 
807 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
808 
809 // ************************ vim: set sw=4 sts=4 et: ************************ //
810