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/conversion/polyDualMesh/meshDualiser.C
2019-02-04 21:58:51 +01:00

1482 lines
42 KiB
C++

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / 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/>.
Class
meshDualiser
\*---------------------------------------------------------------------------*/
#include "meshDualiser.H"
#include "meshTools.H"
#include "polyMesh.H"
#include "directTopoChange.H"
#include "mapPolyMesh.H"
#include "edgeFaceCirculator.H"
#include "mergePoints.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::meshDualiser, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::meshDualiser::checkPolyTopoChange(const directTopoChange& meshMod)
{
// Assume no removed points
pointField points(meshMod.points().size());
forAll(meshMod.points(), i)
{
points[i] = meshMod.points()[i];
}
labelList oldToNew;
pointField newPoints;
bool hasMerged = mergePoints
(
points,
1E-6,
false,
oldToNew,
newPoints
);
if (hasMerged)
{
labelListList newToOld(invertOneToMany(newPoints.size(), oldToNew));
forAll(newToOld, newI)
{
if (newToOld[newI].size() != 1)
{
FatalErrorIn
(
"meshDualiser::checkPolyTopoChange(const directTopoChange&)"
) << "duplicate verts:" << newToOld[newI]
<< " coords:"
<< UIndirectList<point>(points, newToOld[newI])()
<< abort(FatalError);
}
}
}
}
// Dump state so far.
void Foam::meshDualiser::dumpPolyTopoChange
(
const directTopoChange& meshMod,
const fileName& prefix
)
{
OFstream str1(prefix + "Faces.obj");
OFstream str2(prefix + "Edges.obj");
Info<< "Dumping current directTopoChange. Faces to " << str1.name()
<< " , points and edges to " << str2.name() << endl;
const DynamicList<point>& points = meshMod.points();
forAll(points, pointI)
{
meshTools::writeOBJ(str1, points[pointI]);
meshTools::writeOBJ(str2, points[pointI]);
}
const DynamicList<face>& faces = meshMod.faces();
forAll(faces, faceI)
{
const face& f = faces[faceI];
str1<< 'f';
forAll(f, fp)
{
str1<< ' ' << f[fp]+1;
}
str1<< nl;
str2<< 'l';
forAll(f, fp)
{
str2<< ' ' << f[fp]+1;
}
str2<< ' ' << f[0]+1 << nl;
}
}
//- Given cell and point on mesh finds the corresponding dualCell. Most
// points become only one cell but the feature points might be split.
Foam::label Foam::meshDualiser::findDualCell
(
const label cellI,
const label pointI
) const
{
const labelList& dualCells = pointToDualCells_[pointI];
if (dualCells.size() == 1)
{
return dualCells[0];
}
else
{
label index = findIndex(mesh_.pointCells()[pointI], cellI);
return dualCells[index];
}
}
// Helper function to generate dualpoints on all boundary edges emanating
// from (boundary & feature) point
void Foam::meshDualiser::generateDualBoundaryEdges
(
const PackedBoolList& isBoundaryEdge,
const label pointI,
directTopoChange& meshMod
)
{
const labelList& pEdges = mesh_.pointEdges()[pointI];
forAll(pEdges, pEdgeI)
{
label edgeI = pEdges[pEdgeI];
if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.get(edgeI) == 1)
{
const edge& e = mesh_.edges()[edgeI];
edgeToDualPoint_[edgeI] = meshMod.addPoint
(
e.centre(mesh_.points()),
pointI, // masterPoint
-1, // zoneID
true // inCell
);
}
}
}
// Return true if point on face has same dual cells on both owner and neighbour
// sides.
bool Foam::meshDualiser::sameDualCell
(
const label faceI,
const label pointI
) const
{
if (!mesh_.isInternalFace(faceI))
{
FatalErrorIn("sameDualCell(const label, const label)")
<< "face:" << faceI << " is not internal face."
<< abort(FatalError);
}
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
return findDualCell(own, pointI) == findDualCell(nei, pointI);
}
Foam::label Foam::meshDualiser::addInternalFace
(
const label masterPointI,
const label masterEdgeI,
const label masterFaceI,
const bool edgeOrder,
const label dualCell0,
const label dualCell1,
const dynamicLabelList& verts,
directTopoChange& meshMod
) const
{
face newFace(verts);
if (edgeOrder != (dualCell0 < dualCell1))
{
reverse(newFace);
}
if (debug)
{
pointField facePoints(meshMod.points(), newFace);
labelList oldToNew;
pointField newPoints;
bool hasMerged = mergePoints
(
facePoints,
1E-6,
false,
oldToNew,
newPoints
);
if (hasMerged)
{
FatalErrorIn("addInternalFace(..)")
<< "verts:" << verts << " newFace:" << newFace
<< " face points:" << facePoints
<< abort(FatalError);
}
}
label zoneID = -1;
bool zoneFlip = false;
if (masterFaceI != -1)
{
zoneID = mesh_.faceZones().whichZone(masterFaceI);
if (zoneID != -1)
{
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
}
}
label dualFaceI;
if (dualCell0 < dualCell1)
{
dualFaceI = meshMod.addFace
(
newFace,
dualCell0, // own
dualCell1, // nei
masterPointI, // masterPointID
masterEdgeI, // masterEdgeID
masterFaceI, // masterFaceID
false, // flipFaceFlux
-1, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
//pointField dualPoints(meshMod.points());
//vector n(newFace.normal(dualPoints));
//n /= mag(n);
//Pout<< "Generated internal dualFace:" << dualFaceI
// << " verts:" << newFace
// << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
// << " n:" << n
// << " between dualowner:" << dualCell0
// << " dualneigbour:" << dualCell1
// << endl;
}
else
{
dualFaceI = meshMod.addFace
(
newFace,
dualCell1, // own
dualCell0, // nei
masterPointI, // masterPointID
masterEdgeI, // masterEdgeID
masterFaceI, // masterFaceID
false, // flipFaceFlux
-1, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
//pointField dualPoints(meshMod.points());
//vector n(newFace.normal(dualPoints));
//n /= mag(n);
//Pout<< "Generated internal dualFace:" << dualFaceI
// << " verts:" << newFace
// << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
// << " n:" << n
// << " between dualowner:" << dualCell1
// << " dualneigbour:" << dualCell0
// << endl;
}
return dualFaceI;
}
Foam::label Foam::meshDualiser::addBoundaryFace
(
const label masterPointI,
const label masterEdgeI,
const label masterFaceI,
const label dualCellI,
const label patchI,
const dynamicLabelList& verts,
directTopoChange& meshMod
) const
{
face newFace(verts);
label zoneID = -1;
bool zoneFlip = false;
if (masterFaceI != -1)
{
zoneID = mesh_.faceZones().whichZone(masterFaceI);
if (zoneID != -1)
{
const faceZone& fZone = mesh_.faceZones()[zoneID];
zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
}
}
label dualFaceI = meshMod.addFace
(
newFace,
dualCellI, // own
-1, // nei
masterPointI, // masterPointID
masterEdgeI, // masterEdgeID
masterFaceI, // masterFaceID
false, // flipFaceFlux
patchI, // patchID
zoneID, // zoneID
zoneFlip // zoneFlip
);
//pointField dualPoints(meshMod.points());
//vector n(newFace.normal(dualPoints));
//n /= mag(n);
//Pout<< "Generated boundary dualFace:" << dualFaceI
// << " verts:" << newFace
// << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
// << " n:" << n
// << " on dualowner:" << dualCellI
// << endl;
return dualFaceI;
}
// Walks around edgeI.
// splitFace=true : creates multiple faces
// splitFace=false: creates single face if same dual cells on both sides,
// multiple faces otherwise.
void Foam::meshDualiser::createFacesAroundEdge
(
const bool splitFace,
const PackedBoolList& isBoundaryEdge,
const label edgeI,
const label startFaceI,
directTopoChange& meshMod,
boolList& doneEFaces
) const
{
const edge& e = mesh_.edges()[edgeI];
const labelList& eFaces = mesh_.edgeFaces()[edgeI];
label fp = edgeFaceCirculator::getMinIndex
(
mesh_.faces()[startFaceI],
e[0],
e[1]
);
edgeFaceCirculator ie
(
mesh_,
startFaceI, // face
true, // ownerSide
fp, // fp
isBoundaryEdge.get(edgeI) == 1 // isBoundaryEdge
);
ie.setCanonical();
bool edgeOrder = ie.sameOrder(e[0], e[1]);
label startFaceLabel = ie.faceLabel();
//Pout<< "At edge:" << edgeI << " verts:" << e
// << " points:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
// << " started walking at face:" << ie.faceLabel()
// << " verts:" << mesh_.faces()[ie.faceLabel()]
// << " edgeOrder:" << edgeOrder
// << " in direction of cell:" << ie.cellLabel()
// << endl;
// Walk and collect face.
dynamicLabelList verts(100);
if (edgeToDualPoint_[edgeI] != -1)
{
verts.append(edgeToDualPoint_[edgeI]);
}
if (faceToDualPoint_[ie.faceLabel()] != -1)
{
doneEFaces[findIndex(eFaces, ie.faceLabel())] = true;
verts.append(faceToDualPoint_[ie.faceLabel()]);
}
if (cellToDualPoint_[ie.cellLabel()] != -1)
{
verts.append(cellToDualPoint_[ie.cellLabel()]);
}
label currentDualCell0 = findDualCell(ie.cellLabel(), e[0]);
label currentDualCell1 = findDualCell(ie.cellLabel(), e[1]);
++ie;
while (true)
{
label faceI = ie.faceLabel();
// Mark face as visited.
doneEFaces[findIndex(eFaces, faceI)] = true;
if (faceToDualPoint_[faceI] != -1)
{
verts.append(faceToDualPoint_[faceI]);
}
label cellI = ie.cellLabel();
if (cellI == -1)
{
// At ending boundary face. We've stored the face point above
// so this is the whole face.
break;
}
label dualCell0 = findDualCell(cellI, e[0]);
label dualCell1 = findDualCell(cellI, e[1]);
// Generate face. (always if splitFace=true; only if needed to
// separate cells otherwise)
if
(
splitFace
|| (
dualCell0 != currentDualCell0
|| dualCell1 != currentDualCell1
)
)
{
// Close current face.
addInternalFace
(
-1, // masterPointI
edgeI, // masterEdgeI
-1, // masterFaceI
edgeOrder,
currentDualCell0,
currentDualCell1,
verts.shrink(),
meshMod
);
// Restart
currentDualCell0 = dualCell0;
currentDualCell1 = dualCell1;
verts.clear();
if (edgeToDualPoint_[edgeI] != -1)
{
verts.append(edgeToDualPoint_[edgeI]);
}
if (faceToDualPoint_[faceI] != -1)
{
verts.append(faceToDualPoint_[faceI]);
}
}
if (cellToDualPoint_[cellI] != -1)
{
verts.append(cellToDualPoint_[cellI]);
}
++ie;
if (ie == ie.end())
{
// Back at start face (for internal edge only). See if this needs
// adding.
if (isBoundaryEdge.get(edgeI) == 0)
{
label startDual = faceToDualPoint_[startFaceLabel];
if (startDual != -1 && findIndex(verts, startDual) == -1)
{
verts.append(startDual);
}
}
break;
}
}
verts.shrink();
addInternalFace
(
-1, // masterPointI
edgeI, // masterEdgeI
-1, // masterFaceI
edgeOrder,
currentDualCell0,
currentDualCell1,
verts,
meshMod
);
}
// Walks around circumference of faceI. Creates single face. Gets given
// starting (feature) edge to start from. Returns ending edge. (all edges
// in form of index in faceEdges)
void Foam::meshDualiser::createFaceFromInternalFace
(
const label faceI,
label& fp,
directTopoChange& meshMod
) const
{
const face& f = mesh_.faces()[faceI];
const labelList& fEdges = mesh_.faceEdges()[faceI];
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
//Pout<< "createFaceFromInternalFace : At face:" << faceI
// << " verts:" << f
// << " points:" << UIndirectList<point>(mesh_.points(), f)()
// << " started walking at edge:" << fEdges[fp]
// << " verts:" << mesh_.edges()[fEdges[fp]]
// << endl;
// Walk and collect face.
dynamicLabelList verts(100);
verts.append(faceToDualPoint_[faceI]);
verts.append(edgeToDualPoint_[fEdges[fp]]);
// Step to vertex after edge mid
fp = f.fcIndex(fp);
// Get cells on either side of face at that point
label currentDualCell0 = findDualCell(own, f[fp]);
label currentDualCell1 = findDualCell(nei, f[fp]);
forAll(f, i)
{
// Check vertex
if (pointToDualPoint_[f[fp]] != -1)
{
verts.append(pointToDualPoint_[f[fp]]);
}
// Edge between fp and fp+1
label edgeI = fEdges[fp];
if (edgeToDualPoint_[edgeI] != -1)
{
verts.append(edgeToDualPoint_[edgeI]);
}
// Next vertex on edge
label nextFp = f.fcIndex(fp);
// Get dual cells on nextFp to check whether face needs closing.
label dualCell0 = findDualCell(own, f[nextFp]);
label dualCell1 = findDualCell(nei, f[nextFp]);
if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
{
// Check: make sure that there is a midpoint on the edge.
if (edgeToDualPoint_[edgeI] == -1)
{
FatalErrorIn("createFacesFromInternalFace(..)")
<< "face:" << faceI << " verts:" << f
<< " points:" << UIndirectList<point>(mesh_.points(), f)()
<< " no feature edge between " << f[fp]
<< " and " << f[nextFp] << " although have different"
<< " dual cells." << endl
<< "point " << f[fp] << " has dual cells "
<< currentDualCell0 << " and " << currentDualCell1
<< " ; point "<< f[nextFp] << " has dual cells "
<< dualCell0 << " and " << dualCell1
<< abort(FatalError);
}
// Close current face.
verts.shrink();
addInternalFace
(
-1, // masterPointI
-1, // masterEdgeI
faceI, // masterFaceI
true, // edgeOrder,
currentDualCell0,
currentDualCell1,
verts,
meshMod
);
break;
}
fp = nextFp;
}
}
// Given a point on a face converts the faces around the point.
// (pointFaces()). Gets starting face and marks off visited faces in donePFaces.
void Foam::meshDualiser::createFacesAroundBoundaryPoint
(
const label patchI,
const label patchPointI,
const label startFaceI,
directTopoChange& meshMod,
boolList& donePFaces // pFaces visited
) const
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const polyPatch& pp = patches[patchI];
const labelList& pFaces = pp.pointFaces()[patchPointI];
const labelList& own = mesh_.faceOwner();
label pointI = pp.meshPoints()[patchPointI];
if (pointToDualPoint_[pointI] == -1)
{
// Not a feature point. Loop over all connected
// pointFaces.
// Starting face
label faceI = startFaceI;
dynamicLabelList verts(4);
while (true)
{
label index = findIndex(pFaces, faceI-pp.start());
// Has face been visited already?
if (donePFaces[index])
{
break;
}
donePFaces[index] = true;
// Insert face centre
verts.append(faceToDualPoint_[faceI]);
label dualCellI = findDualCell(own[faceI], pointI);
// Get the edge before the patchPointI
const face& f = mesh_.faces()[faceI];
label fp = findIndex(f, pointI);
label prevFp = f.rcIndex(fp);
label edgeI = mesh_.faceEdges()[faceI][prevFp];
if (edgeToDualPoint_[edgeI] != -1)
{
verts.append(edgeToDualPoint_[edgeI]);
}
// Get next boundary face (whilst staying on edge).
edgeFaceCirculator circ
(
mesh_,
faceI,
true, // ownerSide
prevFp, // index of edge in face
true // isBoundaryEdge
);
do
{
++circ;
}
while (mesh_.isInternalFace(circ.faceLabel()));
// Step to next face
faceI = circ.faceLabel();
if (faceI < pp.start() || faceI >= pp.start()+pp.size())
{
FatalErrorIn
(
"meshDualiser::createFacesAroundBoundaryPoint(..)"
) << "Walked from face on patch:" << patchI
<< " to face:" << faceI
<< " fc:" << mesh_.faceCentres()[faceI]
<< " on patch:" << patches.whichPatch(faceI)
<< abort(FatalError);
}
// Check if different cell.
if (dualCellI != findDualCell(own[faceI], pointI))
{
FatalErrorIn
(
"meshDualiser::createFacesAroundBoundaryPoint(..)"
) << "Different dual cells but no feature edge"
<< " inbetween point:" << pointI
<< " coord:" << mesh_.points()[pointI]
<< abort(FatalError);
}
}
verts.shrink();
label dualCellI = findDualCell(own[faceI], pointI);
//Bit dodgy: create dualface from the last face (instead of from
// the central point). This will also use the original faceZone to
// put the new face (which might span multiple original faces) in.
addBoundaryFace
(
//pointI, // masterPointI
-1, // masterPointI
-1, // masterEdgeI
faceI, // masterFaceI
dualCellI,
patchI,
verts,
meshMod
);
}
else
{
label faceI = startFaceI;
// Storage for face
dynamicLabelList verts(mesh_.faces()[faceI].size());
// Starting point.
verts.append(pointToDualPoint_[pointI]);
// Find edge between pointI and next point on face.
const labelList& fEdges = mesh_.faceEdges()[faceI];
label nextEdgeI = fEdges[findIndex(mesh_.faces()[faceI], pointI)];
if (edgeToDualPoint_[nextEdgeI] != -1)
{
verts.append(edgeToDualPoint_[nextEdgeI]);
}
do
{
label index = findIndex(pFaces, faceI-pp.start());
// Has face been visited already?
if (donePFaces[index])
{
break;
}
donePFaces[index] = true;
// Face centre
verts.append(faceToDualPoint_[faceI]);
// Find edge before pointI on faceI
const labelList& fEdges = mesh_.faceEdges()[faceI];
const face& f = mesh_.faces()[faceI];
label prevFp = f.rcIndex(findIndex(f, pointI));
label edgeI = fEdges[prevFp];
if (edgeToDualPoint_[edgeI] != -1)
{
// Feature edge. Close any face so far. Note: uses face to
// create dualFace from. Could use pointI instead.
verts.append(edgeToDualPoint_[edgeI]);
addBoundaryFace
(
-1, // masterPointI
-1, // masterEdgeI
faceI, // masterFaceI
findDualCell(own[faceI], pointI),
patchI,
verts.shrink(),
meshMod
);
verts.clear();
verts.append(pointToDualPoint_[pointI]);
verts.append(edgeToDualPoint_[edgeI]);
}
// Cross edgeI to next boundary face
edgeFaceCirculator circ
(
mesh_,
faceI,
true, // ownerSide
prevFp, // index of edge in face
true // isBoundaryEdge
);
do
{
++circ;
}
while (mesh_.isInternalFace(circ.faceLabel()));
// Step to next face. Quit if not on same patch.
faceI = circ.faceLabel();
}
while
(
faceI != startFaceI
&& faceI >= pp.start()
&& faceI < pp.start()+pp.size()
);
if (verts.size() > 2)
{
// Note: face created from face, not from pointI
addBoundaryFace
(
-1, // masterPointI
-1, // masterEdgeI
startFaceI, // masterFaceI
findDualCell(own[faceI], pointI),
patchI,
verts.shrink(),
meshMod
);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::meshDualiser::meshDualiser(const polyMesh& mesh)
:
mesh_(mesh),
pointToDualCells_(mesh_.nPoints()),
pointToDualPoint_(mesh_.nPoints(), -1),
cellToDualPoint_(mesh_.nCells()),
faceToDualPoint_(mesh_.nFaces(), -1),
edgeToDualPoint_(mesh_.nEdges(), -1)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::meshDualiser::setRefinement
(
const bool splitFace,
const labelList& featureFaces,
const labelList& featureEdges,
const labelList& singleCellFeaturePoints,
const labelList& multiCellFeaturePoints,
directTopoChange& meshMod
)
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const vectorField& cellCentres = mesh_.cellCentres();
// Mark boundary edges and points.
// (Note: in 1.4.2 we can use the built-in mesh point ordering
// facility instead)
PackedBoolList isBoundaryEdge(mesh_.nEdges());
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
{
const labelList& fEdges = mesh_.faceEdges()[faceI];
forAll(fEdges, i)
{
isBoundaryEdge.set(fEdges[i], 1);
}
}
if (splitFace)
{
// This is a special mode where whenever we are walking around an edge
// every area through a cell becomes a separate dualface. So two
// dual cells will probably have more than one dualface between them!
// This mode implies that
// - all faces have to be feature faces since there has to be a
// dualpoint at the face centre.
// - all edges have to be feature edges ,,
boolList featureFaceSet(mesh_.nFaces(), false);
forAll(featureFaces, i)
{
featureFaceSet[featureFaces[i]] = true;
}
label faceI = findIndex(featureFaceSet, false);
if (faceI != -1)
{
FatalErrorIn
(
"meshDualiser::setRefinement"
"(const labelList&, const labelList&, const labelList&"
", const labelList, directTopoChange&)"
) << "In split-face-mode (splitFace=true) but not all faces"
<< " marked as feature faces." << endl
<< "First conflicting face:" << faceI
<< " centre:" << mesh_.faceCentres()[faceI]
<< abort(FatalError);
}
boolList featureEdgeSet(mesh_.nEdges(), false);
forAll(featureEdges, i)
{
featureEdgeSet[featureEdges[i]] = true;
}
label edgeI = findIndex(featureEdgeSet, false);
if (edgeI != -1)
{
const edge& e = mesh_.edges()[edgeI];
FatalErrorIn
(
"meshDualiser::setRefinement"
"(const labelList&, const labelList&, const labelList&"
", const labelList, directTopoChange&)"
) << "In split-face-mode (splitFace=true) but not all edges"
<< " marked as feature edges." << endl
<< "First conflicting edge:" << edgeI
<< " verts:" << e
<< " coords:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
<< abort(FatalError);
}
}
else
{
// Check that all boundary faces are feature faces.
boolList featureFaceSet(mesh_.nFaces(), false);
forAll(featureFaces, i)
{
featureFaceSet[featureFaces[i]] = true;
}
for
(
label faceI = mesh_.nInternalFaces();
faceI < mesh_.nFaces();
faceI++
)
{
if (!featureFaceSet[faceI])
{
FatalErrorIn
(
"meshDualiser::setRefinement"
"(const labelList&, const labelList&, const labelList&"
", const labelList, directTopoChange&)"
) << "Not all boundary faces marked as feature faces."
<< endl
<< "First conflicting face:" << faceI
<< " centre:" << mesh_.faceCentres()[faceI]
<< abort(FatalError);
}
}
}
// Start creating cells, points, and faces (in that order)
// 1. Mark which cells to create
// Mostly every point becomes one cell but sometimes (for feature points)
// all cells surrounding a feature point become cells. Also a non-manifold
// point can create two cells! So a dual cell is uniquely defined by a
// mesh point + cell (as in pointCells index)
// 2. Mark which face centres to create
// 3. Internal faces can now consist of
// - only cell centres of walk around edge
// - cell centres + face centres of walk around edge
// - same but now other side is not a single cell
// 4. Boundary faces (or internal faces between cell zones!) now consist of
// - walk around boundary point.
autoPtr<OFstream> dualCcStr;
if (debug)
{
dualCcStr.reset(new OFstream("dualCc.obj"));
Pout<< "Dumping centres of dual cells to " << dualCcStr().name()
<< endl;
}
// Dual cells (from points)
// ~~~~~~~~~~~~~~~~~~~~~~~~
// pointToDualCells_[pointI]
// - single entry : all cells surrounding point all become the same
// cell.
// - multiple entries: in order of pointCells.
// feature points that become single cell
forAll(singleCellFeaturePoints, i)
{
label pointI = singleCellFeaturePoints[i];
pointToDualPoint_[pointI] = meshMod.addPoint
(
mesh_.points()[pointI],
pointI, // masterPoint
mesh_.pointZones().whichZone(pointI), // zoneID
true // inCell
);
// Generate single cell
pointToDualCells_[pointI].setSize(1);
pointToDualCells_[pointI][0] = meshMod.addCell
(
pointI, //masterPointID,
-1, //masterEdgeID,
-1, //masterFaceID,
-1, //masterCellID,
-1 //zoneID
);
if (dualCcStr.valid())
{
meshTools::writeOBJ(dualCcStr(), mesh_.points()[pointI]);
}
}
// feature points that become multiple cells
forAll(multiCellFeaturePoints, i)
{
label pointI = multiCellFeaturePoints[i];
if (pointToDualCells_[pointI].size() > 0)
{
FatalErrorIn
(
"meshDualiser::setRefinement"
"(const labelList&, const labelList&, const labelList&"
", const labelList, directTopoChange&)"
) << "Point " << pointI << " at:" << mesh_.points()[pointI]
<< " is both in singleCellFeaturePoints"
<< " and multiCellFeaturePoints."
<< abort(FatalError);
}
pointToDualPoint_[pointI] = meshMod.addPoint
(
mesh_.points()[pointI],
pointI, // masterPoint
mesh_.pointZones().whichZone(pointI), // zoneID
true // inCell
);
// Create dualcell for every cell connected to dual point
const labelList& pCells = mesh_.pointCells()[pointI];
pointToDualCells_[pointI].setSize(pCells.size());
forAll(pCells, pCellI)
{
pointToDualCells_[pointI][pCellI] = meshMod.addCell
(
pointI, //masterPointID
-1, //masterEdgeID
-1, //masterFaceID
-1, //masterCellID
mesh_.cellZones().whichZone(pCells[pCellI]) //zoneID
);
if (dualCcStr.valid())
{
meshTools::writeOBJ
(
dualCcStr(),
0.5*(mesh_.points()[pointI]+cellCentres[pCells[pCellI]])
);
}
}
}
// Normal points
forAll(mesh_.points(), pointI)
{
if (pointToDualCells_[pointI].empty())
{
pointToDualCells_[pointI].setSize(1);
pointToDualCells_[pointI][0] = meshMod.addCell
(
pointI, //masterPointID,
-1, //masterEdgeID,
-1, //masterFaceID,
-1, //masterCellID,
-1 //zoneID
);
if (dualCcStr.valid())
{
meshTools::writeOBJ(dualCcStr(), mesh_.points()[pointI]);
}
}
}
// Dual points (from cell centres, feature faces, feature edges)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(cellToDualPoint_, cellI)
{
cellToDualPoint_[cellI] = meshMod.addPoint
(
cellCentres[cellI],
mesh_.faces()[mesh_.cells()[cellI][0]][0], // masterPoint
-1, // zoneID
true // inCell
);
}
// From face to dual point
forAll(featureFaces, i)
{
label faceI = featureFaces[i];
faceToDualPoint_[faceI] = meshMod.addPoint
(
mesh_.faceCentres()[faceI],
mesh_.faces()[faceI][0], // masterPoint
-1, // zoneID
true // inCell
);
}
// Detect whether different dual cells on either side of a face. This
// would neccesitate having a dual face built from the face and thus a
// dual point at the face centre.
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
if (faceToDualPoint_[faceI] == -1)
{
const face& f = mesh_.faces()[faceI];
forAll(f, fp)
{
label ownDualCell = findDualCell(own[faceI], f[fp]);
label neiDualCell = findDualCell(nei[faceI], f[fp]);
if (ownDualCell != neiDualCell)
{
faceToDualPoint_[faceI] = meshMod.addPoint
(
mesh_.faceCentres()[faceI],
f[fp], // masterPoint
-1, // zoneID
true // inCell
);
break;
}
}
}
}
// From edge to dual point
forAll(featureEdges, i)
{
label edgeI = featureEdges[i];
const edge& e = mesh_.edges()[edgeI];
edgeToDualPoint_[edgeI] = meshMod.addPoint
(
e.centre(mesh_.points()),
e[0], // masterPoint
-1, // zoneID
true // inCell
);
}
// Detect whether different dual cells on either side of an edge. This
// would neccesitate having a dual face built perpendicular to the edge
// and thus a dual point at the mid of the edge.
// Note: not really true - the face can be built without the edge centre!
const labelListList& edgeCells = mesh_.edgeCells();
forAll(edgeCells, edgeI)
{
if (edgeToDualPoint_[edgeI] == -1)
{
const edge& e = mesh_.edges()[edgeI];
// We need a point on the edge if not all cells on both sides
// are the same.
const labelList& eCells = mesh_.edgeCells()[edgeI];
label dualE0 = findDualCell(eCells[0], e[0]);
label dualE1 = findDualCell(eCells[0], e[1]);
for (label i = 1; i < eCells.size(); i++)
{
label newDualE0 = findDualCell(eCells[i], e[0]);
if (dualE0 != newDualE0)
{
edgeToDualPoint_[edgeI] = meshMod.addPoint
(
e.centre(mesh_.points()),
e[0], // masterPoint
mesh_.pointZones().whichZone(e[0]), // zoneID
true // inCell
);
break;
}
label newDualE1 = findDualCell(eCells[i], e[1]);
if (dualE1 != newDualE1)
{
edgeToDualPoint_[edgeI] = meshMod.addPoint
(
e.centre(mesh_.points()),
e[1], // masterPoint
mesh_.pointZones().whichZone(e[1]), // zoneID
true // inCell
);
break;
}
}
}
}
// Make sure all boundary edges emanating from feature points are
// feature edges as well.
forAll(singleCellFeaturePoints, i)
{
generateDualBoundaryEdges
(
isBoundaryEdge,
singleCellFeaturePoints[i],
meshMod
);
}
forAll(multiCellFeaturePoints, i)
{
generateDualBoundaryEdges
(
isBoundaryEdge,
multiCellFeaturePoints[i],
meshMod
);
}
// Check for duplicate points
if (debug)
{
dumpPolyTopoChange(meshMod, "generatedPoints_");
checkPolyTopoChange(meshMod);
}
// Now we have all points and cells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - pointToDualCells_ : per point a single dualCell or multiple dualCells
// - pointToDualPoint_ : per point -1 or the dual point at the coordinate
// - edgeToDualPoint_ : per edge -1 or the edge centre
// - faceToDualPoint_ : per face -1 or the face centre
// - cellToDualPoint_ : per cell the cell centre
// Now we have to walk all edges and construct faces. Either single face
// per edge or multiple (-if nonmanifold edge -if different dualcells)
const edgeList& edges = mesh_.edges();
forAll(edges, edgeI)
{
const labelList& eFaces = mesh_.edgeFaces()[edgeI];
boolList doneEFaces(eFaces.size(), false);
forAll(eFaces, i)
{
if (!doneEFaces[i])
{
// We found a face that hasn't yet been visited. This might
// happen for non-manifold edges where a single edge can
// become multiple faces.
label startFaceI = eFaces[i];
//Pout<< "Walking edge:" << edgeI
// << " points:" << mesh_.points()[e[0]]
// << mesh_.points()[e[1]]
// << " startFace:" << startFaceI
// << " at:" << mesh_.faceCentres()[startFaceI]
// << endl;
createFacesAroundEdge
(
splitFace,
isBoundaryEdge,
edgeI,
startFaceI,
meshMod,
doneEFaces
);
}
}
}
if (debug)
{
dumpPolyTopoChange(meshMod, "generatedFacesFromEdges_");
}
// Create faces from feature faces. These can be internal or external faces.
// - feature face : centre needs to be included.
// - single cells on either side: triangulate
// - multiple cells: create single face between unique cell pair. Only
// create face where cells differ on either side.
// - non-feature face : inbetween cell zones.
forAll(faceToDualPoint_, faceI)
{
if (faceToDualPoint_[faceI] != -1 && mesh_.isInternalFace(faceI))
{
const face& f = mesh_.faces()[faceI];
const labelList& fEdges = mesh_.faceEdges()[faceI];
// Starting edge
label fp = 0;
do
{
// Find edge that is in dual mesh and where the
// next point (fp+1) has different dual cells on either side.
bool foundStart = false;
do
{
if
(
edgeToDualPoint_[fEdges[fp]] != -1
&& !sameDualCell(faceI, f.nextLabel(fp))
)
{
foundStart = true;
break;
}
fp = f.fcIndex(fp);
}
while (fp != 0);
if (!foundStart)
{
break;
}
// Walk from edge fp and generate a face.
createFaceFromInternalFace
(
faceI,
fp,
meshMod
);
}
while(fp != 0);
}
}
if (debug)
{
dumpPolyTopoChange(meshMod, "generatedFacesFromFeatFaces_");
}
// Create boundary faces. Every boundary point has one or more dualcells.
// These need to be closed.
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelListList& pointFaces = pp.pointFaces();
forAll(pointFaces, patchPointI)
{
const labelList& pFaces = pointFaces[patchPointI];
boolList donePFaces(pFaces.size(), false);
forAll(pFaces, i)
{
if (!donePFaces[i])
{
// Starting face
label startFaceI = pp.start()+pFaces[i];
//Pout<< "Walking around point:" << pointI
// << " coord:" << mesh_.points()[pointI]
// << " on patch:" << patchI
// << " startFace:" << startFaceI
// << " at:" << mesh_.faceCentres()[startFaceI]
// << endl;
createFacesAroundBoundaryPoint
(
patchI,
patchPointI,
startFaceI,
meshMod,
donePFaces // pFaces visited
);
}
}
}
}
if (debug)
{
dumpPolyTopoChange(meshMod, "generatedFacesFromBndFaces_");
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //