FreeFOAM The Cross-Platform CFD Toolkit
orientedSurface.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 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "orientedSurface.H"
30 #include <meshTools/treeBoundBox.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
35 
36 
37 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
38 
39 // Return true if face uses edge from start to end.
40 bool Foam::orientedSurface::edgeOrder
41 (
42  const labelledTri& f,
43  const edge& e
44 )
45 {
46  if
47  (
48  (f[0] == e[0] && f[1] == e[1])
49  || (f[1] == e[0] && f[2] == e[1])
50  || (f[2] == e[0] && f[0] == e[1])
51  )
52  {
53  return true;
54  }
55  else
56  {
57  return false;
58  }
59 }
60 
61 
62 // Return true if edge is used in opposite order in faces
63 bool Foam::orientedSurface::consistentEdge
64 (
65  const edge& e,
66  const labelledTri& f0,
67  const labelledTri& f1
68 )
69 {
70  return edgeOrder(f0, e) ^ edgeOrder(f1, e);
71 }
72 
73 
74 Foam::labelList Foam::orientedSurface::faceToEdge
75 (
76  const triSurface& s,
77  const labelList& changedFaces
78 )
79 {
80  labelList changedEdges(3*changedFaces.size());
81  label changedI = 0;
82 
83  forAll(changedFaces, i)
84  {
85  const labelList& fEdges = s.faceEdges()[changedFaces[i]];
86 
87  forAll(fEdges, j)
88  {
89  changedEdges[changedI++] = fEdges[j];
90  }
91  }
92  changedEdges.setSize(changedI);
93 
94  return changedEdges;
95 }
96 
97 
98 Foam::labelList Foam::orientedSurface::edgeToFace
99 (
100  const triSurface& s,
101  const labelList& changedEdges,
102  labelList& flip
103 )
104 {
105  labelList changedFaces(2*changedEdges.size());
106  label changedI = 0;
107 
108  forAll(changedEdges, i)
109  {
110  label edgeI = changedEdges[i];
111 
112  const labelList& eFaces = s.edgeFaces()[edgeI];
113 
114  if (eFaces.size() < 2)
115  {
116  // Do nothing, faces was already visited.
117  }
118  else if (eFaces.size() == 2)
119  {
120  label face0 = eFaces[0];
121  label face1 = eFaces[1];
122 
123  const labelledTri& f0 = s.localFaces()[face0];
124  const labelledTri& f1 = s.localFaces()[face1];
125 
126  if (flip[face0] == UNVISITED)
127  {
128  if (flip[face1] == UNVISITED)
129  {
130  FatalErrorIn("orientedSurface::edgeToFace") << "Problem"
131  << abort(FatalError);
132  }
133  else
134  {
135  // Face1 has a flip state, face0 hasn't
136  if (consistentEdge(s.edges()[edgeI], f0, f1))
137  {
138  // Take over flip status
139  flip[face0] = (flip[face1] == FLIP ? FLIP : NOFLIP);
140  }
141  else
142  {
143  // Invert
144  flip[face0] = (flip[face1] == FLIP ? NOFLIP : FLIP);
145  }
146  changedFaces[changedI++] = face0;
147  }
148  }
149  else
150  {
151  if (flip[face1] == UNVISITED)
152  {
153  // Face0 has a flip state, face1 hasn't
154  if (consistentEdge(s.edges()[edgeI], f0, f1))
155  {
156  flip[face1] = (flip[face0] == FLIP ? FLIP : NOFLIP);
157  }
158  else
159  {
160  flip[face1] = (flip[face0] == FLIP ? NOFLIP : FLIP);
161  }
162  changedFaces[changedI++] = face1;
163  }
164  }
165  }
166  else
167  {
168  // Multiply connected. Do what?
169  }
170  }
171  changedFaces.setSize(changedI);
172 
173  return changedFaces;
174 }
175 
176 
177 void Foam::orientedSurface::walkSurface
178 (
179  const triSurface& s,
180  const label startFaceI,
181  labelList& flipState
182 )
183 {
184  // List of faces that were changed in the last iteration.
185  labelList changedFaces(1, startFaceI);
186  // List of edges that were changed in the last iteration.
187  labelList changedEdges;
188 
189  while(true)
190  {
191  changedEdges = faceToEdge(s, changedFaces);
192 
193  if (debug)
194  {
195  Pout<< "From changedFaces:" << changedFaces.size()
196  << " to changedEdges:" << changedEdges.size()
197  << endl;
198  }
199 
200  if (changedEdges.empty())
201  {
202  break;
203  }
204 
205  changedFaces = edgeToFace(s, changedEdges, flipState);
206 
207  if (debug)
208  {
209  Pout<< "From changedEdges:" << changedEdges.size()
210  << " to changedFaces:" << changedFaces.size()
211  << endl;
212  }
213 
214  if (changedFaces.empty())
215  {
216  break;
217  }
218  }
219 }
220 
221 
222 void Foam::orientedSurface::propagateOrientation
223 (
224  const triSurface& s,
225  const point& samplePoint,
226  const bool orientOutside,
227  const label nearestFaceI,
228  const point& nearestPt,
229  labelList& flipState
230 )
231 {
232  //
233  // Determine orientation to normal on nearest face
234  //
236  (
237  s,
238  samplePoint,
239  nearestFaceI,
240  nearestPt,
241  10*SMALL
242  );
243 
244  if (side == triSurfaceTools::UNKNOWN)
245  {
246  // Non-closed surface. Do what? For now behave as if no flipping
247  // nessecary
248  flipState[nearestFaceI] = NOFLIP;
249  }
250  else if ((side == triSurfaceTools::OUTSIDE) == orientOutside)
251  {
252  // outside & orientOutside or inside & !orientOutside
253  // Normals on surface pointing correctly. No need to flip normals
254  flipState[nearestFaceI] = NOFLIP;
255  }
256  else
257  {
258  // Need to flip normals.
259  flipState[nearestFaceI] = FLIP;
260  }
261 
262  if (debug)
263  {
264  vector n = triSurfaceTools::surfaceNormal(s, nearestFaceI, nearestPt);
265 
266  Pout<< "orientedSurface::propagateOrientation : starting face"
267  << " orientation:" << nl
268  << " for samplePoint:" << samplePoint << nl
269  << " starting from point:" << nearestPt << nl
270  << " on face:" << nearestFaceI << nl
271  << " with normal:" << n << nl
272  << " decided side:" << label(side)
273  << endl;
274  }
275 
276  // Walk the surface from nearestFaceI, changing the flipstate.
277  walkSurface(s, nearestFaceI, flipState);
278 }
279 
280 
281 bool Foam::orientedSurface::flipSurface
282 (
283  triSurface& s,
284  const labelList& flipState
285 )
286 {
287  bool hasFlipped = false;
288 
289  // Flip tris in s
290  forAll(flipState, faceI)
291  {
292  if (flipState[faceI] == UNVISITED)
293  {
295  (
296  "orientSurface(const point&, const label, const point&)"
297  ) << "unvisited face " << faceI
298  << abort(FatalError);
299  }
300  else if (flipState[faceI] == FLIP)
301  {
302  labelledTri& tri = s[faceI];
303 
304  label tmp = tri[0];
305 
306  tri[0] = tri[1];
307  tri[1] = tmp;
308 
309  hasFlipped = true;
310  }
311  }
312  // Recalculate normals
313  s.clearOut();
314 
315  return hasFlipped;
316 }
317 
318 
319 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
320 
321 // Null constructor
323 :
324  triSurface()
325 {}
326 
327 
328 // Construct from surface and point which defines outside
330 (
331  const triSurface& surf,
332  const point& samplePoint,
333  const bool orientOutside
334 )
335 :
336  triSurface(surf)
337 {
338  orient(*this, samplePoint, orientOutside);
339 }
340 
341 
342 // Construct from surface. Calculate outside point.
344 (
345  const triSurface& surf,
346  const bool orientOutside
347 )
348 :
349  triSurface(surf)
350 {
351  // BoundBox calculation without localPoints
352  treeBoundBox bb(surf.points(), surf.meshPoints());
353 
354  point outsidePoint = bb.max() + bb.span();
355 
356  orient(*this, outsidePoint, orientOutside);
357 }
358 
359 
360 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
361 
363 (
364  triSurface& s,
365  const point& samplePoint,
366  const bool orientOutside
367 )
368 {
369  bool anyFlipped = false;
370 
371  // Do initial flipping to make triangles consistent. Otherwise if the
372  // nearest is e.g. on an edge inbetween inconsistent triangles it might
373  // make the wrong decision.
374  if (s.size() > 0)
375  {
376  // Whether face has to be flipped.
377  // UNVISITED: unvisited
378  // NOFLIP: no need to flip
379  // FLIP: need to flip
380  labelList flipState(s.size(), UNVISITED);
381 
382  flipState[0] = NOFLIP;
383  walkSurface(s, 0, flipState);
384 
385  anyFlipped = flipSurface(s, flipState);
386  }
387 
388 
389  // Whether face has to be flipped.
390  // UNVISITED: unvisited
391  // NOFLIP: no need to flip
392  // FLIP: need to flip
393  labelList flipState(s.size(), UNVISITED);
394 
395 
396  while (true)
397  {
398  // Linear search for nearest unvisited point on surface.
399 
400  scalar minDist = GREAT;
401  point minPoint;
402  label minFaceI = -1;
403 
404  forAll(s, faceI)
405  {
406  if (flipState[faceI] == UNVISITED)
407  {
408  const labelledTri& f = s[faceI];
409 
410  pointHit curHit =
412  (
413  s.points()[f[0]],
414  s.points()[f[1]],
415  s.points()[f[2]]
416  ).nearestPoint(samplePoint);
417 
418  if (curHit.distance() < minDist)
419  {
420  minDist = curHit.distance();
421  minPoint = curHit.rawPoint();
422  minFaceI = faceI;
423  }
424  }
425  }
426 
427  // Did we find anything?
428  if (minFaceI == -1)
429  {
430  break;
431  }
432 
433  // From this nearest face see if needs to be flipped and then
434  // go outwards.
435  propagateOrientation
436  (
437  s,
438  samplePoint,
439  orientOutside,
440  minFaceI,
441  minPoint,
442  flipState
443  );
444  }
445 
446  // Now finally flip triangles according to flipState.
447  bool geomFlipped = flipSurface(s, flipState);
448 
449  return anyFlipped || geomFlipped;
450 }
451 
452 
453 // ************************ vim: set sw=4 sts=4 et: ************************ //