src/engine/engineTopoChangerMesh/regionSide/regionSide.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "regionSide.H"
  29 #include "meshTools.H"
  30 #include "primitiveMesh.H"
  31 #include "IndirectList.H"
  32
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 defineTypeNameAndDebug(regionSide, 0);
  40
  41 }
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 // Step across edge onto other face on cell
  46 Foam::label Foam::regionSide::otherFace
  47 (
  48     const primitiveMesh& mesh,
  49     const label cellI,
  50     const label faceI,
  51     const label edgeI
  52 )
  53 {
  54     label f0I, f1I;
  55
  56     meshTools::getEdgeFaces(mesh, cellI, edgeI, f0I, f1I);
  57
  58     if (f0I == faceI)
  59     {
  60         return f1I;
  61     }
  62     else
  63     {
  64         return f0I;
  65     }
  66 }
  67
  68
  69 // Step across point to other edge on face
  70 Foam::label Foam::regionSide::otherEdge
  71 (
  72     const primitiveMesh& mesh,
  73     const label faceI,
  74     const label edgeI,
  75     const label pointI
  76 )
  77 {
  78     const edge& e = mesh.edges()[edgeI];
  79
  80     // Get other point on edge.
  81     label freePointI = e.otherVertex(pointI);
  82
  83     const labelList& fEdges = mesh.faceEdges()[faceI];
  84
  85     forAll(fEdges, fEdgeI)
  86     {
  87         label otherEdgeI = fEdges[fEdgeI];
  88
  89         const edge& otherE = mesh.edges()[otherEdgeI];
  90
  91         if
  92         (
  93             (
  94                 otherE.start() == pointI
  95              && otherE.end() != freePointI
  96             )
  97          || (
  98                 otherE.end() == pointI
  99              && otherE.start() != freePointI
 100             )
 101         )
 102         {
 103             // otherE shares one (but not two) points with e.
 104             return otherEdgeI;
 105         }
 106     }
 107
 108     FatalErrorIn
 109     (
 110         "regionSide::otherEdge(const primitiveMesh&, const label, const label"
 111         ", const label)"
 112     )   << "Cannot find other edge on face " << faceI << " that uses point "
 113         << pointI << " but not point " << freePointI << endl
 114         << "Edges on face:" << fEdges
 115         << " verts:" << IndirectList<edge>(mesh.edges(), fEdges)()
 116         << " Vertices on face:"
 117         << mesh.faces()[faceI]
 118         << " Vertices on original edge:" << e << abort(FatalError);
 119
 120     return -1;
 121 }
 122
 123
 124 // Step from faceI (on side cellI) to connected face & cell without crossing
 125 // fenceEdges.
 126 void Foam::regionSide::visitConnectedFaces
 127 (
 128     const primitiveMesh& mesh,
 129     const labelHashSet& region,
 130     const labelHashSet& fenceEdges,
 131     const label cellI,
 132     const label faceI,
 133     labelHashSet& visitedFace
 134 )
 135 {
 136     if (!visitedFace.found(faceI))
 137     {
 138         if (debug)
 139         {
 140             Info<< "visitConnectedFaces : cellI:" << cellI << " faceI:"
 141                 << faceI << "  isOwner:" << (cellI == mesh.faceOwner()[faceI])
 142                 << endl;
 143         }
 144
 145         // Mark as visited
 146         visitedFace.insert(faceI);
 147
 148         // Mark which side of face was visited.
 149         if (cellI == mesh.faceOwner()[faceI])
 150         {
 151             sideOwner_.insert(faceI);
 152         }
 153
 154
 155         // Visit all neighbouring faces on faceSet. Stay on this 'side' of
 156         // face by doing edge-face-cell walk.
 157         const labelList& fEdges = mesh.faceEdges()[faceI];
 158
 159         forAll(fEdges, fEdgeI)
 160         {
 161             label edgeI = fEdges[fEdgeI];
 162
 163             if (!fenceEdges.found(edgeI))
 164             {
 165                 // Step along faces on edge from cell to cell until
 166                 // we hit face on faceSet.
 167
 168                 // Find face reachable from edge
 169                 label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
 170
 171                 if (mesh.isInternalFace(otherFaceI))
 172                 {
 173                     label otherCellI = cellI;
 174
 175                     // Keep on crossing faces/cells until back on face on
 176                     // surface
 177                     while (!region.found(otherFaceI))
 178                     {
 179                         visitedFace.insert(otherFaceI);
 180
 181                         if (debug)
 182                         {
 183                             Info<< "visitConnectedFaces : cellI:" << cellI
 184                                 << " found insideEdgeFace:" << otherFaceI
 185                                 << endl;
 186                         }
 187
 188
 189                         // Cross otherFaceI into neighbouring cell
 190                         otherCellI =
 191                             meshTools::otherCell
 192                             (
 193                                 mesh,
 194                                 otherCellI,
 195                                 otherFaceI
 196                             );
 197
 198                         otherFaceI =
 199                                 otherFace
 200                                 (
 201                                     mesh,
 202                                     otherCellI,
 203                                     otherFaceI,
 204                                     edgeI
 205                                 );
 206                     }
 207
 208                     visitConnectedFaces
 209                     (
 210                         mesh,
 211                         region,
 212                         fenceEdges,
 213                         otherCellI,
 214                         otherFaceI,
 215                         visitedFace
 216                     );
 217                 }
 218             }
 219         }
 220     }
 221 }
 222
 223
 224 // From edge on face connected to point on region (regionPointI) cross
 225 // to all other edges using this point by walking across faces
 226 // Does not cross regionEdges so stays on one side
 227 // of region
 228 void Foam::regionSide::walkPointConnectedFaces
 229 (
 230     const primitiveMesh& mesh,
 231     const labelHashSet& regionEdges,
 232     const label regionPointI,
 233     const label startFaceI,
 234     const label startEdgeI,
 235     labelHashSet& visitedEdges
 236 )
 237 {
 238     // Mark as visited
 239     insidePointFaces_.insert(startFaceI);
 240
 241     if (debug)
 242     {
 243         Info<< "walkPointConnectedFaces : regionPointI:" << regionPointI
 244             << " faceI:" << startFaceI
 245             << " edgeI:" << startEdgeI << " verts:"
 246             << mesh.edges()[startEdgeI]
 247             << endl;
 248     }
 249
 250     // Cross faceI i.e. get edge not startEdgeI which uses regionPointI
 251     label edgeI = otherEdge(mesh, startFaceI, startEdgeI, regionPointI);
 252
 253     if (!regionEdges.found(edgeI))
 254     {
 255         if (!visitedEdges.found(edgeI))
 256         {
 257             visitedEdges.insert(edgeI);
 258
 259             if (debug)
 260             {
 261                 Info<< "Crossed face from "
 262                     << " edgeI:" << startEdgeI << " verts:"
 263                     << mesh.edges()[startEdgeI]
 264                     << " to edge:" << edgeI << " verts:"
 265                     << mesh.edges()[edgeI]
 266                     << endl;
 267             }
 268
 269             // Cross edge to all faces connected to it.
 270
 271             const labelList& eFaces = mesh.edgeFaces()[edgeI];
 272
 273             forAll(eFaces, eFaceI)
 274             {
 275                 label faceI = eFaces[eFaceI];
 276
 277                 walkPointConnectedFaces
 278                 (
 279                     mesh,
 280                     regionEdges,
 281                     regionPointI,
 282                     faceI,
 283                     edgeI,
 284                     visitedEdges
 285                 );
 286             }
 287         }
 288     }
 289 }
 290
 291
 292 // Find all faces reachable from all non-fence points and staying on
 293 // regionFaces side.
 294 void Foam::regionSide::walkAllPointConnectedFaces
 295 (
 296     const primitiveMesh& mesh,
 297     const labelHashSet& regionFaces,
 298     const labelHashSet& fencePoints
 299 )
 300 {
 301     //
 302     // Get all (internal and external) edges on region.
 303     //
 304     labelHashSet regionEdges(4*regionFaces.size());
 305
 306     for
 307     (
 308         labelHashSet::const_iterator iter = regionFaces.begin();
 309         iter != regionFaces.end();
 310         ++iter
 311     )
 312     {
 313         label faceI = iter.key();
 314
 315         const labelList& fEdges = mesh.faceEdges()[faceI];
 316
 317         forAll(fEdges, fEdgeI)
 318         {
 319             regionEdges.insert(fEdges[fEdgeI]);
 320         }
 321     }
 322
 323
 324     //
 325     // Visit all internal points on surface.
 326     //
 327
 328     // Storage for visited points
 329     labelHashSet visitedPoint(4*regionFaces.size());
 330
 331     // Insert fence points so we don't visit them
 332     for
 333     (
 334         labelHashSet::const_iterator iter = fencePoints.begin();
 335         iter != fencePoints.end();
 336         ++iter
 337     )
 338     {
 339         visitedPoint.insert(iter.key());
 340     }
 341
 342     labelHashSet visitedEdges(2*fencePoints.size());
 343
 344
 345     if (debug)
 346     {
 347         Info<< "Excluding visit of points:" << visitedPoint << endl;
 348     }
 349
 350     for
 351     (
 352         labelHashSet::const_iterator iter = regionFaces.begin();
 353         iter != regionFaces.end();
 354         ++iter
 355     )
 356     {
 357         label faceI = iter.key();
 358
 359         // Get side of face.
 360         label cellI;
 361
 362         if (sideOwner_.found(faceI))
 363         {
 364             cellI = mesh.faceOwner()[faceI];
 365         }
 366         else
 367         {
 368             cellI = mesh.faceNeighbour()[faceI];
 369         }
 370
 371         // Find starting point and edge on face.
 372         const labelList& fEdges = mesh.faceEdges()[faceI];
 373
 374         forAll(fEdges, fEdgeI)
 375         {
 376             label edgeI = fEdges[fEdgeI];
 377
 378             // Get the face 'perpendicular' to faceI on region.
 379             label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
 380
 381             // Edge
 382             const edge& e = mesh.edges()[edgeI];
 383
 384             if (!visitedPoint.found(e.start()))
 385             {
 386                 Info<< "Determining visibility from point " << e.start()
 387                     << endl;
 388
 389                 visitedPoint.insert(e.start());
 390
 391                 //edgeI = otherEdge(mesh, otherFaceI, edgeI, e.start());
 392
 393                 walkPointConnectedFaces
 394                 (
 395                     mesh,
 396                     regionEdges,
 397                     e.start(),
 398                     otherFaceI,
 399                     edgeI,
 400                     visitedEdges
 401                 );
 402             }
 403             if (!visitedPoint.found(e.end()))
 404             {
 405                 Info<< "Determining visibility from point " << e.end()
 406                     << endl;
 407
 408                 visitedPoint.insert(e.end());
 409
 410                 //edgeI = otherEdge(mesh, otherFaceI, edgeI, e.end());
 411
 412                 walkPointConnectedFaces
 413                 (
 414                     mesh,
 415                     regionEdges,
 416                     e.end(),
 417                     otherFaceI,
 418                     edgeI,
 419                     visitedEdges
 420                 );
 421             }
 422         }
 423     }
 424 }
 425
 426
 427 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 428
 429 // Construct from components
 430 Foam::regionSide::regionSide
 431 (
 432     const primitiveMesh& mesh,
 433     const labelHashSet& region,         // faces of region
 434     const labelHashSet& fenceEdges,     // outside edges
 435     const label startCellI,
 436     const label startFaceI
 437 )
 438 :
 439     sideOwner_(region.size()),
 440     insidePointFaces_(region.size())
 441 {
 442     // Storage for visited faces
 443     labelHashSet visitedFace(region.size());
 444
 445     // Visit all faces on this side of region.
 446     // Mark for each face whether owner (or neighbour) has been visited.
 447     // Sets sideOwner_
 448     visitConnectedFaces
 449     (
 450         mesh,
 451         region,
 452         fenceEdges,
 453         startCellI,
 454         startFaceI,
 455         visitedFace
 456     );
 457
 458     //
 459     // Visit all internal points on region and mark faces visitable from these.
 460     // Sets insidePointFaces_.
 461     //
 462
 463     labelHashSet fencePoints(fenceEdges.size());
 464
 465     for
 466     (
 467         labelHashSet::const_iterator iter = fenceEdges.begin();
 468         iter != fenceEdges.end();
 469         ++iter
 470     )
 471     {
 472         const edge& e = mesh.edges()[iter.key()];
 473
 474         fencePoints.insert(e.start());
 475         fencePoints.insert(e.end());
 476     }
 477
 478     walkAllPointConnectedFaces(mesh, region, fencePoints);
 479 }
 480
 481
 482 // ************************************************************************* //