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