This repository has been archived on 2023-11-20. You can view files and clone it, but cannot push or open issues or pull requests.
foam-extend4.1-coherent-io/applications/utilities/mesh/manipulation/splitMesh/regionSide.C
2016-06-21 15:04:12 +02:00

479 lines
12 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "regionSide.H"
#include "meshTools.H"
#include "primitiveMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(regionSide, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Step across edge onto other face on cell
Foam::label Foam::regionSide::otherFace
(
const primitiveMesh& mesh,
const label cellI,
const label faceI,
const label edgeI
)
{
label f0I, f1I;
meshTools::getEdgeFaces(mesh, cellI, edgeI, f0I, f1I);
if (f0I == faceI)
{
return f1I;
}
else
{
return f0I;
}
}
// Step across point to other edge on face
Foam::label Foam::regionSide::otherEdge
(
const primitiveMesh& mesh,
const label faceI,
const label edgeI,
const label pointI
)
{
const edge& e = mesh.edges()[edgeI];
// Get other point on edge.
label freePointI = e.otherVertex(pointI);
const labelList& fEdges = mesh.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
label otherEdgeI = fEdges[fEdgeI];
const edge& otherE = mesh.edges()[otherEdgeI];
if
(
(
otherE.start() == pointI
&& otherE.end() != freePointI
)
|| (
otherE.end() == pointI
&& otherE.start() != freePointI
)
)
{
// otherE shares one (but not two) points with e.
return otherEdgeI;
}
}
FatalErrorIn
(
"regionSide::otherEdge(const primitiveMesh&, const label, const label"
", const label)"
) << "Cannot find other edge on face " << faceI << " that uses point "
<< pointI << " but not point " << freePointI << endl
<< "Edges on face:" << fEdges
<< " verts:" << UIndirectList<edge>(mesh.edges(), fEdges)()
<< " Vertices on face:"
<< mesh.faces()[faceI]
<< " Vertices on original edge:" << e << abort(FatalError);
return -1;
}
// Step from faceI (on side cellI) to connected face & cell without crossing
// fenceEdges.
void Foam::regionSide::visitConnectedFaces
(
const primitiveMesh& mesh,
const labelHashSet& region,
const labelHashSet& fenceEdges,
const label cellI,
const label faceI,
labelHashSet& visitedFace
)
{
if (!visitedFace.found(faceI))
{
if (debug)
{
Info<< "visitConnectedFaces : cellI:" << cellI << " faceI:"
<< faceI << " isOwner:" << (cellI == mesh.faceOwner()[faceI])
<< endl;
}
// Mark as visited
visitedFace.insert(faceI);
// Mark which side of face was visited.
if (cellI == mesh.faceOwner()[faceI])
{
sideOwner_.insert(faceI);
}
// Visit all neighbouring faces on faceSet. Stay on this 'side' of
// face by doing edge-face-cell walk.
const labelList& fEdges = mesh.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
label edgeI = fEdges[fEdgeI];
if (!fenceEdges.found(edgeI))
{
// Step along faces on edge from cell to cell until
// we hit face on faceSet.
// Find face reachable from edge
label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
if (mesh.isInternalFace(otherFaceI))
{
label otherCellI = cellI;
// Keep on crossing faces/cells until back on face on
// surface
while (!region.found(otherFaceI))
{
visitedFace.insert(otherFaceI);
if (debug)
{
Info<< "visitConnectedFaces : cellI:" << cellI
<< " found insideEdgeFace:" << otherFaceI
<< endl;
}
// Cross otherFaceI into neighbouring cell
otherCellI =
meshTools::otherCell
(
mesh,
otherCellI,
otherFaceI
);
otherFaceI =
otherFace
(
mesh,
otherCellI,
otherFaceI,
edgeI
);
}
visitConnectedFaces
(
mesh,
region,
fenceEdges,
otherCellI,
otherFaceI,
visitedFace
);
}
}
}
}
}
// From edge on face connected to point on region (regionPointI) cross
// to all other edges using this point by walking across faces
// Does not cross regionEdges so stays on one side
// of region
void Foam::regionSide::walkPointConnectedFaces
(
const primitiveMesh& mesh,
const labelHashSet& regionEdges,
const label regionPointI,
const label startFaceI,
const label startEdgeI,
labelHashSet& visitedEdges
)
{
// Mark as visited
insidePointFaces_.insert(startFaceI);
if (debug)
{
Info<< "walkPointConnectedFaces : regionPointI:" << regionPointI
<< " faceI:" << startFaceI
<< " edgeI:" << startEdgeI << " verts:"
<< mesh.edges()[startEdgeI]
<< endl;
}
// Cross faceI i.e. get edge not startEdgeI which uses regionPointI
label edgeI = otherEdge(mesh, startFaceI, startEdgeI, regionPointI);
if (!regionEdges.found(edgeI))
{
if (!visitedEdges.found(edgeI))
{
visitedEdges.insert(edgeI);
if (debug)
{
Info<< "Crossed face from "
<< " edgeI:" << startEdgeI << " verts:"
<< mesh.edges()[startEdgeI]
<< " to edge:" << edgeI << " verts:"
<< mesh.edges()[edgeI]
<< endl;
}
// Cross edge to all faces connected to it.
const labelList& eFaces = mesh.edgeFaces()[edgeI];
forAll(eFaces, eFaceI)
{
label faceI = eFaces[eFaceI];
walkPointConnectedFaces
(
mesh,
regionEdges,
regionPointI,
faceI,
edgeI,
visitedEdges
);
}
}
}
}
// Find all faces reachable from all non-fence points and staying on
// regionFaces side.
void Foam::regionSide::walkAllPointConnectedFaces
(
const primitiveMesh& mesh,
const labelHashSet& regionFaces,
const labelHashSet& fencePoints
)
{
//
// Get all (internal and external) edges on region.
//
labelHashSet regionEdges(4*regionFaces.size());
for
(
labelHashSet::const_iterator iter = regionFaces.begin();
iter != regionFaces.end();
++iter
)
{
label faceI = iter.key();
const labelList& fEdges = mesh.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
regionEdges.insert(fEdges[fEdgeI]);
}
}
//
// Visit all internal points on surface.
//
// Storage for visited points
labelHashSet visitedPoint(4*regionFaces.size());
// Insert fence points so we don't visit them
for
(
labelHashSet::const_iterator iter = fencePoints.begin();
iter != fencePoints.end();
++iter
)
{
visitedPoint.insert(iter.key());
}
labelHashSet visitedEdges(2*fencePoints.size());
if (debug)
{
Info<< "Excluding visit of points:" << visitedPoint << endl;
}
for
(
labelHashSet::const_iterator iter = regionFaces.begin();
iter != regionFaces.end();
++iter
)
{
label faceI = iter.key();
// Get side of face.
label cellI;
if (sideOwner_.found(faceI))
{
cellI = mesh.faceOwner()[faceI];
}
else
{
cellI = mesh.faceNeighbour()[faceI];
}
// Find starting point and edge on face.
const labelList& fEdges = mesh.faceEdges()[faceI];
forAll(fEdges, fEdgeI)
{
label edgeI = fEdges[fEdgeI];
// Get the face 'perpendicular' to faceI on region.
label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
// Edge
const edge& e = mesh.edges()[edgeI];
if (!visitedPoint.found(e.start()))
{
Info<< "Determining visibility from point " << e.start()
<< endl;
visitedPoint.insert(e.start());
//edgeI = otherEdge(mesh, otherFaceI, edgeI, e.start());
walkPointConnectedFaces
(
mesh,
regionEdges,
e.start(),
otherFaceI,
edgeI,
visitedEdges
);
}
if (!visitedPoint.found(e.end()))
{
Info<< "Determining visibility from point " << e.end()
<< endl;
visitedPoint.insert(e.end());
//edgeI = otherEdge(mesh, otherFaceI, edgeI, e.end());
walkPointConnectedFaces
(
mesh,
regionEdges,
e.end(),
otherFaceI,
edgeI,
visitedEdges
);
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::regionSide::regionSide
(
const primitiveMesh& mesh,
const labelHashSet& region, // faces of region
const labelHashSet& fenceEdges, // outside edges
const label startCellI,
const label startFaceI
)
:
sideOwner_(region.size()),
insidePointFaces_(region.size())
{
// Storage for visited faces
labelHashSet visitedFace(region.size());
// Visit all faces on this side of region.
// Mark for each face whether owner (or neighbour) has been visited.
// Sets sideOwner_
visitConnectedFaces
(
mesh,
region,
fenceEdges,
startCellI,
startFaceI,
visitedFace
);
//
// Visit all internal points on region and mark faces visitable from these.
// Sets insidePointFaces_.
//
labelHashSet fencePoints(fenceEdges.size());
for
(
labelHashSet::const_iterator iter = fenceEdges.begin();
iter != fenceEdges.end();
++iter
)
{
const edge& e = mesh.edges()[iter.key()];
fencePoints.insert(e.start());
fencePoints.insert(e.end());
}
walkAllPointConnectedFaces(mesh, region, fencePoints);
}
// ************************************************************************* //