orientedSurface.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 "orientedSurface.H"
27 #include "triSurfaceTools.H"
28 #include "triSurfaceSearch.H"
29 #include "treeBoundBox.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 defineTypeNameAndDebug(orientedSurface, 0);
36 }
37 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 // Return true if edge is used in opposite order in faces
42 bool Foam::orientedSurface::consistentEdge
43 (
44  const edge& e,
45  const triSurface::FaceType& f0,
46  const triSurface::FaceType& f1
47 )
48 {
49  return (f0.edgeDirection(e) > 0) ^ (f1.edgeDirection(e) > 0);
50 }
51 
52 
53 Foam::labelList Foam::orientedSurface::faceToEdge
54 (
55  const triSurface& s,
56  const labelList& changedFaces
57 )
58 {
59  labelList changedEdges(3*changedFaces.size());
60  label changedI = 0;
61 
62  forAll(changedFaces, i)
63  {
64  const labelList& fEdges = s.faceEdges()[changedFaces[i]];
65 
66  forAll(fEdges, j)
67  {
68  changedEdges[changedI++] = fEdges[j];
69  }
70  }
71  changedEdges.setSize(changedI);
72 
73  return changedEdges;
74 }
75 
76 
77 Foam::labelList Foam::orientedSurface::edgeToFace
78 (
79  const triSurface& s,
80  const labelList& changedEdges,
81  labelList& flip
82 )
83 {
84  labelList changedFaces(2*changedEdges.size());
85  label changedI = 0;
86 
87  forAll(changedEdges, i)
88  {
89  label edgeI = changedEdges[i];
90 
91  const labelList& eFaces = s.edgeFaces()[edgeI];
92 
93  if (eFaces.size() < 2)
94  {
95  // Do nothing, faces was already visited.
96  }
97  else if (eFaces.size() == 2)
98  {
99  label face0 = eFaces[0];
100  label face1 = eFaces[1];
101 
102  const triSurface::FaceType& f0 = s.localFaces()[face0];
103  const triSurface::FaceType& f1 = s.localFaces()[face1];
104 
105  if (flip[face0] == UNVISITED)
106  {
107  if (flip[face1] == UNVISITED)
108  {
110  << abort(FatalError);
111  }
112  else
113  {
114  // Face1 has a flip state, face0 hasn't
115  if (consistentEdge(s.edges()[edgeI], f0, f1))
116  {
117  // Take over flip status
118  flip[face0] = (flip[face1] == FLIP ? FLIP : NOFLIP);
119  }
120  else
121  {
122  // Invert
123  flip[face0] = (flip[face1] == FLIP ? NOFLIP : FLIP);
124  }
125  changedFaces[changedI++] = face0;
126  }
127  }
128  else
129  {
130  if (flip[face1] == UNVISITED)
131  {
132  // Face0 has a flip state, face1 hasn't
133  if (consistentEdge(s.edges()[edgeI], f0, f1))
134  {
135  flip[face1] = (flip[face0] == FLIP ? FLIP : NOFLIP);
136  }
137  else
138  {
139  flip[face1] = (flip[face0] == FLIP ? NOFLIP : FLIP);
140  }
141  changedFaces[changedI++] = face1;
142  }
143  }
144  }
145  else
146  {
147  // Multiply connected. Do what?
148  }
149  }
150  changedFaces.setSize(changedI);
151 
152  return changedFaces;
153 }
154 
155 
156 void Foam::orientedSurface::walkSurface
157 (
158  const triSurface& s,
159  const label startFacei,
160  labelList& flipState
161 )
162 {
163  // List of faces that were changed in the last iteration.
164  labelList changedFaces(1, startFacei);
165  // List of edges that were changed in the last iteration.
166  labelList changedEdges;
167 
168  while (true)
169  {
170  changedEdges = faceToEdge(s, changedFaces);
171 
172  if (changedEdges.empty())
173  {
174  break;
175  }
176 
177  changedFaces = edgeToFace(s, changedEdges, flipState);
178 
179  if (changedFaces.empty())
180  {
181  break;
182  }
183  }
184 }
185 
186 
187 void Foam::orientedSurface::propagateOrientation
188 (
189  const triSurface& s,
190  const point& samplePoint,
191  const bool orientOutside,
192  const label nearestFacei,
193  const point& nearestPt,
194  labelList& flipState
195 )
196 {
197  //
198  // Determine orientation to normal on nearest face
199  //
201  (
202  s,
203  samplePoint,
204  nearestFacei
205  );
206 
207  if (side == triSurfaceTools::UNKNOWN)
208  {
209  // Non-closed surface. Do what? For now behave as if no flipping
210  // necessary
211  flipState[nearestFacei] = NOFLIP;
212  }
213  else if ((side == triSurfaceTools::OUTSIDE) == orientOutside)
214  {
215  // outside & orientOutside or inside & !orientOutside
216  // Normals on surface pointing correctly. No need to flip normals
217  flipState[nearestFacei] = NOFLIP;
218  }
219  else
220  {
221  // Need to flip normals.
222  flipState[nearestFacei] = FLIP;
223  }
224 
225  if (debug)
226  {
227  vector n = triSurfaceTools::surfaceNormal(s, nearestFacei, nearestPt);
228 
229  Pout<< "orientedSurface::propagateOrientation : starting face"
230  << " orientation:" << nl
231  << " for samplePoint:" << samplePoint << nl
232  << " starting from point:" << nearestPt << nl
233  << " on face:" << nearestFacei << nl
234  << " with normal:" << n << nl
235  << " decided side:" << label(side)
236  << endl;
237  }
238 
239  // Walk the surface from nearestFacei, changing the flipstate.
240  walkSurface(s, nearestFacei, flipState);
241 }
242 
243 
244 // Find side for zoneI only by counting the number of intersections. Determines
245 // if face is oriented consistent with outwards pointing normals.
246 void Foam::orientedSurface::findZoneSide
247 (
248  const triSurfaceSearch& surfSearches,
249  const labelList& faceZone,
250  const label zoneI,
251  const point& outsidePoint,
252  label& zoneFacei,
253  bool& isOutside
254 )
255 {
256  const triSurface& s = surfSearches.surface();
257 
258  zoneFacei = -1;
259  isOutside = false;
260 
261  pointField start(1, outsidePoint);
262  List<List<pointIndexHit>> hits(1, List<pointIndexHit>());
263 
264  forAll(faceZone, facei)
265  {
266  if (faceZone[facei] == zoneI)
267  {
268  const point& fc = s.faceCentres()[facei];
269  const vector& n = s.faceNormals()[facei];
270 
271  const vector d = fc - outsidePoint;
272  const scalar magD = mag(d);
273 
274  // Check if normal different enough to decide upon
275  if (magD > small && (mag(n & d/magD) > 1e-6))
276  {
277  pointField end(1, fc + d);
278 
279  // Info<< "Zone " << zoneI << " : Shooting ray"
280  // << " from " << outsidePoint
281  // << " to " << end
282  // << " to pierce triangle " << facei
283  // << " with centre " << fc << endl;
284 
285 
286  surfSearches.findLineAll(start, end, hits);
287 
288  label zoneIndex = -1;
289  forAll(hits[0], i)
290  {
291  if (hits[0][i].index() == facei)
292  {
293  zoneIndex = i;
294  break;
295  }
296  }
297 
298  if (zoneIndex != -1)
299  {
300  zoneFacei = facei;
301 
302  if ((zoneIndex%2) == 0)
303  {
304  // Odd number of intersections. Check if normal points
305  // in direction of ray
306  isOutside = ((n & d) < 0);
307  }
308  else
309  {
310  isOutside = ((n & d) > 0);
311  }
312  break;
313  }
314  }
315  }
316  }
317 }
318 
319 
320 bool Foam::orientedSurface::flipSurface
321 (
322  triSurface& s,
323  const labelList& flipState
324 )
325 {
326  bool hasFlipped = false;
327 
328  // Flip tris in s
329  forAll(flipState, facei)
330  {
331  if (flipState[facei] == UNVISITED)
332  {
334  << "unvisited face " << facei
335  << abort(FatalError);
336  }
337  else if (flipState[facei] == FLIP)
338  {
339  labelledTri& tri = s[facei];
340 
341  label tmp = tri[0];
342 
343  tri[0] = tri[1];
344  tri[1] = tmp;
345 
346  hasFlipped = true;
347  }
348  }
349  // Recalculate normals
350  if (hasFlipped)
351  {
352  s.clearOut();
353  }
354  return hasFlipped;
355 }
356 
357 
358 bool Foam::orientedSurface::orientConsistent(triSurface& s)
359 {
360  bool anyFlipped = false;
361 
362  // Do initial flipping to make triangles consistent. Otherwise if the
363  // nearest is e.g. on an edge in between inconsistent triangles it might
364  // make the wrong decision.
365  if (s.size() > 0)
366  {
367  // Whether face has to be flipped.
368  // UNVISITED: unvisited
369  // NOFLIP: no need to flip
370  // FLIP: need to flip
371  labelList flipState(s.size(), UNVISITED);
372 
373  label facei = 0;
374  while (true)
375  {
376  label startFacei = -1;
377  while (facei < s.size())
378  {
379  if (flipState[facei] == UNVISITED)
380  {
381  startFacei = facei;
382  break;
383  }
384  facei++;
385  }
386 
387  if (startFacei == -1)
388  {
389  break;
390  }
391 
392  flipState[startFacei] = NOFLIP;
393  walkSurface(s, startFacei, flipState);
394  }
395 
396  anyFlipped = flipSurface(s, flipState);
397  }
398  return anyFlipped;
399 }
400 
401 
402 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
403 
404 // Null constructor
406 :
407  triSurface()
408 {}
409 
410 
411 // Construct from surface and point which defines outside
413 (
414  const triSurface& surf,
415  const point& samplePoint,
416  const bool orientOutside
417 )
418 :
419  triSurface(surf)
420 {
421  orient(*this, samplePoint, orientOutside);
422 }
423 
424 
425 // Construct from surface. Calculate outside point.
427 (
428  const triSurface& surf,
429  const bool orientOutside
430 )
431 :
432  triSurface(surf)
433 {
434  // BoundBox calculation without localPoints
435  treeBoundBox bb(surf.points(), surf.meshPoints());
436 
437  point outsidePoint = bb.max() + bb.span();
438 
439  orient(*this, outsidePoint, orientOutside);
440 }
441 
442 
443 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
444 
446 (
447  triSurface& s,
448  const point& samplePoint,
449  const bool orientOutside
450 )
451 {
452  // Do initial flipping to make triangles consistent. Otherwise if the
453  // nearest is e.g. on an edge in between inconsistent triangles it might
454  // make the wrong decision.
455  bool topoFlipped = orientConsistent(s);
456 
457 
458  // Whether face has to be flipped.
459  // UNVISITED: unvisited
460  // NOFLIP: no need to flip
461  // FLIP: need to flip
462  labelList flipState(s.size(), UNVISITED);
463 
464 
465  while (true)
466  {
467  // Linear search for nearest unvisited point on surface.
468 
469  scalar minDist = great;
470  point minPoint;
471  label minFacei = -1;
472 
473  forAll(s, facei)
474  {
475  if (flipState[facei] == UNVISITED)
476  {
477  pointHit curHit =
478  s[facei].nearestPoint(samplePoint, s.points());
479 
480  if (curHit.distance() < minDist)
481  {
482  minDist = curHit.distance();
483  minPoint = curHit.rawPoint();
484  minFacei = facei;
485  }
486  }
487  }
488 
489  // Did we find anything?
490  if (minFacei == -1)
491  {
492  break;
493  }
494 
495  // From this nearest face see if needs to be flipped and then
496  // go outwards.
497  propagateOrientation
498  (
499  s,
500  samplePoint,
501  orientOutside,
502  minFacei,
503  minPoint,
504  flipState
505  );
506  }
507 
508  // Now finally flip triangles according to flipState.
509  bool geomFlipped = flipSurface(s, flipState);
510 
511  return topoFlipped || geomFlipped;
512 }
513 
514 
516 (
517  triSurface& s,
518  const triSurfaceSearch& querySurf,
519  const point& samplePoint,
520  const bool orientOutside
521 )
522 {
523  // Do initial flipping to make triangles consistent. Otherwise if the
524  // nearest is e.g. on an edge in between inconsistent triangles it might
525  // make the wrong decision.
526  bool topoFlipped = orientConsistent(s);
527 
528  // Determine disconnected parts of surface
529  boolList borderEdge(s.nEdges(), false);
530  forAll(s.edgeFaces(), edgeI)
531  {
532  if (s.edgeFaces()[edgeI].size() != 2)
533  {
534  borderEdge[edgeI] = true;
535  }
536  }
537  labelList faceZone;
538  label nZones = s.markZones(borderEdge, faceZone);
539 
540  // Check intersection with one face per zone.
541 
542  labelList flipState(s.size(), UNVISITED);
543  for (label zoneI = 0; zoneI < nZones; zoneI++)
544  {
545  label zoneFacei = -1;
546  bool isOutside;
547  findZoneSide
548  (
549  querySurf,
550  faceZone,
551  zoneI,
552  samplePoint,
553 
554  zoneFacei,
555  isOutside
556  );
557 
558  if (isOutside == orientOutside)
559  {
560  flipState[zoneFacei] = NOFLIP;
561  }
562  else
563  {
564  flipState[zoneFacei] = FLIP;
565  }
566  walkSurface(s, zoneFacei, flipState);
567  }
568 
569  // Now finally flip triangles according to flipState.
570  bool geomFlipped = flipSurface(s, flipState);
571 
572  return topoFlipped || geomFlipped;
573 }
574 
575 
576 // ************************************************************************* //
#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
static const Form max
Definition: VectorSpace.H:115
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool orient(triSurface &, const point &, const bool orientOutside)
Flip faces such that normals are consistent with point:
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static vector surfaceNormal(const triSurface &surf, const label nearestFacei, const point &nearestPt)
Triangle (unit) normal. If nearest point to triangle on edge use.
scalar minDist(const List< pointIndexHit > &hitList)
scalar f1
Definition: createFields.H:15
orientedSurface()
Construct null.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
Helper class to search on triSurface.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const Field< PointType > & points() const
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
label nEdges() const
Return number of edges in patch.
sideType
On which side of surface.
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
triSurface()
Construct null.
Definition: triSurface.C:534
std::remove_reference< ::Foam::List< labelledTri > >::type::value_type FaceType
static sideType surfaceSide(const triSurface &surf, const point &sample, const label nearestFacei)
Given nearest point (to sample) on surface determines which side.
vector point
Point is a vector.
Definition: point.H:41
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
label markZones(const boolList &borderEdge, labelList &faceZone) const
(size and) fills faceZone with zone of face. Zone is area
Definition: triSurface.C:999
Triangulated surface description with patch information.
Definition: triSurface.H:66
Namespace for OpenFOAM.