mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
3257 lines
80 KiB
C
3257 lines
80 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011-2017 OpenFOAM Foundation
|
|
Copyright (C) 2017-2019 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM 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.
|
|
|
|
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "cellCuts.H"
|
|
#include "polyMesh.H"
|
|
#include "Time.H"
|
|
#include "ListOps.H"
|
|
#include "cellLooper.H"
|
|
#include "refineCell.H"
|
|
#include "meshTools.H"
|
|
#include "geomCellLooper.H"
|
|
#include "OFstream.H"
|
|
#include "plane.H"
|
|
#include "syncTools.H"
|
|
#include "dummyTransform.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
defineTypeNameAndDebug(cellCuts, 0);
|
|
|
|
//- Template specialization for pTraits<edge> so we can use syncTools
|
|
// functionality
|
|
template<>
|
|
class pTraits<edge>
|
|
{
|
|
public:
|
|
|
|
//- Component type
|
|
typedef edge cmptType;
|
|
};
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Static Functions * * * * * * * * * * * //
|
|
|
|
Foam::label Foam::cellCuts::findPartIndex
|
|
(
|
|
const labelList& elems,
|
|
const label nElems,
|
|
const label val
|
|
)
|
|
{
|
|
for (label i = 0; i < nElems; i++)
|
|
{
|
|
if (elems[i] == val)
|
|
{
|
|
return i;
|
|
}
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
|
|
Foam::boolList Foam::cellCuts::expand
|
|
(
|
|
const label size,
|
|
const labelList& labels
|
|
)
|
|
{
|
|
boolList result(size, false);
|
|
|
|
forAll(labels, labelI)
|
|
{
|
|
result[labels[labelI]] = true;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
Foam::scalarField Foam::cellCuts::expand
|
|
(
|
|
const label size,
|
|
const labelList& labels,
|
|
const scalarField& weights
|
|
)
|
|
{
|
|
scalarField result(size, -GREAT);
|
|
|
|
forAll(labels, labelI)
|
|
{
|
|
result[labels[labelI]] = weights[labelI];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::firstUnique
|
|
(
|
|
const labelList& lst,
|
|
const Map<label>& map
|
|
)
|
|
{
|
|
forAll(lst, i)
|
|
{
|
|
if (!map.found(lst[i]))
|
|
{
|
|
return i;
|
|
}
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
void Foam::cellCuts::syncProc()
|
|
{
|
|
if (!Pstream::parRun())
|
|
{
|
|
return;
|
|
}
|
|
|
|
syncTools::syncPointList(mesh(), pointIsCut_, orEqOp<bool>(), false);
|
|
syncTools::syncEdgeList(mesh(), edgeIsCut_, orEqOp<bool>(), false);
|
|
syncTools::syncEdgeList(mesh(), edgeWeight_, maxEqOp<scalar>(), -GREAT);
|
|
|
|
{
|
|
const label nBnd = mesh().nBoundaryFaces();
|
|
|
|
// Convert faceSplitCut into face-local data: vertex and edge w.r.t.
|
|
// vertex 0: (since this is same on both sides)
|
|
//
|
|
// Sending-side vertex Receiving-side vertex
|
|
// 0 0
|
|
// 1 3
|
|
// 2 2
|
|
// 3 1
|
|
//
|
|
// Sending-side edge Receiving side edge
|
|
// 0-1 3-0
|
|
// 1-2 2-3
|
|
// 2-3 1-2
|
|
// 3-0 0-1
|
|
//
|
|
// Encoding is as index:
|
|
// 0 : not set
|
|
// >0 : vertex, vertex is index-1
|
|
// <0 : edge, edge is -index-1
|
|
|
|
edgeList relCuts(nBnd, edge(0, 0));
|
|
|
|
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
|
|
|
|
forAll(pbm, patchi)
|
|
{
|
|
const polyPatch& pp = pbm[patchi];
|
|
|
|
if (pp.coupled())
|
|
{
|
|
forAll(pp, i)
|
|
{
|
|
label facei = pp.start()+i;
|
|
label bFacei = facei-mesh().nInternalFaces();
|
|
|
|
const auto iter = faceSplitCut_.cfind(facei);
|
|
|
|
if (iter.found())
|
|
{
|
|
const face& f = mesh().faces()[facei];
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
const edge& cuts = iter.val();
|
|
|
|
forAll(cuts, i)
|
|
{
|
|
if (isEdge(cuts[i]))
|
|
{
|
|
label edgei = getEdge(cuts[i]);
|
|
label index = fEdges.find(edgei);
|
|
relCuts[bFacei][i] = -index-1;
|
|
}
|
|
else
|
|
{
|
|
label index = f.find(getVertex(cuts[i]));
|
|
relCuts[bFacei][i] = index+1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Exchange
|
|
syncTools::syncBoundaryFaceList
|
|
(
|
|
mesh(),
|
|
relCuts,
|
|
eqOp<edge>(),
|
|
dummyTransform()
|
|
);
|
|
|
|
// Convert relCuts back into mesh based data
|
|
forAll(pbm, patchi)
|
|
{
|
|
const polyPatch& pp = pbm[patchi];
|
|
|
|
if (pp.coupled())
|
|
{
|
|
forAll(pp, i)
|
|
{
|
|
label facei = pp.start()+i;
|
|
label bFacei = facei-mesh().nInternalFaces();
|
|
|
|
const edge& relCut = relCuts[bFacei];
|
|
if (relCut != edge(0, 0))
|
|
{
|
|
const face& f = mesh().faces()[facei];
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
|
|
edge absoluteCut(0, 0);
|
|
forAll(relCut, i)
|
|
{
|
|
if (relCut[i] < 0)
|
|
{
|
|
label oppFp = -relCut[i]-1;
|
|
label fp = f.size()-1-oppFp;
|
|
absoluteCut[i] = edgeToEVert(fEdges[fp]);
|
|
}
|
|
else
|
|
{
|
|
label oppFp = relCut[i]-1;
|
|
label fp =
|
|
(
|
|
oppFp == 0
|
|
? 0
|
|
: f.size()-oppFp
|
|
);
|
|
absoluteCut[i] = vertToEVert(f[fp]);
|
|
}
|
|
}
|
|
|
|
if
|
|
(
|
|
!faceSplitCut_.insert(facei, absoluteCut)
|
|
&& faceSplitCut_[facei] != absoluteCut
|
|
)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Cut " << faceSplitCut_[facei]
|
|
<< " on face " << mesh().faceCentres()[facei]
|
|
<< " of coupled patch " << pp.name()
|
|
<< " is not consistent with coupled cut "
|
|
<< absoluteCut
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::writeUncutOBJ
|
|
(
|
|
const fileName& dir,
|
|
const label celli
|
|
) const
|
|
{
|
|
// Cell edges
|
|
OFstream cutsStream(dir / "cell_" + name(celli) + ".obj");
|
|
|
|
Pout<< "Writing cell for time " << mesh().time().timeName()
|
|
<< " to " << cutsStream.name() << nl;
|
|
|
|
meshTools::writeOBJ
|
|
(
|
|
cutsStream,
|
|
mesh().cells(),
|
|
mesh().faces(),
|
|
mesh().points(),
|
|
labelList(1, celli)
|
|
);
|
|
|
|
// Loop cutting cell in two
|
|
OFstream cutStream(dir / "cellCuts_" + name(celli) + ".obj");
|
|
|
|
Pout<< "Writing raw cuts on cell for time " << mesh().time().timeName()
|
|
<< " to " << cutStream.name() << nl;
|
|
|
|
const labelList& cPoints = mesh().cellPoints()[celli];
|
|
|
|
forAll(cPoints, i)
|
|
{
|
|
label pointi = cPoints[i];
|
|
if (pointIsCut_[pointi])
|
|
{
|
|
meshTools::writeOBJ(cutStream, mesh().points()[pointi]);
|
|
}
|
|
}
|
|
|
|
const pointField& pts = mesh().points();
|
|
|
|
const labelList& cEdges = mesh().cellEdges()[celli];
|
|
|
|
forAll(cEdges, i)
|
|
{
|
|
label edgeI = cEdges[i];
|
|
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
const edge& e = mesh().edges()[edgeI];
|
|
|
|
const scalar w = edgeWeight_[edgeI];
|
|
|
|
meshTools::writeOBJ(cutStream, w*pts[e[1]] + (1-w)*pts[e[0]]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::writeOBJ
|
|
(
|
|
const fileName& dir,
|
|
const label celli,
|
|
const pointField& loopPoints,
|
|
const labelList& anchors
|
|
) const
|
|
{
|
|
// Cell edges
|
|
OFstream cutsStream(dir / "cell_" + name(celli) + ".obj");
|
|
|
|
Pout<< "Writing cell for time " << mesh().time().timeName()
|
|
<< " to " << cutsStream.name() << nl;
|
|
|
|
meshTools::writeOBJ
|
|
(
|
|
cutsStream,
|
|
mesh().cells(),
|
|
mesh().faces(),
|
|
mesh().points(),
|
|
labelList(1, celli)
|
|
);
|
|
|
|
|
|
// Loop cutting cell in two
|
|
OFstream loopStream(dir / "cellLoop_" + name(celli) + ".obj");
|
|
|
|
Pout<< "Writing loop for time " << mesh().time().timeName()
|
|
<< " to " << loopStream.name() << nl;
|
|
|
|
label vertI = 0;
|
|
|
|
writeOBJ(loopStream, loopPoints, vertI);
|
|
|
|
|
|
// Anchors for cell
|
|
OFstream anchorStream(dir / "anchors_" + name(celli) + ".obj");
|
|
|
|
Pout<< "Writing anchors for time " << mesh().time().timeName()
|
|
<< " to " << anchorStream.name() << endl;
|
|
|
|
forAll(anchors, i)
|
|
{
|
|
meshTools::writeOBJ(anchorStream, mesh().points()[anchors[i]]);
|
|
}
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::edgeEdgeToFace
|
|
(
|
|
const label celli,
|
|
const label edgeA,
|
|
const label edgeB
|
|
) const
|
|
{
|
|
const labelList& cFaces = mesh().cells()[celli];
|
|
|
|
forAll(cFaces, cFacei)
|
|
{
|
|
label facei = cFaces[cFacei];
|
|
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
|
|
if (fEdges.found(edgeA) && fEdges.found(edgeB))
|
|
{
|
|
return facei;
|
|
}
|
|
}
|
|
|
|
// Coming here means the loop is illegal since the two edges
|
|
// are not shared by a face. We just mark loop as invalid and continue.
|
|
|
|
WarningInFunction
|
|
<< "cellCuts : Cannot find face on cell "
|
|
<< celli << " that has both edges " << edgeA << ' ' << edgeB << endl
|
|
<< "faces : " << cFaces << endl
|
|
<< "edgeA : " << mesh().edges()[edgeA] << endl
|
|
<< "edgeB : " << mesh().edges()[edgeB] << endl
|
|
<< "Marking the loop across this cell as invalid" << endl;
|
|
|
|
return -1;
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::edgeVertexToFace
|
|
(
|
|
const label celli,
|
|
const label edgeI,
|
|
const label vertI
|
|
) const
|
|
{
|
|
const labelList& cFaces = mesh().cells()[celli];
|
|
|
|
forAll(cFaces, cFacei)
|
|
{
|
|
label facei = cFaces[cFacei];
|
|
|
|
const face& f = mesh().faces()[facei];
|
|
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
|
|
if (fEdges.found(edgeI) && f.found(vertI))
|
|
{
|
|
return facei;
|
|
}
|
|
}
|
|
|
|
WarningInFunction
|
|
<< "cellCuts : Cannot find face on cell "
|
|
<< celli << " that has both edge " << edgeI << " and vertex "
|
|
<< vertI << endl
|
|
<< "faces : " << cFaces << endl
|
|
<< "edge : " << mesh().edges()[edgeI] << endl
|
|
<< "Marking the loop across this cell as invalid" << endl;
|
|
|
|
return -1;
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::vertexVertexToFace
|
|
(
|
|
const label celli,
|
|
const label vertA,
|
|
const label vertB
|
|
) const
|
|
{
|
|
const labelList& cFaces = mesh().cells()[celli];
|
|
|
|
forAll(cFaces, cFacei)
|
|
{
|
|
label facei = cFaces[cFacei];
|
|
|
|
const face& f = mesh().faces()[facei];
|
|
|
|
if (f.found(vertA) && f.found(vertB))
|
|
{
|
|
return facei;
|
|
}
|
|
}
|
|
|
|
WarningInFunction
|
|
<< "cellCuts : Cannot find face on cell "
|
|
<< celli << " that has vertex " << vertA << " and vertex "
|
|
<< vertB << endl
|
|
<< "faces : " << cFaces << endl
|
|
<< "Marking the loop across this cell as invalid" << endl;
|
|
|
|
return -1;
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::calcFaceCuts() const
|
|
{
|
|
if (faceCutsPtr_.valid())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "faceCuts already calculated" << abort(FatalError);
|
|
}
|
|
|
|
const faceList& faces = mesh().faces();
|
|
|
|
faceCutsPtr_.reset(new labelListList(mesh().nFaces()));
|
|
labelListList& faceCuts = faceCutsPtr_();
|
|
|
|
for (label facei = 0; facei < mesh().nFaces(); facei++)
|
|
{
|
|
const face& f = faces[facei];
|
|
|
|
// Big enough storage (possibly all points and all edges cut). Shrink
|
|
// later on.
|
|
labelList& cuts = faceCuts[facei];
|
|
|
|
cuts.setSize(2*f.size());
|
|
|
|
label cutI = 0;
|
|
|
|
// Do point-edge-point walk over face and collect all cuts.
|
|
// Problem is that we want to start from one of the endpoints of a
|
|
// string of connected cuts; we don't want to start somewhere in the
|
|
// middle.
|
|
|
|
// Pass1: find first point cut not preceded by a cut.
|
|
label startFp = -1;
|
|
|
|
forAll(f, fp)
|
|
{
|
|
label v0 = f[fp];
|
|
|
|
if (pointIsCut_[v0])
|
|
{
|
|
label vMin1 = f[f.rcIndex(fp)];
|
|
|
|
if
|
|
(
|
|
!pointIsCut_[vMin1]
|
|
&& !edgeIsCut_[findEdge(facei, v0, vMin1)]
|
|
)
|
|
{
|
|
cuts[cutI++] = vertToEVert(v0);
|
|
startFp = f.fcIndex(fp);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Pass2: first edge cut not preceded by point cut
|
|
if (startFp == -1)
|
|
{
|
|
forAll(f, fp)
|
|
{
|
|
label fp1 = f.fcIndex(fp);
|
|
|
|
label v0 = f[fp];
|
|
label v1 = f[fp1];
|
|
|
|
label edgeI = findEdge(facei, v0, v1);
|
|
|
|
if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
|
|
{
|
|
cuts[cutI++] = edgeToEVert(edgeI);
|
|
startFp = fp1;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Pass3: nothing found so far. Either face is not cut at all or all
|
|
// vertices are cut. Start from 0.
|
|
if (startFp == -1)
|
|
{
|
|
startFp = 0;
|
|
}
|
|
|
|
// Store all cuts starting from startFp;
|
|
label fp = startFp;
|
|
|
|
bool allVerticesCut = true;
|
|
|
|
forAll(f, i)
|
|
{
|
|
label fp1 = f.fcIndex(fp);
|
|
|
|
// Get the three items: current vertex, next vertex and edge
|
|
// inbetween
|
|
label v0 = f[fp];
|
|
label v1 = f[fp1];
|
|
label edgeI = findEdge(facei, v0, v1);
|
|
|
|
if (pointIsCut_[v0])
|
|
{
|
|
cuts[cutI++] = vertToEVert(v0);
|
|
}
|
|
else
|
|
{
|
|
// For check if all vertices have been cut (= illegal)
|
|
allVerticesCut = false;
|
|
}
|
|
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
cuts[cutI++] = edgeToEVert(edgeI);
|
|
}
|
|
|
|
fp = fp1;
|
|
}
|
|
|
|
|
|
if (allVerticesCut)
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Face " << facei << " vertices " << f
|
|
<< " has all its vertices cut. Not cutting face." << endl;
|
|
}
|
|
|
|
cutI = 0;
|
|
}
|
|
|
|
// Remove duplicate starting point
|
|
if (cutI > 1 && cuts[cutI-1] == cuts[0])
|
|
{
|
|
cutI--;
|
|
}
|
|
cuts.setSize(cutI);
|
|
}
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::findEdge
|
|
(
|
|
const label facei,
|
|
const label v0,
|
|
const label v1
|
|
) const
|
|
{
|
|
const edgeList& edges = mesh().edges();
|
|
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
|
|
forAll(fEdges, i)
|
|
{
|
|
const edge& e = edges[fEdges[i]];
|
|
|
|
if
|
|
(
|
|
(e[0] == v0 && e[1] == v1)
|
|
|| (e[0] == v1 && e[1] == v0)
|
|
)
|
|
{
|
|
return fEdges[i];
|
|
}
|
|
}
|
|
|
|
return -1;
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::loopFace
|
|
(
|
|
const label celli,
|
|
const labelList& loop
|
|
) const
|
|
{
|
|
const cell& cFaces = mesh().cells()[celli];
|
|
|
|
forAll(cFaces, cFacei)
|
|
{
|
|
label facei = cFaces[cFacei];
|
|
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
const face& f = mesh().faces()[facei];
|
|
|
|
bool allOnFace = true;
|
|
|
|
forAll(loop, i)
|
|
{
|
|
label cut = loop[i];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
if (!fEdges.found(getEdge(cut)))
|
|
{
|
|
// Edge not on face. Skip face.
|
|
allOnFace = false;
|
|
break;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (!f.found(getVertex(cut)))
|
|
{
|
|
// Vertex not on face. Skip face.
|
|
allOnFace = false;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (allOnFace)
|
|
{
|
|
// Found face where all elements of loop are on the face.
|
|
return facei;
|
|
}
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::walkPoint
|
|
(
|
|
const label celli,
|
|
const label startCut,
|
|
|
|
const label exclude0,
|
|
const label exclude1,
|
|
|
|
const label otherCut,
|
|
|
|
label& nVisited,
|
|
labelList& visited
|
|
) const
|
|
{
|
|
label vertI = getVertex(otherCut);
|
|
|
|
const labelList& pFaces = mesh().pointFaces()[vertI];
|
|
|
|
forAll(pFaces, pFacei)
|
|
{
|
|
label otherFacei = pFaces[pFacei];
|
|
|
|
if
|
|
(
|
|
otherFacei != exclude0
|
|
&& otherFacei != exclude1
|
|
&& meshTools::faceOnCell(mesh(), celli, otherFacei)
|
|
)
|
|
{
|
|
label oldNVisited = nVisited;
|
|
|
|
bool foundLoop =
|
|
walkCell
|
|
(
|
|
celli,
|
|
startCut,
|
|
otherFacei,
|
|
otherCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
|
|
if (foundLoop)
|
|
{
|
|
return true;
|
|
}
|
|
|
|
// No success. Restore state and continue
|
|
nVisited = oldNVisited;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::crossEdge
|
|
(
|
|
const label celli,
|
|
const label startCut,
|
|
const label facei,
|
|
const label otherCut,
|
|
|
|
label& nVisited,
|
|
labelList& visited
|
|
) const
|
|
{
|
|
// Cross edge to other face
|
|
label edgeI = getEdge(otherCut);
|
|
|
|
label otherFacei = meshTools::otherFace(mesh(), celli, facei, edgeI);
|
|
|
|
// Store old state
|
|
label oldNVisited = nVisited;
|
|
|
|
// Recurse to otherCut
|
|
bool foundLoop =
|
|
walkCell
|
|
(
|
|
celli,
|
|
startCut,
|
|
otherFacei,
|
|
otherCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
|
|
if (foundLoop)
|
|
{
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
// No success. Restore state (i.e. backtrack)
|
|
nVisited = oldNVisited;
|
|
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::addCut
|
|
(
|
|
const label celli,
|
|
const label cut,
|
|
label& nVisited,
|
|
labelList& visited
|
|
) const
|
|
{
|
|
// Check for duplicate cuts.
|
|
if (findPartIndex(visited, nVisited, cut) != -1)
|
|
{
|
|
// Truncate (copy of) visited for ease of printing.
|
|
labelList truncVisited(visited);
|
|
truncVisited.setSize(nVisited);
|
|
|
|
if (verbose_ || debug)
|
|
{
|
|
Pout<< "For cell " << celli << " : trying to add duplicate cut "
|
|
<< cut;
|
|
labelList cuts(1, cut);
|
|
writeCuts(Pout, cuts, loopWeights(cuts));
|
|
|
|
Pout<< " to path:";
|
|
writeCuts(Pout, truncVisited, loopWeights(truncVisited));
|
|
Pout<< endl;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
else
|
|
{
|
|
visited[nVisited++] = cut;
|
|
|
|
return true;
|
|
}
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::walkFace
|
|
(
|
|
const label celli,
|
|
const label startCut,
|
|
const label facei,
|
|
const label cut,
|
|
|
|
label& lastCut,
|
|
label& beforeLastCut,
|
|
label& nVisited,
|
|
labelList& visited
|
|
) const
|
|
{
|
|
const labelList& fCuts = faceCuts()[facei];
|
|
|
|
if (fCuts.size() < 2)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
// Easy case : two cuts.
|
|
if (fCuts.size() == 2)
|
|
{
|
|
if (fCuts[0] == cut)
|
|
{
|
|
if (!addCut(celli, cut, nVisited, visited))
|
|
{
|
|
return false;
|
|
}
|
|
|
|
beforeLastCut = cut;
|
|
lastCut = fCuts[1];
|
|
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
if (!addCut(celli, cut, nVisited, visited))
|
|
{
|
|
return false;
|
|
}
|
|
|
|
beforeLastCut = cut;
|
|
lastCut = fCuts[0];
|
|
|
|
return true;
|
|
}
|
|
}
|
|
|
|
// Harder case: more than 2 cuts on face.
|
|
// Should be start or end of string of cuts. Store all but last cut.
|
|
if (fCuts[0] == cut)
|
|
{
|
|
// Walk forward
|
|
for (label i = 0; i < fCuts.size()-1; i++)
|
|
{
|
|
if (!addCut(celli, fCuts[i], nVisited, visited))
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
beforeLastCut = fCuts[fCuts.size()-2];
|
|
lastCut = fCuts[fCuts.size()-1];
|
|
}
|
|
else if (fCuts[fCuts.size()-1] == cut)
|
|
{
|
|
for (label i = fCuts.size()-1; i >= 1; --i)
|
|
{
|
|
if (!addCut(celli, fCuts[i], nVisited, visited))
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
beforeLastCut = fCuts[1];
|
|
lastCut = fCuts[0];
|
|
}
|
|
else
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "In middle of cut. cell:" << celli << " face:" << facei
|
|
<< " cuts:" << fCuts << " current cut:" << cut << endl;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
|
|
bool Foam::cellCuts::walkCell
|
|
(
|
|
const label celli,
|
|
const label startCut,
|
|
const label facei,
|
|
const label cut,
|
|
|
|
label& nVisited,
|
|
labelList& visited
|
|
) const
|
|
{
|
|
// Walk across face, storing cuts. Return last two cuts visited.
|
|
label lastCut = -1;
|
|
label beforeLastCut = -1;
|
|
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< "For cell:" << celli << " walked across face " << facei
|
|
<< " from cut ";
|
|
labelList cuts(1, cut);
|
|
writeCuts(Pout, cuts, loopWeights(cuts));
|
|
Pout<< endl;
|
|
}
|
|
|
|
bool validWalk = walkFace
|
|
(
|
|
celli,
|
|
startCut,
|
|
facei,
|
|
cut,
|
|
|
|
lastCut,
|
|
beforeLastCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
|
|
if (!validWalk)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< " to last cut ";
|
|
labelList cuts(1, lastCut);
|
|
writeCuts(Pout, cuts, loopWeights(cuts));
|
|
Pout<< endl;
|
|
}
|
|
|
|
// Check if starting point reached.
|
|
if (lastCut == startCut)
|
|
{
|
|
if (nVisited >= 3)
|
|
{
|
|
if (debug & 2)
|
|
{
|
|
// Truncate (copy of) visited for ease of printing.
|
|
labelList truncVisited(visited);
|
|
truncVisited.setSize(nVisited);
|
|
|
|
Pout<< "For cell " << celli << " : found closed path:";
|
|
writeCuts(Pout, truncVisited, loopWeights(truncVisited));
|
|
Pout<< " closed by " << lastCut << endl;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
// Loop but too short. Probably folded back on itself.
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
// Check last cut and one before that to determine type.
|
|
|
|
// From beforeLastCut to lastCut:
|
|
// - from edge to edge
|
|
// (- always not along existing edge)
|
|
// - from edge to vertex
|
|
// - not along existing edge
|
|
// - along existing edge: not allowed?
|
|
// - from vertex to edge
|
|
// - not along existing edge
|
|
// - along existing edge. Not allowed. See above.
|
|
// - from vertex to vertex
|
|
// - not along existing edge
|
|
// - along existing edge
|
|
|
|
if (isEdge(beforeLastCut))
|
|
{
|
|
if (isEdge(lastCut))
|
|
{
|
|
// beforeLastCut=edge, lastCut=edge.
|
|
|
|
// Cross lastCut (=edge) into face not facei
|
|
return crossEdge
|
|
(
|
|
celli,
|
|
startCut,
|
|
facei,
|
|
lastCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
}
|
|
else
|
|
{
|
|
// beforeLastCut=edge, lastCut=vertex.
|
|
|
|
// Go from lastCut to all connected faces (except facei)
|
|
return walkPoint
|
|
(
|
|
celli,
|
|
startCut,
|
|
facei, // exclude0: don't cross back on current face
|
|
-1, // exclude1
|
|
lastCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (isEdge(lastCut))
|
|
{
|
|
// beforeLastCut=vertex, lastCut=edge.
|
|
return crossEdge
|
|
(
|
|
celli,
|
|
startCut,
|
|
facei,
|
|
lastCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
}
|
|
else
|
|
{
|
|
// beforeLastCut=vertex, lastCut=vertex. Check if along existing
|
|
// edge.
|
|
label edgeI =
|
|
findEdge
|
|
(
|
|
facei,
|
|
getVertex(beforeLastCut),
|
|
getVertex(lastCut)
|
|
);
|
|
|
|
if (edgeI != -1)
|
|
{
|
|
// Cut along existing edge. So is in fact on two faces.
|
|
// Get faces on both sides of the edge to make
|
|
// sure we do not fold back on to those.
|
|
|
|
label f0, f1;
|
|
meshTools::getEdgeFaces(mesh(), celli, edgeI, f0, f1);
|
|
|
|
return walkPoint
|
|
(
|
|
celli,
|
|
startCut,
|
|
f0,
|
|
f1,
|
|
lastCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
|
|
}
|
|
else
|
|
{
|
|
// Cross cut across face.
|
|
return walkPoint
|
|
(
|
|
celli,
|
|
startCut,
|
|
facei, // face to exclude
|
|
-1, // face to exclude
|
|
lastCut,
|
|
nVisited,
|
|
visited
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::calcCellLoops(const labelList& cutCells)
|
|
{
|
|
// Determine for every cut cell the loop (= face) it is cut by. Done by
|
|
// starting from a cut edge or cut vertex and walking across faces, from
|
|
// cut to cut, until starting cut hit.
|
|
// If multiple loops are possible across a cell circumference takes the
|
|
// first one found.
|
|
|
|
// Calculate cuts per face.
|
|
const labelListList& allFaceCuts = faceCuts();
|
|
|
|
// Per cell the number of faces with valid cuts. Is used as quick
|
|
// rejection to see if cell can be cut.
|
|
labelList nCutFaces(mesh().nCells(), Zero);
|
|
|
|
forAll(allFaceCuts, facei)
|
|
{
|
|
const labelList& fCuts = allFaceCuts[facei];
|
|
|
|
if (fCuts.size() == mesh().faces()[facei].size())
|
|
{
|
|
// Too many cuts on face. WalkCell would get very upset so disable.
|
|
nCutFaces[mesh().faceOwner()[facei]] = labelMin;
|
|
|
|
if (mesh().isInternalFace(facei))
|
|
{
|
|
nCutFaces[mesh().faceNeighbour()[facei]] = labelMin;
|
|
}
|
|
}
|
|
else if (fCuts.size() >= 2)
|
|
{
|
|
// Could be valid cut. Update count for owner and neighbour.
|
|
nCutFaces[mesh().faceOwner()[facei]]++;
|
|
|
|
if (mesh().isInternalFace(facei))
|
|
{
|
|
nCutFaces[mesh().faceNeighbour()[facei]]++;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Stack of visited cuts (nVisited used as stack pointer)
|
|
// Size big enough.
|
|
labelList visited(mesh().nPoints());
|
|
|
|
forAll(cutCells, i)
|
|
{
|
|
label celli = cutCells[i];
|
|
|
|
bool validLoop = false;
|
|
|
|
// Quick rejection: has enough faces that are cut?
|
|
if (nCutFaces[celli] >= 1)
|
|
{
|
|
const labelList& cFaces = mesh().cells()[celli];
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< "cell:" << celli << " cut faces:" << endl;
|
|
forAll(cFaces, i)
|
|
{
|
|
label facei = cFaces[i];
|
|
const labelList& fCuts = allFaceCuts[facei];
|
|
|
|
Pout<< " face:" << facei << " cuts:";
|
|
writeCuts(Pout, fCuts, loopWeights(fCuts));
|
|
Pout<< endl;
|
|
}
|
|
}
|
|
|
|
label nVisited = 0;
|
|
|
|
// Determine the first cut face to start walking from.
|
|
forAll(cFaces, cFacei)
|
|
{
|
|
label facei = cFaces[cFacei];
|
|
|
|
const labelList& fCuts = allFaceCuts[facei];
|
|
|
|
// Take first or last cut of multiple on face.
|
|
// Note that in calcFaceCuts
|
|
// we have already made sure this is the start or end of a
|
|
// string of cuts.
|
|
if (fCuts.size() >= 2)
|
|
{
|
|
// Try walking from start of fCuts.
|
|
nVisited = 0;
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< "cell:" << celli
|
|
<< " start walk at face:" << facei
|
|
<< " cut:";
|
|
labelList cuts(1, fCuts[0]);
|
|
writeCuts(Pout, cuts, loopWeights(cuts));
|
|
Pout<< endl;
|
|
}
|
|
|
|
validLoop =
|
|
walkCell
|
|
(
|
|
celli,
|
|
fCuts[0],
|
|
facei,
|
|
fCuts[0],
|
|
|
|
nVisited,
|
|
visited
|
|
);
|
|
|
|
if (validLoop)
|
|
{
|
|
break;
|
|
}
|
|
|
|
// No need to try and walk from end of cuts since
|
|
// all paths already tried by walkCell.
|
|
}
|
|
}
|
|
|
|
if (validLoop)
|
|
{
|
|
// Copy nVisited elements out of visited (since visited is
|
|
// never truncated for efficiency reasons)
|
|
|
|
labelList& loop = cellLoops_[celli];
|
|
|
|
loop.setSize(nVisited);
|
|
|
|
for (label i = 0; i < nVisited; i++)
|
|
{
|
|
loop[i] = visited[i];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Invalid loop. Leave cellLoops_[celli] zero size which
|
|
// flags this.
|
|
if (verbose_ || debug)
|
|
{
|
|
Pout<< "calcCellLoops(const labelList&) :"
|
|
<< " did not find valid"
|
|
<< " loop for cell " << celli << endl;
|
|
// Dump cell and cuts on cell.
|
|
writeUncutOBJ(".", celli);
|
|
}
|
|
cellLoops_[celli].setSize(0);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
//Pout<< "calcCellLoops(const labelList&) : did not find valid"
|
|
// << " loop for cell " << celli << " since not enough cut faces"
|
|
// << endl;
|
|
cellLoops_[celli].setSize(0);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::walkEdges
|
|
(
|
|
const label celli,
|
|
const label pointi,
|
|
const label status,
|
|
|
|
Map<label>& edgeStatus,
|
|
Map<label>& pointStatus
|
|
) const
|
|
{
|
|
if (pointStatus.insert(pointi, status))
|
|
{
|
|
// First visit to pointi
|
|
|
|
const labelList& pEdges = mesh().pointEdges()[pointi];
|
|
|
|
forAll(pEdges, pEdgeI)
|
|
{
|
|
label edgeI = pEdges[pEdgeI];
|
|
|
|
if
|
|
(
|
|
meshTools::edgeOnCell(mesh(), celli, edgeI)
|
|
&& edgeStatus.insert(edgeI, status)
|
|
)
|
|
{
|
|
// First visit to edgeI so recurse.
|
|
|
|
label v2 = mesh().edges()[edgeI].otherVertex(pointi);
|
|
|
|
walkEdges(celli, v2, status, edgeStatus, pointStatus);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::labelList Foam::cellCuts::nonAnchorPoints
|
|
(
|
|
const labelList& cellPoints,
|
|
const labelList& anchorPoints,
|
|
const labelList& loop
|
|
) const
|
|
{
|
|
labelList newElems(cellPoints.size());
|
|
label newElemI = 0;
|
|
|
|
forAll(cellPoints, i)
|
|
{
|
|
label pointi = cellPoints[i];
|
|
|
|
if
|
|
(
|
|
!anchorPoints.found(pointi)
|
|
&& !loop.found(vertToEVert(pointi))
|
|
)
|
|
{
|
|
newElems[newElemI++] = pointi;
|
|
}
|
|
}
|
|
|
|
newElems.setSize(newElemI);
|
|
|
|
return newElems;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::loopAnchorConsistent
|
|
(
|
|
const label celli,
|
|
const pointField& loopPts,
|
|
const labelList& anchorPoints
|
|
) const
|
|
{
|
|
// Create identity face for ease of calculation of normal etc.
|
|
face f(identity(loopPts.size()));
|
|
|
|
const vector areaNorm = f.areaNormal(loopPts);
|
|
const point ctr = f.centre(loopPts);
|
|
|
|
|
|
// Get average position of anchor points.
|
|
vector avg(Zero);
|
|
|
|
for (const label pointi : anchorPoints)
|
|
{
|
|
avg += mesh().points()[pointi];
|
|
}
|
|
avg /= anchorPoints.size();
|
|
|
|
|
|
if (((avg - ctr) & areaNorm) > 0)
|
|
{
|
|
return true;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::calcAnchors
|
|
(
|
|
const label celli,
|
|
const labelList& loop,
|
|
const pointField& loopPts,
|
|
|
|
labelList& anchorPoints
|
|
) const
|
|
{
|
|
const edgeList& edges = mesh().edges();
|
|
|
|
const labelList& cPoints = mesh().cellPoints()[celli];
|
|
const labelList& cEdges = mesh().cellEdges()[celli];
|
|
const cell& cFaces = mesh().cells()[celli];
|
|
|
|
// Points on loop
|
|
|
|
// Status of point:
|
|
// 0 - on loop
|
|
// 1 - point set 1
|
|
// 2 - point set 2
|
|
Map<label> pointStatus(2*cPoints.size());
|
|
Map<label> edgeStatus(2*cEdges.size());
|
|
|
|
// Mark loop vertices
|
|
forAll(loop, i)
|
|
{
|
|
label cut = loop[i];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
edgeStatus.insert(getEdge(cut), 0);
|
|
}
|
|
else
|
|
{
|
|
pointStatus.insert(getVertex(cut), 0);
|
|
}
|
|
}
|
|
// Since edges between two cut vertices have not been marked, mark them
|
|
// explicitly
|
|
forAll(cEdges, i)
|
|
{
|
|
label edgeI = cEdges[i];
|
|
const edge& e = edges[edgeI];
|
|
|
|
if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
|
|
{
|
|
edgeStatus.insert(edgeI, 0);
|
|
}
|
|
}
|
|
|
|
|
|
// Find uncut starting vertex
|
|
label uncutIndex = firstUnique(cPoints, pointStatus);
|
|
|
|
if (uncutIndex == -1)
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell " << celli << endl
|
|
<< "Can not find point on cell which is not cut by loop."
|
|
<< endl;
|
|
|
|
writeOBJ(".", celli, loopPts, labelList(0));
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
// Walk unset vertices and edges and mark with 1 in pointStatus, edgeStatus
|
|
walkEdges(celli, cPoints[uncutIndex], 1, edgeStatus, pointStatus);
|
|
|
|
// Find new uncut starting vertex
|
|
uncutIndex = firstUnique(cPoints, pointStatus);
|
|
|
|
if (uncutIndex == -1)
|
|
{
|
|
// All vertices either in loop or in anchor. So split is along single
|
|
// face.
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell " << celli << endl
|
|
<< "All vertices of cell are either in loop or in anchor set"
|
|
<< endl;
|
|
|
|
writeOBJ(".", celli, loopPts, labelList(0));
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
// Walk unset vertices and edges and mark with 2. These are the
|
|
// pointset 2.
|
|
walkEdges(celli, cPoints[uncutIndex], 2, edgeStatus, pointStatus);
|
|
|
|
// Collect both sets in lists.
|
|
DynamicList<label> connectedPoints(cPoints.size());
|
|
DynamicList<label> otherPoints(cPoints.size());
|
|
|
|
forAllConstIters(pointStatus, iter)
|
|
{
|
|
if (iter.val() == 1)
|
|
{
|
|
connectedPoints.append(iter.key());
|
|
}
|
|
else if (iter.val() == 2)
|
|
{
|
|
otherPoints.append(iter.key());
|
|
}
|
|
}
|
|
connectedPoints.shrink();
|
|
otherPoints.shrink();
|
|
|
|
// Check that all points have been used.
|
|
uncutIndex = firstUnique(cPoints, pointStatus);
|
|
|
|
if (uncutIndex != -1)
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell " << celli
|
|
<< " since it splits the cell into more than two cells" << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, connectedPoints);
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
// Check that both parts (connectedPoints, otherPoints) have enough faces.
|
|
labelHashSet connectedFaces(2*cFaces.size());
|
|
labelHashSet otherFaces(2*cFaces.size());
|
|
|
|
forAllConstIters(pointStatus, iter)
|
|
{
|
|
const label pointi = iter.key();
|
|
|
|
const labelList& pFaces = mesh().pointFaces()[pointi];
|
|
|
|
if (iter.val() == 1)
|
|
{
|
|
forAll(pFaces, pFacei)
|
|
{
|
|
if (meshTools::faceOnCell(mesh(), celli, pFaces[pFacei]))
|
|
{
|
|
connectedFaces.insert(pFaces[pFacei]);
|
|
}
|
|
}
|
|
}
|
|
else if (iter.val() == 2)
|
|
{
|
|
forAll(pFaces, pFacei)
|
|
{
|
|
if (meshTools::faceOnCell(mesh(), celli, pFaces[pFacei]))
|
|
{
|
|
otherFaces.insert(pFaces[pFacei]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (connectedFaces.size() < 3)
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell " << celli
|
|
<< " since would have too few faces on one side." << nl
|
|
<< "All faces:" << cFaces << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, connectedPoints);
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
if (otherFaces.size() < 3)
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell " << celli
|
|
<< " since would have too few faces on one side." << nl
|
|
<< "All faces:" << cFaces << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, otherPoints);
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
// Check that faces are split into two regions and not more.
|
|
// When walking across the face the only transition of pointStatus is
|
|
// from set1 to loop to set2 or back. Not allowed is from set1 to loop to
|
|
// set1.
|
|
{
|
|
forAll(cFaces, i)
|
|
{
|
|
label facei = cFaces[i];
|
|
|
|
const face& f = mesh().faces()[facei];
|
|
|
|
bool hasSet1 = false;
|
|
bool hasSet2 = false;
|
|
|
|
label prevStat = pointStatus[f[0]];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
label v0 = f[fp];
|
|
label pStat = pointStatus[v0];
|
|
|
|
if (pStat == prevStat)
|
|
{
|
|
}
|
|
else if (pStat == 0)
|
|
{
|
|
// Loop.
|
|
}
|
|
else if (pStat == 1)
|
|
{
|
|
if (hasSet1)
|
|
{
|
|
// Second occurrence of set1.
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell "
|
|
<< celli
|
|
<< " since face " << f << " would be split into"
|
|
<< " more than two faces" << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, otherPoints);
|
|
}
|
|
return false;
|
|
}
|
|
|
|
hasSet1 = true;
|
|
}
|
|
else if (pStat == 2)
|
|
{
|
|
if (hasSet2)
|
|
{
|
|
// Second occurrence of set1.
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell "
|
|
<< celli
|
|
<< " since face " << f << " would be split into"
|
|
<< " more than two faces" << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, otherPoints);
|
|
}
|
|
return false;
|
|
}
|
|
|
|
hasSet2 = true;
|
|
}
|
|
else
|
|
{
|
|
FatalErrorInFunction
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
prevStat = pStat;
|
|
|
|
|
|
label v1 = f.nextLabel(fp);
|
|
label edgeI = findEdge(facei, v0, v1);
|
|
|
|
label eStat = edgeStatus[edgeI];
|
|
|
|
if (eStat == prevStat)
|
|
{
|
|
}
|
|
else if (eStat == 0)
|
|
{
|
|
// Loop.
|
|
}
|
|
else if (eStat == 1)
|
|
{
|
|
if (hasSet1)
|
|
{
|
|
// Second occurrence of set1.
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell "
|
|
<< celli
|
|
<< " since face " << f << " would be split into"
|
|
<< " more than two faces" << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, otherPoints);
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
hasSet1 = true;
|
|
}
|
|
else if (eStat == 2)
|
|
{
|
|
if (hasSet2)
|
|
{
|
|
// Second occurrence of set1.
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Invalid loop " << loop << " for cell "
|
|
<< celli
|
|
<< " since face " << f << " would be split into"
|
|
<< " more than two faces" << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, otherPoints);
|
|
}
|
|
return false;
|
|
}
|
|
|
|
hasSet2 = true;
|
|
}
|
|
prevStat = eStat;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
// Check which one of point sets to use.
|
|
bool loopOk = loopAnchorConsistent(celli, loopPts, connectedPoints);
|
|
|
|
//if (debug)
|
|
{
|
|
// Additional check: are non-anchor points on other side?
|
|
bool otherLoopOk = loopAnchorConsistent(celli, loopPts, otherPoints);
|
|
|
|
if (loopOk == otherLoopOk)
|
|
{
|
|
// Both sets of points are supposedly on the same side as the
|
|
// loop normal. Oops.
|
|
|
|
WarningInFunction
|
|
<< " For cell:" << celli
|
|
<< " achorpoints and nonanchorpoints are geometrically"
|
|
<< " on same side!" << endl
|
|
<< "cellPoints:" << cPoints << endl
|
|
<< "loop:" << loop << endl
|
|
<< "anchors:" << connectedPoints << endl
|
|
<< "otherPoints:" << otherPoints << endl;
|
|
|
|
writeOBJ(".", celli, loopPts, connectedPoints);
|
|
}
|
|
}
|
|
|
|
if (loopOk)
|
|
{
|
|
// connectedPoints on 'outside' of loop so these are anchor points
|
|
anchorPoints.transfer(connectedPoints);
|
|
connectedPoints.clear();
|
|
}
|
|
else
|
|
{
|
|
anchorPoints.transfer(otherPoints);
|
|
otherPoints.clear();
|
|
}
|
|
return true;
|
|
}
|
|
|
|
|
|
Foam::pointField Foam::cellCuts::loopPoints
|
|
(
|
|
const labelList& loop,
|
|
const scalarField& loopWeights
|
|
) const
|
|
{
|
|
pointField loopPts(loop.size());
|
|
|
|
forAll(loop, fp)
|
|
{
|
|
loopPts[fp] = coord(loop[fp], loopWeights[fp]);
|
|
}
|
|
return loopPts;
|
|
}
|
|
|
|
|
|
Foam::scalarField Foam::cellCuts::loopWeights(const labelList& loop) const
|
|
{
|
|
scalarField weights(loop.size());
|
|
|
|
forAll(loop, fp)
|
|
{
|
|
label cut = loop[fp];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
label edgeI = getEdge(cut);
|
|
|
|
weights[fp] = edgeWeight_[edgeI];
|
|
}
|
|
else
|
|
{
|
|
weights[fp] = -GREAT;
|
|
}
|
|
}
|
|
return weights;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::validEdgeLoop
|
|
(
|
|
const labelList& loop,
|
|
const scalarField& loopWeights
|
|
) const
|
|
{
|
|
forAll(loop, fp)
|
|
{
|
|
label cut = loop[fp];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
label edgeI = getEdge(cut);
|
|
|
|
// Check: cut compatible only if can be snapped to existing one.
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
scalar edgeLen = mesh().edges()[edgeI].mag(mesh().points());
|
|
|
|
if
|
|
(
|
|
mag(loopWeights[fp] - edgeWeight_[edgeI])
|
|
> geomCellLooper::snapTol()*edgeLen
|
|
)
|
|
{
|
|
// Positions differ too much->would create two cuts.
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
|
|
Foam::label Foam::cellCuts::countFaceCuts
|
|
(
|
|
const label facei,
|
|
const labelList& loop
|
|
) const
|
|
{
|
|
// Includes cuts through vertices and through edges.
|
|
// Assumes that if edge is cut both in edgeIsCut and in loop that the
|
|
// position of the cut is the same.
|
|
|
|
label nCuts = 0;
|
|
|
|
// Count cut vertices
|
|
const face& f = mesh().faces()[facei];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
label vertI = f[fp];
|
|
|
|
// Vertex already cut or mentioned in current loop.
|
|
if
|
|
(
|
|
pointIsCut_[vertI]
|
|
|| loop.found(vertToEVert(vertI))
|
|
)
|
|
{
|
|
nCuts++;
|
|
}
|
|
}
|
|
|
|
// Count cut edges.
|
|
const labelList& fEdges = mesh().faceEdges()[facei];
|
|
|
|
forAll(fEdges, fEdgeI)
|
|
{
|
|
label edgeI = fEdges[fEdgeI];
|
|
|
|
// Edge already cut or mentioned in current loop.
|
|
if
|
|
(
|
|
edgeIsCut_[edgeI]
|
|
|| loop.found(edgeToEVert(edgeI))
|
|
)
|
|
{
|
|
nCuts++;
|
|
}
|
|
}
|
|
|
|
return nCuts;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::conservativeValidLoop
|
|
(
|
|
const label celli,
|
|
const labelList& loop
|
|
) const
|
|
{
|
|
|
|
if (loop.size() < 2)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
forAll(loop, cutI)
|
|
{
|
|
if (isEdge(loop[cutI]))
|
|
{
|
|
label edgeI = getEdge(loop[cutI]);
|
|
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
// edge compatibility already checked.
|
|
}
|
|
else
|
|
{
|
|
// Quick rejection: vertices of edge itself cannot be cut.
|
|
const edge& e = mesh().edges()[edgeI];
|
|
|
|
if (pointIsCut_[e.start()] || pointIsCut_[e.end()])
|
|
{
|
|
return false;
|
|
}
|
|
|
|
|
|
// Check faces using this edge
|
|
const labelList& eFaces = mesh().edgeFaces()[edgeI];
|
|
|
|
forAll(eFaces, eFacei)
|
|
{
|
|
label nCuts = countFaceCuts(eFaces[eFacei], loop);
|
|
|
|
if (nCuts > 2)
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Vertex cut
|
|
|
|
label vertI = getVertex(loop[cutI]);
|
|
|
|
if (!pointIsCut_[vertI])
|
|
{
|
|
// New cut through vertex.
|
|
|
|
// Check edges using vertex.
|
|
const labelList& pEdges = mesh().pointEdges()[vertI];
|
|
|
|
forAll(pEdges, pEdgeI)
|
|
{
|
|
label edgeI = pEdges[pEdgeI];
|
|
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
// Check faces using vertex.
|
|
const labelList& pFaces = mesh().pointFaces()[vertI];
|
|
|
|
forAll(pFaces, pFacei)
|
|
{
|
|
label nCuts = countFaceCuts(pFaces[pFacei], loop);
|
|
|
|
if (nCuts > 2)
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// All ok.
|
|
return true;
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::validLoop
|
|
(
|
|
const label celli,
|
|
const labelList& loop,
|
|
const scalarField& loopWeights,
|
|
|
|
Map<edge>& newFaceSplitCut,
|
|
labelList& anchorPoints
|
|
) const
|
|
{
|
|
// Determine compatibility of loop with existing cut pattern. Does not use
|
|
// derived cut-addressing (faceCuts), only pointIsCut, edgeIsCut.
|
|
// Adds any cross-cuts found to newFaceSplitCut and sets cell points on
|
|
// one side of the loop in anchorPoints.
|
|
|
|
if (loop.size() < 2)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
if (debug & 4)
|
|
{
|
|
// Allow as fallback the 'old' loop checking where only a single
|
|
// cut per face is allowed.
|
|
if (!conservativeValidLoop(celli, loop))
|
|
{
|
|
Info << "Invalid conservative loop: " << loop << endl;
|
|
return false;
|
|
}
|
|
}
|
|
|
|
forAll(loop, fp)
|
|
{
|
|
label cut = loop[fp];
|
|
label nextCut = loop[(fp+1) % loop.size()];
|
|
|
|
// Label (if any) of face cut (so cut not along existing edge)
|
|
label meshFacei = -1;
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
label edgeI = getEdge(cut);
|
|
|
|
// Look one cut ahead to find if it is along existing edge.
|
|
|
|
if (isEdge(nextCut))
|
|
{
|
|
// From edge to edge -> cross cut
|
|
label nextEdgeI = getEdge(nextCut);
|
|
|
|
// Find face and mark as to be split.
|
|
meshFacei = edgeEdgeToFace(celli, edgeI, nextEdgeI);
|
|
|
|
if (meshFacei == -1)
|
|
{
|
|
// Can't find face using both cut edges.
|
|
return false;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// From edge to vertex -> cross cut only if no existing edge.
|
|
|
|
label nextVertI = getVertex(nextCut);
|
|
const edge& e = mesh().edges()[edgeI];
|
|
|
|
if (e.start() != nextVertI && e.end() != nextVertI)
|
|
{
|
|
// New edge. Find face and mark as to be split.
|
|
meshFacei = edgeVertexToFace(celli, edgeI, nextVertI);
|
|
|
|
if (meshFacei == -1)
|
|
{
|
|
// Can't find face. Illegal.
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Cut is vertex.
|
|
label vertI = getVertex(cut);
|
|
|
|
if (isEdge(nextCut))
|
|
{
|
|
// From vertex to edge -> cross cut only if no existing edge.
|
|
label nextEdgeI = getEdge(nextCut);
|
|
|
|
const edge& nextE = mesh().edges()[nextEdgeI];
|
|
|
|
if (nextE.start() != vertI && nextE.end() != vertI)
|
|
{
|
|
// Should be cross cut. Find face.
|
|
meshFacei = edgeVertexToFace(celli, nextEdgeI, vertI);
|
|
|
|
if (meshFacei == -1)
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// From vertex to vertex -> cross cut only if no existing edge.
|
|
label nextVertI = getVertex(nextCut);
|
|
|
|
if (meshTools::findEdge(mesh(), vertI, nextVertI) == -1)
|
|
{
|
|
// New edge. Find face.
|
|
meshFacei = vertexVertexToFace(celli, vertI, nextVertI);
|
|
|
|
if (meshFacei == -1)
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (meshFacei != -1)
|
|
{
|
|
// meshFacei is cut across along cut-nextCut (so not along existing
|
|
// edge). Check if this is compatible with existing pattern.
|
|
edge cutEdge(cut, nextCut);
|
|
|
|
const auto iter = faceSplitCut_.cfind(meshFacei);
|
|
|
|
if (!iter.found())
|
|
{
|
|
// Face not yet cut so insert.
|
|
newFaceSplitCut.insert(meshFacei, cutEdge);
|
|
}
|
|
else
|
|
{
|
|
// Face already cut. Ok if same edge.
|
|
if (iter.val() != cutEdge)
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Is there a face on which all cuts are?
|
|
label faceContainingLoop = loopFace(celli, loop);
|
|
|
|
if (faceContainingLoop != -1)
|
|
{
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Found loop on cell " << celli << " with all points"
|
|
<< " on face " << faceContainingLoop << endl;
|
|
}
|
|
|
|
//writeOBJ(".", celli, loopPoints(loop, loopWeights), labelList(0));
|
|
|
|
return false;
|
|
}
|
|
|
|
// Calculate anchor points
|
|
// Final success is determined by whether anchor points can be determined.
|
|
return calcAnchors
|
|
(
|
|
celli,
|
|
loop,
|
|
loopPoints(loop, loopWeights),
|
|
anchorPoints
|
|
);
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::setFromCellLoops()
|
|
{
|
|
// 'Uncut' edges/vertices that are not used in loops.
|
|
pointIsCut_ = false;
|
|
|
|
edgeIsCut_ = false;
|
|
|
|
faceSplitCut_.clear();
|
|
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
const labelList& loop = cellLoops_[celli];
|
|
|
|
if (loop.size())
|
|
{
|
|
// Storage for cross-face cuts
|
|
Map<edge> faceSplitCuts(loop.size());
|
|
// Storage for points on one side of cell.
|
|
labelList anchorPoints;
|
|
|
|
if
|
|
(
|
|
!validLoop
|
|
(
|
|
celli,
|
|
loop,
|
|
loopWeights(loop),
|
|
faceSplitCuts,
|
|
anchorPoints
|
|
)
|
|
)
|
|
{
|
|
//writeOBJ(".", celli, loopPoints(celli), anchorPoints);
|
|
if (verbose_ || debug)
|
|
{
|
|
WarningInFunction
|
|
<< "Illegal loop " << loop
|
|
<< " when recreating cut-addressing"
|
|
<< " from existing cellLoops for cell " << celli
|
|
<< endl;
|
|
}
|
|
|
|
cellLoops_[celli].setSize(0);
|
|
cellAnchorPoints_[celli].setSize(0);
|
|
}
|
|
else
|
|
{
|
|
// Copy anchor points.
|
|
cellAnchorPoints_[celli].transfer(anchorPoints);
|
|
|
|
// Copy faceSplitCuts into overall faceSplit info.
|
|
forAllConstIters(faceSplitCuts, iter)
|
|
{
|
|
faceSplitCut_.insert(iter.key(), iter.val());
|
|
}
|
|
|
|
// Update edgeIsCut, pointIsCut information
|
|
forAll(loop, cutI)
|
|
{
|
|
label cut = loop[cutI];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
edgeIsCut_[getEdge(cut)] = true;
|
|
}
|
|
else
|
|
{
|
|
pointIsCut_[getVertex(cut)] = true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Reset edge weights
|
|
forAll(edgeIsCut_, edgeI)
|
|
{
|
|
if (!edgeIsCut_[edgeI])
|
|
{
|
|
// Weight not used. Set to illegal value to make any use fall over.
|
|
edgeWeight_[edgeI] = -GREAT;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
bool Foam::cellCuts::setFromCellLoop
|
|
(
|
|
const label celli,
|
|
const labelList& loop,
|
|
const scalarField& loopWeights
|
|
)
|
|
{
|
|
// Update basic cut information from single cellLoop. Returns true if loop
|
|
// was valid.
|
|
|
|
// Dump loop for debugging.
|
|
if (debug)
|
|
{
|
|
OFstream str("last_cell.obj");
|
|
|
|
str<< "# edges of cell " << celli << nl;
|
|
|
|
meshTools::writeOBJ
|
|
(
|
|
str,
|
|
mesh().cells(),
|
|
mesh().faces(),
|
|
mesh().points(),
|
|
labelList(1, celli)
|
|
);
|
|
|
|
|
|
OFstream loopStr("last_loop.obj");
|
|
|
|
loopStr<< "# looppoints for cell " << celli << nl;
|
|
|
|
pointField pointsOfLoop = loopPoints(loop, loopWeights);
|
|
|
|
forAll(pointsOfLoop, i)
|
|
{
|
|
meshTools::writeOBJ(loopStr, pointsOfLoop[i]);
|
|
}
|
|
|
|
str << 'l';
|
|
|
|
forAll(pointsOfLoop, i)
|
|
{
|
|
str << ' ' << i + 1;
|
|
}
|
|
str << ' ' << 1 << nl;
|
|
}
|
|
|
|
bool okLoop = false;
|
|
|
|
if (validEdgeLoop(loop, loopWeights))
|
|
{
|
|
// Storage for cuts across face
|
|
Map<edge> faceSplitCuts(loop.size());
|
|
// Storage for points on one side of cell
|
|
labelList anchorPoints;
|
|
|
|
okLoop =
|
|
validLoop(celli, loop, loopWeights, faceSplitCuts, anchorPoints);
|
|
|
|
if (okLoop)
|
|
{
|
|
// Valid loop. Copy cellLoops and anchorPoints
|
|
cellLoops_[celli] = loop;
|
|
cellAnchorPoints_[celli].transfer(anchorPoints);
|
|
|
|
// Copy split cuts
|
|
forAllConstIters(faceSplitCuts, iter)
|
|
{
|
|
faceSplitCut_.insert(iter.key(), iter.val());
|
|
}
|
|
|
|
|
|
// Update edgeIsCut, pointIsCut information
|
|
forAll(loop, cutI)
|
|
{
|
|
const label cut = loop[cutI];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
const label edgeI = getEdge(cut);
|
|
|
|
edgeIsCut_[edgeI] = true;
|
|
|
|
edgeWeight_[edgeI] = loopWeights[cutI];
|
|
}
|
|
else
|
|
{
|
|
const label vertI = getVertex(cut);
|
|
|
|
pointIsCut_[vertI] = true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return okLoop;
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::setFromCellLoops
|
|
(
|
|
const labelList& cellLabels,
|
|
const labelListList& cellLoops,
|
|
const List<scalarField>& cellLoopWeights
|
|
)
|
|
{
|
|
// 'Uncut' edges/vertices that are not used in loops.
|
|
pointIsCut_ = false;
|
|
|
|
edgeIsCut_ = false;
|
|
|
|
forAll(cellLabels, cellLabelI)
|
|
{
|
|
const label celli = cellLabels[cellLabelI];
|
|
|
|
const labelList& loop = cellLoops[cellLabelI];
|
|
|
|
if (loop.size())
|
|
{
|
|
const scalarField& loopWeights = cellLoopWeights[cellLabelI];
|
|
|
|
if (setFromCellLoop(celli, loop, loopWeights))
|
|
{
|
|
// Valid loop. Call above will have updated all already.
|
|
}
|
|
else
|
|
{
|
|
// Clear cellLoops
|
|
cellLoops_[celli].setSize(0);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::setFromCellCutter
|
|
(
|
|
const cellLooper& cellCutter,
|
|
const List<refineCell>& refCells
|
|
)
|
|
{
|
|
// 'Uncut' edges/vertices that are not used in loops.
|
|
pointIsCut_ = false;
|
|
|
|
edgeIsCut_ = false;
|
|
|
|
// storage for loop of cuts (cut vertices and/or cut edges)
|
|
labelList cellLoop;
|
|
scalarField cellLoopWeights;
|
|
|
|
// For debugging purposes
|
|
DynamicList<label> invalidCutCells(2);
|
|
DynamicList<labelList> invalidCutLoops(2);
|
|
DynamicList<scalarField> invalidCutLoopWeights(2);
|
|
|
|
|
|
forAll(refCells, refCelli)
|
|
{
|
|
const refineCell& refCell = refCells[refCelli];
|
|
|
|
const label celli = refCell.cellNo();
|
|
|
|
const vector& refDir = refCell.direction();
|
|
|
|
|
|
// Cut cell. Determines cellLoop and cellLoopWeights
|
|
bool goodCut =
|
|
cellCutter.cut
|
|
(
|
|
refDir,
|
|
celli,
|
|
|
|
pointIsCut_,
|
|
edgeIsCut_,
|
|
edgeWeight_,
|
|
|
|
cellLoop,
|
|
cellLoopWeights
|
|
);
|
|
|
|
// Check whether edge refinement is on a per face basis compatible with
|
|
// current pattern.
|
|
if (goodCut)
|
|
{
|
|
if (setFromCellLoop(celli, cellLoop, cellLoopWeights))
|
|
{
|
|
// Valid loop. Will have updated all info already.
|
|
}
|
|
else
|
|
{
|
|
cellLoops_[celli].setSize(0);
|
|
|
|
WarningInFunction
|
|
<< "Found loop on cell " << celli
|
|
<< " that resulted in an unexpected bad cut." << nl
|
|
<< " Suggestions:" << nl
|
|
<< " - Turn on the debug switch for 'cellCuts' to get"
|
|
<< " geometry files that identify this cell." << nl
|
|
<< " - Also keep in mind to check the defined"
|
|
<< " reference directions, as these are most likely the"
|
|
<< " origin of the problem."
|
|
<< nl << endl;
|
|
|
|
// Discarded by validLoop
|
|
if (debug)
|
|
{
|
|
invalidCutCells.append(celli);
|
|
invalidCutLoops.append(cellLoop);
|
|
invalidCutLoopWeights.append(cellLoopWeights);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Clear cellLoops
|
|
cellLoops_[celli].setSize(0);
|
|
}
|
|
}
|
|
|
|
if (debug && invalidCutCells.size())
|
|
{
|
|
invalidCutCells.shrink();
|
|
invalidCutLoops.shrink();
|
|
invalidCutLoopWeights.shrink();
|
|
|
|
fileName cutsFile("invalidLoopCells.obj");
|
|
|
|
Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;
|
|
|
|
OFstream cutsStream(cutsFile);
|
|
|
|
meshTools::writeOBJ
|
|
(
|
|
cutsStream,
|
|
mesh().cells(),
|
|
mesh().faces(),
|
|
mesh().points(),
|
|
invalidCutCells
|
|
);
|
|
|
|
fileName loopsFile("invalidLoops.obj");
|
|
|
|
Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;
|
|
|
|
OFstream loopsStream(loopsFile);
|
|
|
|
label vertI = 0;
|
|
|
|
forAll(invalidCutLoops, i)
|
|
{
|
|
writeOBJ
|
|
(
|
|
loopsStream,
|
|
loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
|
|
vertI
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::setFromCellCutter
|
|
(
|
|
const cellLooper& cellCutter,
|
|
const labelList& cellLabels,
|
|
const List<plane>& cellCutPlanes
|
|
)
|
|
{
|
|
// 'Uncut' edges/vertices that are not used in loops.
|
|
pointIsCut_ = false;
|
|
|
|
edgeIsCut_ = false;
|
|
|
|
// storage for loop of cuts (cut vertices and/or cut edges)
|
|
labelList cellLoop;
|
|
scalarField cellLoopWeights;
|
|
|
|
// For debugging purposes
|
|
DynamicList<label> invalidCutCells(2);
|
|
DynamicList<labelList> invalidCutLoops(2);
|
|
DynamicList<scalarField> invalidCutLoopWeights(2);
|
|
|
|
|
|
forAll(cellLabels, i)
|
|
{
|
|
const label celli = cellLabels[i];
|
|
|
|
// Cut cell. Determines cellLoop and cellLoopWeights
|
|
bool goodCut =
|
|
cellCutter.cut
|
|
(
|
|
cellCutPlanes[i],
|
|
celli,
|
|
|
|
pointIsCut_,
|
|
edgeIsCut_,
|
|
edgeWeight_,
|
|
|
|
cellLoop,
|
|
cellLoopWeights
|
|
);
|
|
|
|
// Check whether edge refinement is on a per face basis compatible with
|
|
// current pattern.
|
|
if (goodCut)
|
|
{
|
|
if (setFromCellLoop(celli, cellLoop, cellLoopWeights))
|
|
{
|
|
// Valid loop. Will have updated all info already.
|
|
}
|
|
else
|
|
{
|
|
cellLoops_[celli].setSize(0);
|
|
|
|
// Discarded by validLoop
|
|
if (debug)
|
|
{
|
|
invalidCutCells.append(celli);
|
|
invalidCutLoops.append(cellLoop);
|
|
invalidCutLoopWeights.append(cellLoopWeights);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Clear cellLoops
|
|
cellLoops_[celli].setSize(0);
|
|
}
|
|
}
|
|
|
|
if (debug && invalidCutCells.size())
|
|
{
|
|
invalidCutCells.shrink();
|
|
invalidCutLoops.shrink();
|
|
invalidCutLoopWeights.shrink();
|
|
|
|
fileName cutsFile("invalidLoopCells.obj");
|
|
|
|
Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;
|
|
|
|
OFstream cutsStream(cutsFile);
|
|
|
|
meshTools::writeOBJ
|
|
(
|
|
cutsStream,
|
|
mesh().cells(),
|
|
mesh().faces(),
|
|
mesh().points(),
|
|
invalidCutCells
|
|
);
|
|
|
|
fileName loopsFile("invalidLoops.obj");
|
|
|
|
Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;
|
|
|
|
OFstream loopsStream(loopsFile);
|
|
|
|
label vertI = 0;
|
|
|
|
forAll(invalidCutLoops, i)
|
|
{
|
|
writeOBJ
|
|
(
|
|
loopsStream,
|
|
loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
|
|
vertI
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::orientPlanesAndLoops()
|
|
{
|
|
// Determine anchorPoints if not yet done by validLoop.
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
const labelList& loop = cellLoops_[celli];
|
|
|
|
if (loop.size() && cellAnchorPoints_[celli].empty())
|
|
{
|
|
// Leave anchor points empty if illegal loop.
|
|
calcAnchors
|
|
(
|
|
celli,
|
|
loop,
|
|
loopPoints(celli),
|
|
cellAnchorPoints_[celli]
|
|
);
|
|
}
|
|
}
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< "cellAnchorPoints:" << endl;
|
|
}
|
|
forAll(cellAnchorPoints_, celli)
|
|
{
|
|
if (cellLoops_[celli].size())
|
|
{
|
|
if (cellAnchorPoints_[celli].empty())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "No anchor points for cut cell " << celli << endl
|
|
<< "cellLoop:" << cellLoops_[celli] << abort(FatalError);
|
|
}
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< " cell:" << celli << " anchored at "
|
|
<< cellAnchorPoints_[celli] << endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Calculate number of valid cellLoops
|
|
nLoops_ = 0;
|
|
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
if (cellLoops_[celli].size())
|
|
{
|
|
nLoops_++;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::calcLoopsAndAddressing(const labelList& cutCells)
|
|
{
|
|
// Sanity check on weights
|
|
forAll(edgeIsCut_, edgeI)
|
|
{
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
scalar weight = edgeWeight_[edgeI];
|
|
|
|
if (weight < 0 || weight > 1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Weight out of range [0,1]. Edge " << edgeI
|
|
<< " verts:" << mesh().edges()[edgeI]
|
|
<< " weight:" << weight << abort(FatalError);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Weight not used. Set to illegal value to make any use fall over.
|
|
edgeWeight_[edgeI] = -GREAT;
|
|
}
|
|
}
|
|
|
|
|
|
// Calculate faces that split cells in two
|
|
calcCellLoops(cutCells);
|
|
|
|
if (debug & 2)
|
|
{
|
|
Pout<< "-- cellLoops --" << endl;
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
const labelList& loop = cellLoops_[celli];
|
|
|
|
if (loop.size())
|
|
{
|
|
Pout<< "cell:" << celli << " ";
|
|
writeCuts(Pout, loop, loopWeights(loop));
|
|
Pout<< endl;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Redo basic cut information (pointIsCut, edgeIsCut, faceSplitCut)
|
|
// using cellLoop only.
|
|
setFromCellLoops();
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::check() const
|
|
{
|
|
// Check weights for unsnapped values
|
|
forAll(edgeIsCut_, edgeI)
|
|
{
|
|
if (edgeIsCut_[edgeI])
|
|
{
|
|
if
|
|
(
|
|
edgeWeight_[edgeI] < geomCellLooper::snapTol()
|
|
|| edgeWeight_[edgeI] > (1 - geomCellLooper::snapTol())
|
|
)
|
|
{
|
|
// Should have been snapped.
|
|
WarningInFunction
|
|
<< "edge:" << edgeI << " vertices:"
|
|
<< mesh().edges()[edgeI]
|
|
<< " weight:" << edgeWeight_[edgeI] << " should have been"
|
|
<< " snapped to one of its endpoints"
|
|
<< endl;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (edgeWeight_[edgeI] > - 1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "edge:" << edgeI << " vertices:"
|
|
<< mesh().edges()[edgeI]
|
|
<< " weight:" << edgeWeight_[edgeI] << " is not cut but"
|
|
<< " its weight is not marked invalid"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Check that all elements of cellloop are registered
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
const labelList& loop = cellLoops_[celli];
|
|
|
|
forAll(loop, i)
|
|
{
|
|
label cut = loop[i];
|
|
|
|
if
|
|
(
|
|
(isEdge(cut) && !edgeIsCut_[getEdge(cut)])
|
|
|| (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
|
|
)
|
|
{
|
|
labelList cuts(1, cut);
|
|
writeCuts(Pout, cuts, loopWeights(cuts));
|
|
|
|
FatalErrorInFunction
|
|
<< "cell:" << celli << " loop:"
|
|
<< loop
|
|
<< " cut:" << cut << " is not marked as cut"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Check that no elements of cell loop are anchor point.
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
const labelList& anchors = cellAnchorPoints_[celli];
|
|
|
|
const labelList& loop = cellLoops_[celli];
|
|
|
|
if (loop.size() && anchors.empty())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "cell:" << celli << " loop:" << loop
|
|
<< " has no anchor points"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
|
|
forAll(loop, i)
|
|
{
|
|
label cut = loop[i];
|
|
|
|
if
|
|
(
|
|
!isEdge(cut)
|
|
&& anchors.found(getVertex(cut))
|
|
)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "cell:" << celli << " loop:" << loop
|
|
<< " anchor points:" << anchors
|
|
<< " anchor:" << getVertex(cut) << " is part of loop"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Check that cut faces have a neighbour that is cut.
|
|
boolList nbrCellIsCut;
|
|
{
|
|
boolList cellIsCut(mesh().nCells(), false);
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
cellIsCut[celli] = cellLoops_[celli].size();
|
|
}
|
|
syncTools::swapBoundaryCellList(mesh(), cellIsCut, nbrCellIsCut);
|
|
}
|
|
|
|
forAllConstIters(faceSplitCut_, iter)
|
|
{
|
|
const label facei = iter.key();
|
|
|
|
if (mesh().isInternalFace(facei))
|
|
{
|
|
label own = mesh().faceOwner()[facei];
|
|
label nei = mesh().faceNeighbour()[facei];
|
|
|
|
if (cellLoops_[own].empty() && cellLoops_[nei].empty())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Internal face:" << facei << " cut by " << iter.val()
|
|
<< " has owner:" << own
|
|
<< " and neighbour:" << nei
|
|
<< " that are both uncut"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
label bFacei = facei - mesh().nInternalFaces();
|
|
|
|
label own = mesh().faceOwner()[facei];
|
|
|
|
if (cellLoops_[own].empty() && !nbrCellIsCut[bFacei])
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Boundary face:" << facei << " cut by " << iter.val()
|
|
<< " has owner:" << own
|
|
<< " that is uncut"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::cellCuts::cellCuts
|
|
(
|
|
const polyMesh& mesh,
|
|
const labelList& cutCells,
|
|
const labelList& meshVerts,
|
|
const labelList& meshEdges,
|
|
const scalarField& meshEdgeWeights,
|
|
const bool verbose
|
|
)
|
|
:
|
|
edgeVertex(mesh),
|
|
verbose_(verbose),
|
|
pointIsCut_(expand(mesh.nPoints(), meshVerts)),
|
|
edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
|
|
edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
|
|
faceSplitCut_(cutCells.size()),
|
|
cellLoops_(mesh.nCells()),
|
|
nLoops_(-1),
|
|
cellAnchorPoints_(mesh.nCells())
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : constructor from cut verts and edges" << endl;
|
|
}
|
|
|
|
calcLoopsAndAddressing(cutCells);
|
|
|
|
// Calculate planes and flip cellLoops if necessary
|
|
orientPlanesAndLoops();
|
|
|
|
if (debug)
|
|
{
|
|
check();
|
|
}
|
|
|
|
clearOut();
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : leaving constructor from cut verts and edges"
|
|
<< endl;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::cellCuts::cellCuts
|
|
(
|
|
const polyMesh& mesh,
|
|
const labelList& meshVerts,
|
|
const labelList& meshEdges,
|
|
const scalarField& meshEdgeWeights,
|
|
const bool verbose
|
|
)
|
|
:
|
|
edgeVertex(mesh),
|
|
verbose_(verbose),
|
|
pointIsCut_(expand(mesh.nPoints(), meshVerts)),
|
|
edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
|
|
edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
|
|
faceSplitCut_(mesh.nFaces()/10 + 1),
|
|
cellLoops_(mesh.nCells()),
|
|
nLoops_(-1),
|
|
cellAnchorPoints_(mesh.nCells())
|
|
{
|
|
// Construct from pattern of cuts. Finds out itself which cells are cut.
|
|
// (can go wrong if e.g. all neighbours of cell are refined)
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : constructor from cellLoops" << endl;
|
|
}
|
|
|
|
calcLoopsAndAddressing(identity(mesh.nCells()));
|
|
|
|
// Adds cuts on other side of coupled boundaries
|
|
syncProc();
|
|
|
|
// Calculate planes and flip cellLoops if necessary
|
|
orientPlanesAndLoops();
|
|
|
|
if (debug)
|
|
{
|
|
check();
|
|
}
|
|
|
|
clearOut();
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::cellCuts::cellCuts
|
|
(
|
|
const polyMesh& mesh,
|
|
const labelList& cellLabels,
|
|
const labelListList& cellLoops,
|
|
const List<scalarField>& cellEdgeWeights,
|
|
const bool verbose
|
|
)
|
|
:
|
|
edgeVertex(mesh),
|
|
verbose_(verbose),
|
|
pointIsCut_(mesh.nPoints(), false),
|
|
edgeIsCut_(mesh.nEdges(), false),
|
|
edgeWeight_(mesh.nEdges(), -GREAT),
|
|
faceSplitCut_(cellLabels.size()),
|
|
cellLoops_(mesh.nCells()),
|
|
nLoops_(-1),
|
|
cellAnchorPoints_(mesh.nCells())
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : constructor from cellLoops" << endl;
|
|
}
|
|
|
|
// Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
|
|
// Makes sure cuts are consistent
|
|
setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
|
|
|
|
// Adds cuts on other side of coupled boundaries
|
|
syncProc();
|
|
|
|
// Calculate planes and flip cellLoops if necessary
|
|
orientPlanesAndLoops();
|
|
|
|
if (debug)
|
|
{
|
|
check();
|
|
}
|
|
|
|
clearOut();
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::cellCuts::cellCuts
|
|
(
|
|
const polyMesh& mesh,
|
|
const cellLooper& cellCutter,
|
|
const List<refineCell>& refCells,
|
|
const bool verbose
|
|
)
|
|
:
|
|
edgeVertex(mesh),
|
|
verbose_(verbose),
|
|
pointIsCut_(mesh.nPoints(), false),
|
|
edgeIsCut_(mesh.nEdges(), false),
|
|
edgeWeight_(mesh.nEdges(), -GREAT),
|
|
faceSplitCut_(refCells.size()),
|
|
cellLoops_(mesh.nCells()),
|
|
nLoops_(-1),
|
|
cellAnchorPoints_(mesh.nCells())
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : constructor from cellCutter" << endl;
|
|
}
|
|
|
|
// Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
|
|
// Makes sure cuts are consistent
|
|
setFromCellCutter(cellCutter, refCells);
|
|
|
|
// Adds cuts on other side of coupled boundaries
|
|
syncProc();
|
|
|
|
// Calculate planes and flip cellLoops if necessary
|
|
orientPlanesAndLoops();
|
|
|
|
if (debug)
|
|
{
|
|
check();
|
|
}
|
|
|
|
clearOut();
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : leaving constructor from cellCutter" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::cellCuts::cellCuts
|
|
(
|
|
const polyMesh& mesh,
|
|
const cellLooper& cellCutter,
|
|
const labelList& cellLabels,
|
|
const List<plane>& cutPlanes,
|
|
const bool verbose
|
|
)
|
|
:
|
|
edgeVertex(mesh),
|
|
verbose_(verbose),
|
|
pointIsCut_(mesh.nPoints(), false),
|
|
edgeIsCut_(mesh.nEdges(), false),
|
|
edgeWeight_(mesh.nEdges(), -GREAT),
|
|
faceSplitCut_(cellLabels.size()),
|
|
cellLoops_(mesh.nCells()),
|
|
nLoops_(-1),
|
|
cellAnchorPoints_(mesh.nCells())
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : constructor from cellCutter with prescribed plane"
|
|
<< endl;
|
|
}
|
|
|
|
// Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
|
|
// Makes sure cuts are consistent
|
|
setFromCellCutter(cellCutter, cellLabels, cutPlanes);
|
|
|
|
// Adds cuts on other side of coupled boundaries
|
|
syncProc();
|
|
|
|
// Calculate planes and flip cellLoops if necessary
|
|
orientPlanesAndLoops();
|
|
|
|
if (debug)
|
|
{
|
|
check();
|
|
}
|
|
|
|
clearOut();
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : leaving constructor from cellCutter with prescribed"
|
|
<< " plane" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
Foam::cellCuts::cellCuts
|
|
(
|
|
const polyMesh& mesh,
|
|
const boolList& pointIsCut,
|
|
const boolList& edgeIsCut,
|
|
const scalarField& edgeWeight,
|
|
const Map<edge>& faceSplitCut,
|
|
const labelListList& cellLoops,
|
|
const label nLoops,
|
|
const labelListList& cellAnchorPoints,
|
|
const bool verbose
|
|
)
|
|
:
|
|
edgeVertex(mesh),
|
|
verbose_(verbose),
|
|
pointIsCut_(pointIsCut),
|
|
edgeIsCut_(edgeIsCut),
|
|
edgeWeight_(edgeWeight),
|
|
faceSplitCut_(faceSplitCut),
|
|
cellLoops_(cellLoops),
|
|
nLoops_(nLoops),
|
|
cellAnchorPoints_(cellAnchorPoints)
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "cellCuts : constructor from components" << endl;
|
|
Pout<< "cellCuts : leaving constructor from components" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
|
|
|
Foam::cellCuts::~cellCuts()
|
|
{
|
|
clearOut();
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::clearOut()
|
|
{
|
|
faceCutsPtr_.clear();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
Foam::pointField Foam::cellCuts::loopPoints(const label celli) const
|
|
{
|
|
const labelList& loop = cellLoops_[celli];
|
|
|
|
pointField loopPts(loop.size());
|
|
|
|
forAll(loop, fp)
|
|
{
|
|
const label cut = loop[fp];
|
|
|
|
if (isEdge(cut))
|
|
{
|
|
label edgeI = getEdge(cut);
|
|
|
|
loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
|
|
}
|
|
else
|
|
{
|
|
loopPts[fp] = coord(cut, -GREAT);
|
|
}
|
|
}
|
|
return loopPts;
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::flip(const label celli)
|
|
{
|
|
labelList& loop = cellLoops_[celli];
|
|
|
|
reverse(loop);
|
|
|
|
// Reverse anchor point set.
|
|
cellAnchorPoints_[celli] =
|
|
nonAnchorPoints
|
|
(
|
|
mesh().cellPoints()[celli],
|
|
cellAnchorPoints_[celli],
|
|
loop
|
|
);
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::flipLoopOnly(const label celli)
|
|
{
|
|
labelList& loop = cellLoops_[celli];
|
|
|
|
reverse(loop);
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::writeOBJ
|
|
(
|
|
Ostream& os,
|
|
const pointField& loopPts,
|
|
label& vertI
|
|
) const
|
|
{
|
|
label startVertI = vertI;
|
|
|
|
forAll(loopPts, fp)
|
|
{
|
|
const point& pt = loopPts[fp];
|
|
|
|
os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
|
|
|
|
vertI++;
|
|
}
|
|
|
|
os << 'f';
|
|
forAll(loopPts, fp)
|
|
{
|
|
os << ' ' << startVertI + fp + 1;
|
|
}
|
|
os << endl;
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::writeOBJ(Ostream& os) const
|
|
{
|
|
label vertI = 0;
|
|
|
|
forAll(cellLoops_, celli)
|
|
{
|
|
if (cellLoops_[celli].size())
|
|
{
|
|
writeOBJ(os, loopPoints(celli), vertI);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::cellCuts::writeCellOBJ(const fileName& dir, const label celli) const
|
|
{
|
|
const labelList& anchors = cellAnchorPoints_[celli];
|
|
|
|
writeOBJ(dir, celli, loopPoints(celli), anchors);
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|