src/meshTools/triSurface/orientedSurface/orientedSurface.C

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