applications/utilities/mesh/manipulation/splitMesh/regionSide.C

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