diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/files b/applications/utilities/mesh/advanced/collapseEdges/Make/files
index d89ca6e737c6c84ce878bd4e57566fbd070333a4..a15838abe84da087d0c63693617ed5b410e56b22 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/files
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/files
@@ -1,3 +1,4 @@
collapseEdges.C
+pointEdgeCollapse/pointEdgeCollapse.C
EXE = $(FOAM_APPBIN)/collapseEdges
diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/options b/applications/utilities/mesh/advanced/collapseEdges/Make/options
index 21b17b14c9d7dd8ae8c0c55fdcf698ce54019fcb..d1efa61fd56aef52623d7e352c300f8948503fd2 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/options
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/options
@@ -1,6 +1,10 @@
EXE_INC = \
- -I$(LIB_SRC)/dynamicMesh/lnInclude
+ -I$(LIB_SRC)/dynamicMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -IpointEdgeCollapse
EXE_LIBS = \
-ldynamicMesh \
- -lmeshTools
+ -lmeshTools \
+ -lfiniteVolume
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index f1eb28e5947759fe22aee59bf0c28a3ed3a2c0f0..899159b2a29ba4005f6fdd740260f1fee794674b 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -41,170 +41,634 @@ Description
#include "argList.H"
#include "Time.H"
-#include "edgeCollapser.H"
#include "polyTopoChange.H"
#include "polyTopoChanger.H"
#include "polyMesh.H"
+#include "fvMesh.H"
#include "mapPolyMesh.H"
#include "mathematicalConstants.H"
#include "PackedBoolList.H"
#include "SortableList.H"
#include "unitConversion.H"
+#include "globalMeshData.H"
+#include "globalIndex.H"
+#include "PointEdgeWave.H"
+#include "pointEdgeCollapse.H"
+#include "motionSmoother.H"
+
+#include "OFstream.H"
+#include "meshTools.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Get faceEdges in order of face points, i.e. faceEdges[0] is between
-// f[0] and f[1]
-labelList getSortedEdges
+label findIndex
(
- const edgeList& edges,
- const labelList& f,
- const labelList& edgeLabels
+ const labelList& elems,
+ const label start,
+ const label size,
+ const label val
)
{
- labelList faceEdges(edgeLabels.size(), -1);
+ for (label i = start; i < size; i++)
+ {
+ if (elems[i] == val)
+ {
+ return i;
+ }
+ }
+ return -1;
+}
- // Find starting pos in f for every edgeLabels
- forAll(edgeLabels, i)
+
+void filterFace
+(
+ const label faceI,
+ face& f,
+ const List<pointEdgeCollapse>& allPointInfo,
+ const Map<DynamicList<label> >& collapseStrings
+)
+{
+ label newFp = 0;
+
+ face oldFace = f;
+
+ forAll(f, fp)
{
- label edgeI = edgeLabels[i];
+ label pointI = f[fp];
- const edge& e = edges[edgeI];
+ label collapsePoint = allPointInfo[pointI].collapseIndex();
- label fp = findIndex(f, e[0]);
- label fp1 = f.fcIndex(fp);
+ if (collapseStrings.found(collapsePoint))
+ {
+ collapsePoint = collapseStrings[collapsePoint][0];
+ }
- if (f[fp1] == e[1])
+ if (collapsePoint == -1)
{
- // EdgeI between fp -> fp1
- faceEdges[fp] = edgeI;
+ WarningIn
+ (
+ "filterFace"
+ "(const label, face&, const List<pointEdgeCollapse>&)"
+ )
+ << "Point " << pointI << " was not visited by PointEdgeWave"
+ << endl;
+ }
+ else if (collapsePoint == -2)
+ {
+ f[newFp++] = pointI;
}
else
{
- // EdgeI between fp-1 -> fp
- faceEdges[f.rcIndex(fp)] = edgeI;
+ if (findIndex(f, 0, newFp, collapsePoint) == -1)
+ {
+ f[newFp++] = collapsePoint;
+ }
}
}
- return faceEdges;
+
+ // Check for pinched face. Tries to correct
+ // - consecutive duplicate vertex. Removes duplicate vertex.
+ // - duplicate vertex with one other vertex in between (spike).
+ // Both of these should not really occur! and should be checked before
+ // collapsing edges.
+
+ const label size = newFp;
+
+ newFp = 2;
+
+ for (label fp = 2; fp < size; fp++)
+ {
+ label fp1 = fp-1;
+ label fp2 = fp-2;
+
+ label pointI = f[fp];
+
+ // Search for previous occurrence.
+ label index = findIndex(f, 0, fp, pointI);
+
+ if (index == fp1)
+ {
+ WarningIn
+ (
+ "Foam::edgeCollapser::filterFace(const label faceI, "
+ "face& f) const"
+ ) << "Removing consecutive duplicate vertex in face "
+ << f << endl;
+ // Don't store current pointI
+ }
+ else if (index == fp2)
+ {
+ WarningIn
+ (
+ "Foam::edgeCollapser::filterFace(const label faceI, "
+ "face& f) const"
+ ) << "Removing non-consecutive duplicate vertex in face "
+ << f << endl;
+ // Don't store current pointI and remove previous
+ newFp--;
+ }
+ else if (index != -1)
+ {
+ WarningIn
+ (
+ "Foam::edgeCollapser::filterFace(const label faceI, "
+ "face& f) const"
+ ) << "Pinched face " << f << endl;
+ f[newFp++] = pointI;
+ }
+ else
+ {
+ f[newFp++] = pointI;
+ }
+ }
+
+ f.setSize(newFp);
}
-// Merges edges which are in straight line. I.e. edge split by point.
-label mergeEdges
+bool setRefinement
(
const polyMesh& mesh,
- const scalar maxCos,
- edgeCollapser& collapser
+ polyTopoChange& meshMod,
+ const List<pointEdgeCollapse>& allPointInfo
)
{
- const pointField& points = mesh.points();
- const edgeList& edges = mesh.edges();
- const labelListList& pointEdges = mesh.pointEdges();
- const labelList& region = collapser.pointRegion();
- const labelList& master = collapser.pointRegionMaster();
+ const cellList& cells = mesh.cells();
+ const labelList& faceOwner = mesh.faceOwner();
+ const labelList& faceNeighbour = mesh.faceNeighbour();
+ const labelListList& pointFaces = mesh.pointFaces();
+ const pointZoneMesh& pointZones = mesh.pointZones();
+ globalIndex globalStrings(mesh.nPoints());
+
+ boolList removedPoints(mesh.nPoints(), false);
+
+ // Create strings of edges
+ Map<DynamicList<label> > collapseStrings;
+
+ forAll(allPointInfo, pointI)
+ {
+ const label collapseIndex = allPointInfo[pointI].collapseIndex();
+
+ if (collapseIndex != -1 && collapseIndex != -2)
+ {
+ collapseStrings(collapseIndex).append(pointI);
+ }
+ }
+
+ bool meshChanged = false;
+
+ // Current faces (is also collapseStatus: f.size() < 3)
+ faceList newFaces(mesh.faces());
+
+ // Current cellCollapse status
+ boolList cellRemoved(mesh.nCells(), false);
+
+ label nUnvisited = 0;
+ label nUncollapsed = 0;
label nCollapsed = 0;
- forAll(pointEdges, pointI)
+ forAll(allPointInfo, pI)
{
- const labelList& pEdges = pointEdges[pointI];
+ const pointEdgeCollapse& pec = allPointInfo[pI];
- if (pEdges.size() == 2)
+ if (pec.collapseIndex() == -1)
+ {
+ nUnvisited++;
+ }
+ else if (pec.collapseIndex() == -2)
{
- const edge& leftE = edges[pEdges[0]];
- const edge& rightE = edges[pEdges[1]];
+ nUncollapsed++;
+ }
+ else if (pec.collapseIndex() >= 0)
+ {
+ nCollapsed++;
+ }
+ }
+
+ label nPoints = allPointInfo.size();
+
+ reduce(nPoints, sumOp<label>());
+ reduce(nUnvisited, sumOp<label>());
+ reduce(nUncollapsed, sumOp<label>());
+ reduce(nCollapsed, sumOp<label>());
- // Get the two vertices on both sides of the point
- label leftV = leftE.otherVertex(pointI);
- label rightV = rightE.otherVertex(pointI);
+ Info<< incrIndent;
+ Info<< indent << "Number of points : " << nPoints << nl
+ << indent << "Not visited : " << nUnvisited << nl
+ << indent << "Not collapsed : " << nUncollapsed << nl
+ << indent << "Collapsed : " << nCollapsed << nl
+ << endl;
+ Info<< decrIndent;
- // Collapse only if none of the points part of merge network
- // or all of networks with different masters.
- label midMaster = -1;
- if (region[pointI] != -1)
+ do
+ {
+ // Update face collapse from edge collapses
+ forAll(newFaces, faceI)
+ {
+ filterFace(faceI, newFaces[faceI], allPointInfo, collapseStrings);
+ }
+
+ // Check if faces to be collapsed cause cells to become collapsed.
+ label nCellCollapsed = 0;
+
+ forAll(cells, cellI)
+ {
+ if (!cellRemoved[cellI])
{
- midMaster = master[region[pointI]];
+ const cell& cFaces = cells[cellI];
+
+ label nFaces = cFaces.size();
+
+ forAll(cFaces, i)
+ {
+ label faceI = cFaces[i];
+
+ if (newFaces[faceI].size() < 3)
+ {
+ --nFaces;
+
+ if (nFaces < 4)
+ {
+ Info<< "Cell:" << cellI
+ << " uses faces:" << cFaces
+ << " of which too many are marked for removal:"
+ << endl
+ << " ";
+ forAll(cFaces, j)
+ {
+ if (newFaces[cFaces[j]].size() < 3)
+ {
+ Info<< ' '<< cFaces[j];
+ }
+ }
+ Info<< endl;
+
+ cellRemoved[cellI] = true;
+
+ // Collapse all edges of cell to nothing
+ //collapseEdges(cellEdges[cellI]);
+
+ nCellCollapsed++;
+
+ break;
+ }
+ }
+ }
}
+ }
+
+ if (nCellCollapsed == 0)
+ {
+ break;
+ }
+ } while (true);
+
+
+ // Keep track of faces that have been done already.
+ boolList doneFace(mesh.nFaces(), false);
+
+ {
+ // Mark points used.
+ boolList usedPoint(mesh.nPoints(), false);
- label leftMaster = -2;
- if (region[leftV] != -1)
+ forAll(cellRemoved, cellI)
+ {
+ if (cellRemoved[cellI])
{
- leftMaster = master[region[leftV]];
+ meshMod.removeCell(cellI, -1);
}
+ }
+
+ // Remove faces
+ forAll(newFaces, faceI)
+ {
+ const face& f = newFaces[faceI];
+
+ if (f.size() < 3)
+ {
+ meshMod.removeFace(faceI, -1);
+ meshChanged = true;
- label rightMaster = -3;
- if (region[rightV] != -1)
+ // Mark face as been done.
+ doneFace[faceI] = true;
+ }
+ else
{
- rightMaster = master[region[rightV]];
+ // Kept face. Mark vertices
+ forAll(f, fp)
+ {
+ usedPoint[f[fp]] = true;
+ }
}
+ }
- if
- (
- midMaster != leftMaster
- && midMaster != rightMaster
- && leftMaster != rightMaster
- )
+ // Remove unused vertices that have not been marked for removal already
+ forAll(usedPoint, pointI)
+ {
+ if (!usedPoint[pointI])
{
- // Check if the two edge are in line
- vector leftVec = points[pointI] - points[leftV];
- leftVec /= mag(leftVec) + VSMALL;
+ removedPoints[pointI] = true;
+ meshMod.removePoint(pointI, -1);
+ meshChanged = true;
+ }
+ }
+ }
+
+ // Modify the point location of the remaining points
+ forAll(allPointInfo, pointI)
+ {
+ const label collapseIndex = allPointInfo[pointI].collapseIndex();
+ const point& collapsePoint = allPointInfo[pointI].collapsePoint();
+
+ if
+ (
+ removedPoints[pointI] == false
+ && collapseIndex != -1
+ && collapseIndex != -2
+ )
+ {
+ meshMod.modifyPoint
+ (
+ pointI,
+ collapsePoint,
+ pointZones.whichZone(pointI),
+ false
+ );
+ }
+ }
+
+
+ const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();
+ const faceZoneMesh& faceZones = mesh.faceZones();
+
+ // Renumber faces that use points
+ forAll(allPointInfo, pointI)
+ {
+ if (removedPoints[pointI] == true)
+ {
+ const labelList& changedFaces = pointFaces[pointI];
- vector rightVec = points[rightV] - points[pointI];
- rightVec /= mag(rightVec) + VSMALL;
+ forAll(changedFaces, changedFaceI)
+ {
+ label faceI = changedFaces[changedFaceI];
- if ((leftVec & rightVec) > maxCos)
+ if (!doneFace[faceI])
{
- // Collapse one (left) side of the edge. Make left vertex
- // the master.
- //if (collapser.unaffectedEdge(pEdges[0]))
+ doneFace[faceI] = true;
+
+ // Get current zone info
+ label zoneID = faceZones.whichZone(faceI);
+
+ bool zoneFlip = false;
+
+ if (zoneID >= 0)
{
- collapser.collapseEdge(pEdges[0], leftV);
- nCollapsed++;
+ const faceZone& fZone = faceZones[zoneID];
+
+ zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+ }
+
+ // Get current connectivity
+ label own = faceOwner[faceI];
+ label nei = -1;
+ label patchID = -1;
+
+ if (mesh.isInternalFace(faceI))
+ {
+ nei = faceNeighbour[faceI];
+ }
+ else
+ {
+ patchID = boundaryMesh.whichPatch(faceI);
}
+
+ meshMod.modifyFace
+ (
+ newFaces[faceI], // face
+ faceI, // faceI to change
+ own, // owner
+ nei, // neighbour
+ false, // flipFaceFlux
+ patchID, // patch
+ zoneID,
+ zoneFlip
+ );
+ meshChanged = true;
}
}
}
}
- return nCollapsed;
+ // Print regions:
+// printRegions();
+
+ return meshChanged;
+}
+
+
+// Get faceEdges in order of face points, i.e. faceEdges[0] is between
+// f[0] and f[1]
+labelList getSortedEdges
+(
+ const edgeList& edges,
+ const labelList& f,
+ const labelList& edgeLabels
+)
+{
+ labelList faceEdges(edgeLabels.size(), -1);
+
+ // Find starting pos in f for every edgeLabels
+ forAll(edgeLabels, i)
+ {
+ label edgeI = edgeLabels[i];
+
+ const edge& e = edges[edgeI];
+
+ label fp = findIndex(f, e[0]);
+ label fp1 = f.fcIndex(fp);
+
+ if (f[fp1] == e[1])
+ {
+ // EdgeI between fp -> fp1
+ faceEdges[fp] = edgeI;
+ }
+ else
+ {
+ // EdgeI between fp-1 -> fp
+ faceEdges[f.rcIndex(fp)] = edgeI;
+ }
+ }
+
+ return faceEdges;
}
+// Merges edges which are in straight line. I.e. edge split by point.
+//label mergeEdges
+//(
+// const polyMesh& mesh,
+// const scalar maxCos,
+// List<pointEdgeCollapse>& allPointInfo
+//)
+//{
+// const pointField& points = mesh.points();
+// const edgeList& edges = mesh.edges();
+// const labelListList& pointEdges = mesh.pointEdges();
+// const labelList& region = collapser.pointRegion();
+// const labelList& master = collapser.pointRegionMaster();
+//
+// label nCollapsed = 0;
+//
+// forAll(pointEdges, pointI)
+// {
+// const labelList& pEdges = pointEdges[pointI];
+//
+// if (pEdges.size() == 2)
+// {
+// const edge& leftE = edges[pEdges[0]];
+// const edge& rightE = edges[pEdges[1]];
+//
+// // Get the two vertices on both sides of the point
+// label leftV = leftE.otherVertex(pointI);
+// label rightV = rightE.otherVertex(pointI);
+//
+// // Collapse only if none of the points part of merge network
+// // or all of networks with different masters.
+// label midMaster = -1;
+// if (region[pointI] != -1)
+// {
+// midMaster = master[region[pointI]];
+// }
+//
+// label leftMaster = -2;
+// if (region[leftV] != -1)
+// {
+// leftMaster = master[region[leftV]];
+// }
+//
+// label rightMaster = -3;
+// if (region[rightV] != -1)
+// {
+// rightMaster = master[region[rightV]];
+// }
+//
+// if
+// (
+// midMaster != leftMaster
+// && midMaster != rightMaster
+// && leftMaster != rightMaster
+// )
+// {
+// // Check if the two edge are in line
+// vector leftVec = points[pointI] - points[leftV];
+// leftVec /= mag(leftVec) + VSMALL;
+//
+// vector rightVec = points[rightV] - points[pointI];
+// rightVec /= mag(rightVec) + VSMALL;
+//
+// if ((leftVec & rightVec) > maxCos)
+// {
+// // Collapse one (left) side of the edge. Make left vertex
+// // the master.
+// //if (collapser.unaffectedEdge(pEdges[0]))
+// const edge& e = mesh.edges()[pEdges[0]];
+//
+// if
+// (
+// allPointInfo[e[0]].collapseIndex() < 0
+// && allPointInfo[e[1]].collapseIndex() < 0
+// )
+// {
+// //pointEdgeCollapse pec(mesh.points()[leftV], leftV);
+//
+// allPointInfo[e[0]] = pec;
+// allPointInfo[e[1]] = pec;
+//
+// //collapser.collapseEdge(pEdges[0], leftV);
+// nCollapsed++;
+// }
+// }
+// }
+// }
+// }
+//
+// return nCollapsed;
+//}
+
+
// Return master point edge needs to be collapsed to (or -1)
-label edgeMaster(const PackedBoolList& boundaryPoint, const edge& e)
+label edgeMaster
+(
+ const labelList& boundaryPoint,
+ const bool flipEdge,
+ const edge& e
+)
{
label masterPoint = -1;
+ label e0 = e[0];
+ label e1 = e[1];
+
+ if (flipEdge)
+ {
+ e0 = e[1];
+ e1 = e[0];
+ }
+
+ // Check if one of the points is on a processor
+// if
+// (
+// boundaryPoint[e0] > 0
+// && boundaryPoint[e1] > 0
+// )
+// {
+// if (boundaryPoint[e0] != boundaryPoint[e1])
+// {
+// return -1;
+// }
+// }
+//
+// if (boundaryPoint[e0] > 0)
+// {
+// return e0;
+// }
+// else if (boundaryPoint[e1] > 0)
+// {
+// return e1;
+// }
+
// Collapse edge to boundary point.
- if (boundaryPoint.get(e[0]))
+ if (boundaryPoint[e0] == 0)
{
- if (boundaryPoint.get(e[1]))
+ if (boundaryPoint[e1] == 0)
{
// Both points on boundary. Choose one to collapse to.
// Note: should look at feature edges/points!
- masterPoint = e[0];
+ masterPoint = e0;
}
else
{
- masterPoint = e[0];
+ masterPoint = e0;
}
}
else
{
- if (boundaryPoint.get(e[1]))
+ if (boundaryPoint[e1] == 0)
{
- masterPoint = e[1];
+ masterPoint = e1;
}
else
{
// None on boundary. Choose arbitrary.
// Note: should look at geometry?
- masterPoint = e[0];
+ masterPoint = e0;
}
}
+
return masterPoint;
}
@@ -212,124 +676,266 @@ label edgeMaster(const PackedBoolList& boundaryPoint, const edge& e)
label collapseSmallEdges
(
const polyMesh& mesh,
- const PackedBoolList& boundaryPoint,
+ const scalarList& freezeEdges,
+ const labelList& boundaryPoint,
const scalar minLen,
- edgeCollapser& collapser
+ List<pointEdgeCollapse>& allPointInfo
)
{
const pointField& points = mesh.points();
const edgeList& edges = mesh.edges();
- // Collapse all edges that are too small. Choose intelligently which
- // point to collapse edge to.
-
- label nCollapsed = 0;
+ // Store collapse direction in collapseEdge
+ // -1 -> Do not collapse
+ // 0 -> Collapse to start point
+ // 1 -> Collapse to end point
+ labelList collapseEdge(edges.size(), -1);
forAll(edges, edgeI)
{
const edge& e = edges[edgeI];
- if (e.mag(points) < minLen)
+ if (e.mag(points) < minLen*freezeEdges[edgeI])
{
- label master = edgeMaster(boundaryPoint, e);
-
- if (master != -1) // && collapser.unaffectedEdge(edgeI))
- {
- collapser.collapseEdge(edgeI, master);
- nCollapsed++;
- }
+ collapseEdge[edgeI] = 0;
}
}
- return nCollapsed;
-}
-
-
-// Faces which have edges just larger than collapse length but faces which
-// are very small. This one tries to collapse them if it can be done with
-// edge collapse. For faces where a face gets replace by two edges use
-// collapseFaces
-label collapseHighAspectFaces
-(
- const polyMesh& mesh,
- const PackedBoolList& boundaryPoint,
- const scalar areaFac,
- const scalar edgeRatio,
- edgeCollapser& collapser
-)
-{
- const pointField& points = mesh.points();
- const edgeList& edges = mesh.edges();
- const faceList& faces = mesh.faces();
- const labelListList& faceEdges = mesh.faceEdges();
- scalarField magArea(mag(mesh.faceAreas()));
+ // Check whether edge point order is reversed from mesh to coupledPatch
+// const globalMeshData& globalData = mesh.globalData();
+// const mapDistribute& map = globalData.globalEdgeSlavesMap();
+// const labelList& coupledMeshEdges = globalData.coupledPatchMeshEdges();
+// const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
+// const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation();
+// PackedBoolList meshToPatchSameOrientation(coupledMeshEdges.size(), true);
+//
+// forAll(coupledMeshEdges, eI)
+// {
+// const label meshEdgeIndex = coupledMeshEdges[eI];
+//
+// if (collapseEdge[meshEdgeIndex] != -1)
+// {
+// const edge& meshEdge = edges[meshEdgeIndex];
+// const edge& coupledPatchEdge = coupledPatch.edges()[eI];
+//
+// if
+// (
+// meshEdge[0] == coupledPatch.meshPoints()[coupledPatchEdge[1]]
+// && meshEdge[1] == coupledPatch.meshPoints()[coupledPatchEdge[0]]
+// )
+// {
+// meshToPatchSameOrientation[eI] = false;
+// }
+// }
+// }
+//
+//
+// labelList cppEdgeData(coupledMeshEdges.size(), -1);
+//
+// forAll(coupledMeshEdges, eI)
+// {
+// const label meshEdgeIndex = coupledMeshEdges[eI];
+//
+// if (collapseEdge[meshEdgeIndex] != -1)
+// {
+// if (meshToPatchSameOrientation[eI] == cppOrientation[eI])
+// {
+// cppEdgeData[eI] = 0;
+// }
+// else
+// {
+// cppEdgeData[eI] = 1;
+// }
+// }
+// }
+//
+//
+// // Synchronise cppEdgeData
+// // Use minEqOp reduction, so that edge will only be collapsed on processor
+// // boundary if both processors agree to collapse it
+// globalData.syncData
+// (
+// cppEdgeData,
+// globalData.globalEdgeSlaves(),
+// globalData.globalEdgeTransformedSlaves(),
+// map,
+// minEqOp<label>()
+// );
+//
+//
+// forAll(coupledMeshEdges, eI)
+// {
+// const label meshEdgeIndex = coupledMeshEdges[eI];
+//
+// if (collapseEdge[meshEdgeIndex] != -1)
+// {
+// if (meshToPatchSameOrientation[eI] == cppOrientation[eI])
+// {
+// collapseEdge[meshEdgeIndex] = 0;
+// }
+// else
+// {
+// collapseEdge[meshEdgeIndex] = 1;
+// }
+// }
+// }
- label maxIndex = findMax(magArea);
+ label nCollapsed = 0;
- scalar minArea = areaFac * magArea[maxIndex];
+ DynamicList<label> initPoints(mesh.nPoints());
+ DynamicList<pointEdgeCollapse> initPointInfo(mesh.nPoints());
+ allPointInfo.resize(mesh.nPoints());
- Info<< "Max face area:" << magArea[maxIndex] << endl
- << "Collapse area factor:" << areaFac << endl
- << "Collapse area:" << minArea << endl;
+ globalIndex globalStrings(mesh.nPoints());
- label nCollapsed = 0;
+ List<pointEdgeCollapse> allEdgeInfo(mesh.nEdges());
+ forAll(allEdgeInfo, edgeI)
+ {
+ allEdgeInfo[edgeI] = pointEdgeCollapse(vector::zero, -1);
+ }
- forAll(faces, faceI)
+ forAll(edges, edgeI)
{
- if (magArea[faceI] < minArea)
- {
- const face& f = faces[faceI];
+ const edge& e = edges[edgeI];
- // Get the edges in face point order
- labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
+ if (collapseEdge[edgeI] != -1)
+ {
+ const label master =
+ edgeMaster
+ (
+ boundaryPoint,
+ collapseEdge[edgeI],
+ e
+ );
- SortableList<scalar> lengths(fEdges.size());
- forAll(fEdges, i)
+// if (master != -1)
{
- lengths[i] = edges[fEdges[i]].mag(points);
- }
- lengths.sort();
-
-
- label edgeI = -1;
+ pointEdgeCollapse pec
+ (
+ points[master],
+ globalStrings.toGlobal(master)
+ );
- if (f.size() == 4)
- {
- // Compare second largest to smallest
- if (lengths[2] > edgeRatio*lengths[0])
- {
- // Collapse smallest only. Triangle should be cleared
- // next time around.
- edgeI = fEdges[lengths.indices()[0]];
- }
- }
- else if (f.size() == 3)
- {
- // Compare second largest to smallest
- if (lengths[1] > edgeRatio*lengths[0])
- {
- edgeI = fEdges[lengths.indices()[0]];
- }
- }
+ allEdgeInfo[edgeI] = pec;
+ initPointInfo.append(pec);
+ initPoints.append(e.start());
- if (edgeI != -1)
- {
- label master = edgeMaster(boundaryPoint, edges[edgeI]);
+ initPointInfo.append(pec);
+ initPoints.append(e.end());
- if (master != -1)// && collapser.unaffectedEdge(edgeI))
- {
- collapser.collapseEdge(edgeI, master);
- nCollapsed++;
- }
+ nCollapsed++;
}
}
}
+ PointEdgeWave<pointEdgeCollapse> collapsePropagator
+ (
+ mesh,
+ initPoints,
+ initPointInfo,
+ allPointInfo,
+ allEdgeInfo,
+ mesh.globalData().nTotalPoints() // Maximum number of iterations
+ );
+
return nCollapsed;
}
+// Faces which have edges just larger than collapse length but faces which
+// are very small. This one tries to collapse them if it can be done with
+// edge collapse. For faces where a face gets replace by two edges use
+// collapseFaces
+//label collapseHighAspectFaces
+//(
+// const polyMesh& mesh,
+// const PackedBoolList& boundaryPoint,
+// const Map<label>& processorPoints,
+// const scalar areaFac,
+// const scalar edgeRatio,
+// edgeCollapser& collapser
+//)
+//{
+// const pointField& points = mesh.points();
+// const edgeList& edges = mesh.edges();
+// const faceList& faces = mesh.faces();
+// const labelListList& faceEdges = mesh.faceEdges();
+//
+// scalarField magArea(mag(mesh.faceAreas()));
+//
+// label maxIndex = findMax(magArea);
+//
+// scalar minArea = areaFac * magArea[maxIndex];
+//
+// Info<< "Max face area:" << magArea[maxIndex] << endl
+// << "Collapse area factor:" << areaFac << endl
+// << "Collapse area:" << minArea << endl;
+//
+// label nCollapsed = 0;
+//
+// forAll(faces, faceI)
+// {
+// if (magArea[faceI] < minArea)
+// {
+// const face& f = faces[faceI];
+//
+// // Get the edges in face point order
+// labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
+//
+// SortableList<scalar> lengths(fEdges.size());
+// forAll(fEdges, i)
+// {
+// lengths[i] = edges[fEdges[i]].mag(points);
+// }
+// lengths.sort();
+//
+//
+// label edgeI = -1;
+//
+// if (f.size() == 4)
+// {
+// // Compare second largest to smallest
+// if (lengths[2] > edgeRatio*lengths[0])
+// {
+// // Collapse smallest only. Triangle should be cleared
+// // next time around.
+// edgeI = fEdges[lengths.indices()[0]];
+// }
+// }
+// else if (f.size() == 3)
+// {
+// // Compare second largest to smallest
+// if (lengths[1] > edgeRatio*lengths[0])
+// {
+// edgeI = fEdges[lengths.indices()[0]];
+// }
+// }
+//
+//
+// if (edgeI != -1)
+// {
+// label master = edgeMaster
+// (
+// boundaryPoint,
+// processorPoints,
+// false,
+// edges[edgeI]
+// );
+//
+// if (master != -1)// && collapser.unaffectedEdge(edgeI))
+// {
+// collapser.collapseEdge(edgeI, master);
+// nCollapsed++;
+// }
+// }
+// }
+// }
+//
+// return nCollapsed;
+//}
+
+
void set(const labelList& elems, const bool val, boolList& status)
{
forAll(elems, i)
@@ -340,122 +946,254 @@ void set(const labelList& elems, const bool val, boolList& status)
// Tries to simplify polygons to face of minSize (4=quad, 3=triangle)
-label simplifyFaces
-(
- const polyMesh& mesh,
- const PackedBoolList& boundaryPoint,
- const label minSize,
- const scalar lenGap,
- edgeCollapser& collapser
-)
+//label simplifyFaces
+//(
+// const polyMesh& mesh,
+// const PackedBoolList& boundaryPoint,
+// const Map<label>& processorPoints,
+// const label minSize,
+// const scalar lenGap,
+// edgeCollapser& collapser
+//)
+//{
+// const pointField& points = mesh.points();
+// const edgeList& edges = mesh.edges();
+// const faceList& faces = mesh.faces();
+// const cellList& cells = mesh.cells();
+// const labelListList& faceEdges = mesh.faceEdges();
+// const labelList& faceOwner = mesh.faceOwner();
+// const labelList& faceNeighbour = mesh.faceNeighbour();
+// const labelListList& pointCells = mesh.pointCells();
+// const labelListList& cellEdges = mesh.cellEdges();
+//
+// label nCollapsed = 0;
+//
+// boolList protectedEdge(mesh.nEdges(), false);
+//
+// forAll(faces, faceI)
+// {
+// const face& f = faces[faceI];
+//
+// if
+// (
+// f.size() > minSize
+// && cells[faceOwner[faceI]].size() >= 6
+// && (
+// mesh.isInternalFace(faceI)
+// && cells[faceNeighbour[faceI]].size() >= 6
+// )
+// )
+// {
+// // Get the edges in face point order
+// labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
+//
+// SortableList<scalar> lengths(fEdges.size());
+// forAll(fEdges, i)
+// {
+// lengths[i] = edges[fEdges[i]].mag(points);
+// }
+// lengths.sort();
+//
+//
+// // Now find a gap in length between consecutive elements greater
+// // than lenGap.
+//
+// label gapPos = -1;
+//
+// for (label i = f.size()-1-minSize; i >= 0; --i)
+// {
+// if (lengths[i+1] > lenGap*lengths[i])
+// {
+// gapPos = i;
+//
+// break;
+// }
+// }
+//
+// if (gapPos != -1)
+// {
+// //for (label i = gapPos; i >= 0; --i)
+// label i = 0; // Hack: collapse smallest edge only.
+// {
+// label edgeI = fEdges[lengths.indices()[i]];
+//
+// if (!protectedEdge[edgeI])
+// {
+// const edge& e = edges[edgeI];
+//
+// label master
+// = edgeMaster(boundaryPoint, processorPoints, false, e);
+//
+// if (master != -1)
+// {
+// collapser.collapseEdge(edgeI, master);
+//
+// // Protect all other edges on all cells using edge
+// // points.
+//
+// const labelList& pCells0 = pointCells[e[0]];
+//
+// forAll(pCells0, i)
+// {
+// set(cellEdges[pCells0[i]], true, protectedEdge);
+// }
+// const labelList& pCells1 = pointCells[e[1]];
+//
+// forAll(pCells1, i)
+// {
+// set(cellEdges[pCells1[i]], true, protectedEdge);
+// }
+//
+// nCollapsed++;
+// }
+// }
+// }
+// }
+// }
+// }
+//
+// return nCollapsed;
+//}
+
+
+labelHashSet checkMeshQuality(const polyMesh& mesh)
{
- const pointField& points = mesh.points();
- const edgeList& edges = mesh.edges();
- const faceList& faces = mesh.faces();
- const cellList& cells = mesh.cells();
- const labelListList& faceEdges = mesh.faceEdges();
- const labelList& faceOwner = mesh.faceOwner();
- const labelList& faceNeighbour = mesh.faceNeighbour();
- const labelListList& pointCells = mesh.pointCells();
- const labelListList& cellEdges = mesh.cellEdges();
+ //mesh.checkMesh(true);
+ labelHashSet freezePoints;
- label nCollapsed = 0;
+ IOdictionary meshQualityDict
+ (
+ IOobject
+ (
+ "meshQualityControls",
+ mesh.time().system(),
+ mesh.time(),
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+ );
- boolList protectedEdge(mesh.nEdges(), false);
+ labelHashSet badFaces(mesh.nFaces()/100);
+ DynamicList<label> checkFaces(mesh.nFaces());
- forAll(faces, faceI)
- {
- const face& f = faces[faceI];
+ const vectorField& fAreas = mesh.faceAreas();
- if
- (
- f.size() > minSize
- && cells[faceOwner[faceI]].size() >= 6
- && (
- mesh.isInternalFace(faceI)
- && cells[faceNeighbour[faceI]].size() >= 6
- )
- )
+ scalar faceAreaLimit = SMALL;
+
+ forAll(fAreas, fI)
+ {
+ if (mag(fAreas[fI]) > faceAreaLimit)
{
- // Get the edges in face point order
- labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
+ checkFaces.append(fI);
+ }
+ }
- SortableList<scalar> lengths(fEdges.size());
- forAll(fEdges, i)
- {
- lengths[i] = edges[fEdges[i]].mag(points);
- }
- lengths.sort();
+ motionSmoother::checkMesh
+ (
+ false,
+ mesh,
+ meshQualityDict,
+ checkFaces,
+ badFaces
+ );
+ label nBadFaces = badFaces.size();
+ reduce(nBadFaces, sumOp<label>());
- // Now find a gap in length between consecutive elements greater
- // than lenGap.
+ Info<< nl << "Number of bad faces : " << nBadFaces << endl;
- label gapPos = -1;
+ forAllConstIter(labelHashSet, badFaces, iter)
+ {
+ const face& f = mesh.faces()[iter.key()];
- for (label i = f.size()-1-minSize; i >= 0; --i)
- {
- if (lengths[i+1] > lenGap*lengths[i])
- {
- gapPos = i;
+ forAll(f, pI)
+ {
+ freezePoints.insert(f[pI]);
+ }
+ }
- break;
- }
- }
+// const edgeList& edges = mesh.edges();
+//
+// label nFrozenEdges = 0;
+//
+// OFstream str("frozenEdges.obj");
+//
+// label count = 0;
+// forAll(edges, eI)
+// {
+// const edge& e = edges[eI];
+//
+// if (freezePoints.found(e[0]) && freezePoints.found(e[1]))
+// {
+// freezeEdges[eI] = true;
+// nFrozenEdges++;
+// }
+// }
+
+ return freezePoints;
+}
- if (gapPos != -1)
- {
- //for (label i = gapPos; i >= 0; --i)
- label i = 0; // Hack: collapse smallest edge only.
- {
- label edgeI = fEdges[lengths.indices()[i]];
- if (!protectedEdge[edgeI])
- {
- const edge& e = edges[edgeI];
+// Mark boundary points
+// boundaryPoint:
+// + -1 : point not on boundary
+// + 0 : point on a real boundary
+// + >0 : point on a processor patch with that ID
+labelList findBoundaryPoints(const polyMesh& mesh)
+{
+ const faceList& faces = mesh.faces();
+ const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
- label master = edgeMaster(boundaryPoint, e);
- if (master != -1)
- {
- collapser.collapseEdge(edgeI, master);
+ labelList boundaryPoint(mesh.nPoints(), -1);
- // Protect all other edges on all cells using edge
- // points.
+ // Get all points on a boundary
+ label nIntFaces = mesh.nInternalFaces();
+ for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
+ {
+ const face& f = faces[faceI];
- const labelList& pCells0 = pointCells[e[0]];
+ forAll(f, fp)
+ {
+ boundaryPoint[f[fp]] = 0;
+ }
+ }
- forAll(pCells0, i)
- {
- set(cellEdges[pCells0[i]], true, protectedEdge);
- }
- const labelList& pCells1 = pointCells[e[1]];
+ // Get all processor boundary points and the processor patch label
+ // that they are on.
+ forAll(bMesh, patchI)
+ {
+ const polyPatch& patch = bMesh[patchI];
- forAll(pCells1, i)
- {
- set(cellEdges[pCells1[i]], true, protectedEdge);
- }
+ if (isA<processorPolyPatch>(patch))
+ {
+ const processorPolyPatch& pPatch =
+ refCast<const processorPolyPatch>(patch);
- nCollapsed++;
- }
- }
+ forAll(pPatch, fI)
+ {
+ const face& f = pPatch[fI];
+
+ forAll(f, fp)
+ {
+ boundaryPoint[f[fp]] = patchI;
}
}
}
}
- return nCollapsed;
+ return boundaryPoint;
}
// Main program:
-
int main(int argc, char *argv[])
{
# include "addOverwriteOption.H"
- argList::noParallel();
+
argList::validArgs.append("edge length [m]");
argList::validArgs.append("merge angle (degrees)");
+ argList::addOption("minLenFactor", "scalar", "edge length factor");
# include "setRootCase.H"
# include "createTime.H"
@@ -463,8 +1201,11 @@ int main(int argc, char *argv[])
# include "createPolyMesh.H"
const word oldInstance = mesh.pointsInstance();
- const scalar minLen = args.argRead<scalar>(1);
+ scalar minLen = args.argRead<scalar>(1);
const scalar angle = args.argRead<scalar>(2);
+ const scalar minLenFactor
+ = args.optionLookupOrDefault<scalar>("minLenFactor", 0.5);
+
const bool overwrite = args.optionFound("overwrite");
scalar maxCos = Foam::cos(degToRad(angle));
@@ -475,108 +1216,196 @@ int main(int argc, char *argv[])
<< " degrees" << nl
<< endl;
+ Info<< "If an invalid mesh is generated then the edge length will be " << nl
+ << "multiplied by a factor of " << minLenFactor << " and collapsing "
+ << "will be reattempted" << nl << endl;
bool meshChanged = false;
- while (true)
- {
- const faceList& faces = mesh.faces();
+ checkMeshQuality(mesh);
- // Get all points on the boundary
- PackedBoolList boundaryPoint(mesh.nPoints());
+ autoPtr<fvMesh> fvMeshPtr;
- label nIntFaces = mesh.nInternalFaces();
- for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
- {
- const face& f = faces[faceI];
+ scalarList freezeEdges(mesh.nEdges(), 1.0);
- forAll(f, fp)
- {
- boundaryPoint.set(f[fp], 1);
- }
- }
+ do
+ {
+ label nIterations = 0;
+ label nFrozenEdges = 0;
- // Edge collapsing engine
- edgeCollapser collapser(mesh);
+ fvMeshPtr.reset(new fvMesh(mesh));
+ fvMesh& fvMeshRef = fvMeshPtr();
+ // Contains new point label for original points
+ labelList pointMap(identity(mesh.nPoints()));
- // Collapse all edges that are too small.
- label nCollapsed =
- collapseSmallEdges
- (
- mesh,
- boundaryPoint,
- minLen,
- collapser
- );
- Info<< "Collapsing " << nCollapsed << " small edges" << endl;
+ scalarList tmpFreezeEdges = freezeEdges;
+ autoPtr<mapPolyMesh> morphMap;
- // Remove midpoints on straight edges.
- if (nCollapsed == 0)
+ while (true)
{
- nCollapsed = mergeEdges(mesh, maxCos, collapser);
- Info<< "Collapsing " << nCollapsed << " in line edges" << endl;
- }
+ Info<< "Iteration " << nIterations << incrIndent << endl;
+ labelList boundaryPoint = findBoundaryPoints(fvMeshRef);
- // Remove small sliver faces that can be collapsed to single edge
- if (nCollapsed == 0)
- {
- nCollapsed =
- collapseHighAspectFaces
+ List<pointEdgeCollapse> allPointInfo;
+
+ // Collapse all edges that are too small.
+ label nSmallCollapsed =
+ collapseSmallEdges
(
- mesh,
+ fvMeshRef,
+ tmpFreezeEdges,
boundaryPoint,
- 1e-9, // factor of largest face area
- 5, // factor between smallest and largest edge on
- // face
- collapser
+ minLen,
+ allPointInfo
);
- Info<< "Collapsing " << nCollapsed
+
+ reduce(nSmallCollapsed, sumOp<label>());
+
+ Info<< indent << "Collapsing " << nSmallCollapsed
+ << " small edges" << endl;
+
+
+
+
+
+ label nMerged = 0;
+
+ // Remove midpoints on straight edges.
+ if (nSmallCollapsed == 0)
+ {
+ //nMerged = mergeEdges(fvMeshRef, maxCos, allPointInfo);
+ }
+
+ reduce(nMerged, sumOp<label>());
+
+ Info<< indent << "Collapsing " << nMerged << " in line edges"
+ << endl;
+
+
+
+
+
+
+ label nSliversCollapsed = 0;
+
+ // Remove small sliver faces that can be collapsed to single edge
+// if (nSmallCollapsed == 0 && nMerged == 0)
+// {
+// nSliversCollapsed =
+// collapseHighAspectFaces
+// (
+// mesh,
+// boundaryPoint,
+// processorPoints,
+// 1E-9,
+// 5,
+// collapser
+// );
+// }
+
+ reduce(nSliversCollapsed, sumOp<label>());
+
+ Info<< indent << "Collapsing " << nSliversCollapsed
<< " small high aspect ratio faces" << endl;
- }
- // Simplify faces to quads wherever possible
- //if (nCollapsed == 0)
- //{
- // nCollapsed =
- // simplifyFaces
- // (
- // mesh,
- // boundaryPoint,
- // 4, // minimum size of face
- // 0.2, // gap in edge lengths on face
- // collapser
- // );
- // Info<< "Collapsing " << nCollapsed << " polygonal faces" << endl;
- //}
+ // Simplify faces to quads wherever possible
+ //if (nCollapsed == 0)
+ //{
+ // nCollapsed =
+ // simplifyFaces
+ // (
+ // mesh,
+ // boundaryPoint,
+ // 4, // minimum size of face
+ // 0.2, // gap in edge lengths on face
+ // collapser
+ // );
+ // Info<< "Collapsing " << nCollapsed << " polygonal faces"
+ // << endl;
+ //}
+
+
+ label totalCollapsed =
+ nSmallCollapsed
+ + nMerged
+ + nSliversCollapsed;
+
+ polyTopoChange meshMod(fvMeshRef);
+
+ // Insert mesh refinement into polyTopoChange.
+ setRefinement(fvMeshRef, meshMod, allPointInfo);
+
+ // Do all changes
+ Info<< indent << "Applying changes to the mesh" << nl
+ << decrIndent << endl;
+
+ morphMap = meshMod.changeMesh(fvMeshRef, false);
+
+// // Contains new point label for old points
+// const labelList& reversePointMap = morphMap().reversePointMap();
+//
+// forAll(pointMap, pI)
+// {
+// const label originalPoint = pI;
+// const label currentPoint = pointMap[pI];
+//
+// if (currentPoint < reversePointMap.size())
+// {
+// const label newPoint = reversePointMap[currentPoint];
+//
+// if (newPoint != -1)
+// {
+// pointMap[originalPoint] = newPoint;
+// }
+// }
+// }
+
+ if (totalCollapsed == 0)
+ {
+ labelHashSet freezePoints = checkMeshQuality(fvMeshRef);
- if (nCollapsed == 0)
- {
- break;
- }
+ label nFreezePoints = freezePoints.size();
+ reduce(nFreezePoints, sumOp<label>());
- polyTopoChange meshMod(mesh);
+ nFrozenEdges = nFreezePoints;
- // Insert mesh refinement into polyTopoChange.
- collapser.setRefinement(meshMod);
+ Info<< "Number of frozen points : " << nFreezePoints
+ << endl;
- // Do all changes
- Info<< "Morphing ..." << endl;
+ break;
+ }
- autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
+ if (morphMap().hasMotionPoints())
+ {
+ fvMeshRef.movePoints(morphMap().preMotionPoints());
+ }
- collapser.updateMesh(morphMap());
+ meshChanged = true;
- if (morphMap().hasMotionPoints())
+ nIterations++;
+ }
+
+ if (nFrozenEdges > 0)
{
- mesh.movePoints(morphMap().preMotionPoints());
+ minLen *= minLenFactor;
}
- meshChanged = true;
- }
+ reduce(nFrozenEdges, sumOp<label>());
+
+ Info<< "Number of frozen edges : " << nFrozenEdges << nl
+ << endl;
+
+ if (nFrozenEdges == 0)
+ {
+ break;
+ }
+
+ } while (true);
+
if (meshChanged)
{
@@ -587,14 +1416,21 @@ int main(int argc, char *argv[])
}
else
{
- mesh.setInstance(oldInstance);
+ fvMeshPtr().setInstance(oldInstance);
}
- Info<< "Writing collapsed mesh to time " << runTime.timeName() << endl;
+ Info<< nl << "Writing collapsed mesh to time "
+ << runTime.timeName() << nl << endl;
- mesh.write();
+ fvMeshPtr().write();
}
+ Info<< "Final minimum length : " << minLen << " m" << nl << endl;
+
+ Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+
Info<< "End\n" << endl;
return 0;
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
new file mode 100644
index 0000000000000000000000000000000000000000..6830e1c8f56ed18c006241ffb16388012da4b042
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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 "pointEdgeCollapse.H"
+
+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+ Foam::Ostream& os,
+ const Foam::pointEdgeCollapse& wDist
+)
+{
+ return os
+ << wDist.collapsePoint_ << wDist.collapseIndex_;
+}
+
+Foam::Istream& Foam::operator>>
+(
+ Foam::Istream& is,
+ Foam::pointEdgeCollapse& wDist
+)
+{
+ return is
+ >> wDist.collapsePoint_ >> wDist.collapseIndex_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
new file mode 100644
index 0000000000000000000000000000000000000000..cd8383b8b33b6ed6ddb0b03a172aa632205f3cee
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+ Foam::pointEdgeCollapse
+
+Description
+ Determines length of string of edges walked to point.
+
+SourceFiles
+ pointEdgeCollapseI.H
+ pointEdgeCollapse.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointEdgeCollapse_H
+#define pointEdgeCollapse_H
+
+#include "point.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+ Class pointEdgeCollapse Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointEdgeCollapse
+{
+ // Private data
+
+ //- Collapse location
+ point collapsePoint_;
+
+ //- Collapse string index
+ label collapseIndex_;
+
+
+ // Private Member Functions
+
+ //- Evaluate distance to point.
+ template<class TrackingData>
+ inline bool update
+ (
+ const pointEdgeCollapse& w2,
+ const scalar tol,
+ TrackingData& td
+ );
+
+
+ //- Check for same coordinate
+ inline bool samePoint(const point& pt) const;
+
+public:
+
+ // Constructors
+
+ //- Construct null
+ inline pointEdgeCollapse();
+
+ //- Construct from components
+ inline pointEdgeCollapse
+ (
+ const point& collapsePoint,
+ const label collapseIndex
+ );
+
+
+ // Member Functions
+
+ // Access
+
+ inline const point& collapsePoint() const;
+
+ inline label collapseIndex() const;
+
+
+ // Needed by meshWave
+
+ //- Check whether origin has been changed at all or
+ // still contains original (invalid) value.
+ template<class TrackingData>
+ inline bool valid(TrackingData& td) const;
+
+ //- Convert origin to relative vector to leaving point
+ // (= point coordinate)
+ template<class TrackingData>
+ inline void leaveDomain
+ (
+ const polyPatch& patch,
+ const label patchPointI,
+ const point& pos,
+ TrackingData& td
+ );
+
+ //- Convert relative origin to absolute by adding entering point
+ template<class TrackingData>
+ inline void enterDomain
+ (
+ const polyPatch& patch,
+ const label patchPointI,
+ const point& pos,
+ TrackingData& td
+ );
+
+ //- Apply rotation matrix to origin
+ template<class TrackingData>
+ inline void transform
+ (
+ const tensor& rotTensor,
+ TrackingData& td
+ );
+
+ //- Influence of edge on point
+ template<class TrackingData>
+ inline bool updatePoint
+ (
+ const polyMesh& mesh,
+ const label pointI,
+ const label edgeI,
+ const pointEdgeCollapse& edgeInfo,
+ const scalar tol,
+ TrackingData& td
+ );
+
+ //- Influence of different value on same point.
+ // Merge new and old info.
+ template<class TrackingData>
+ inline bool updatePoint
+ (
+ const polyMesh& mesh,
+ const label pointI,
+ const pointEdgeCollapse& newPointInfo,
+ const scalar tol,
+ TrackingData& td
+ );
+
+ //- Influence of different value on same point.
+ // No information about current position whatsoever.
+ template<class TrackingData>
+ inline bool updatePoint
+ (
+ const pointEdgeCollapse& newPointInfo,
+ const scalar tol,
+ TrackingData& td
+ );
+
+ //- Influence of point on edge.
+ template<class TrackingData>
+ inline bool updateEdge
+ (
+ const polyMesh& mesh,
+ const label edgeI,
+ const label pointI,
+ const pointEdgeCollapse& pointInfo,
+ const scalar tol,
+ TrackingData& td
+ );
+
+ //- Same (like operator==)
+ template<class TrackingData>
+ inline bool equal(const pointEdgeCollapse&, TrackingData&)
+ const;
+
+
+ // Member Operators
+
+ //Note: Used to determine whether to call update.
+ inline bool operator==(const pointEdgeCollapse&) const;
+ inline bool operator!=(const pointEdgeCollapse&) const;
+
+
+ // IOstream Operators
+
+ friend Ostream& operator<<(Ostream&, const pointEdgeCollapse&);
+ friend Istream& operator>>(Istream&, pointEdgeCollapse&);
+};
+
+
+//- Data associated with pointEdgeCollapse type are contiguous
+template<>
+inline bool contiguous<pointEdgeCollapse>()
+{
+ return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointEdgeCollapseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
new file mode 100644
index 0000000000000000000000000000000000000000..3fb1468a88f799f92c872f0ecc9fe0b60bf23747
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
@@ -0,0 +1,313 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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 "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+// Update this with w2.
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::update
+(
+ const pointEdgeCollapse& w2,
+ const scalar tol,
+ TrackingData& td
+)
+{
+ if (!w2.valid(td))
+ {
+ FatalErrorIn("pointEdgeCollapse::update(..)")
+ << "problem." << abort(FatalError);
+ }
+
+ if (w2.collapseIndex_ == -1)
+ {
+ // Not marked for collapse; only happens on edges.
+ return false;
+ }
+
+ if (!valid(td))
+ {
+ operator=(w2);
+ return true;
+ }
+ else
+ {
+ // Same coordinate. Same string?
+ if (w2.collapseIndex_ < collapseIndex_)
+ {
+ // Take over string index from w2 (and also coordinate but this
+ // was same)
+ operator=(w2);
+ return true;
+ }
+ else if (w2.collapseIndex_ == collapseIndex_)
+ {
+ bool identicalPoint = samePoint(w2.collapsePoint_);
+ bool nearer = magSqr(w2.collapsePoint_) < magSqr(collapsePoint_);
+ if (nearer)
+ {
+ operator=(w2);
+ }
+ if (identicalPoint)
+ {
+ return false;
+ }
+ else
+ {
+ return nearer;
+ }
+ }
+ else
+ {
+ return false;
+ }
+
+// if (samePoint(w2.collapsePoint_))
+// {
+// // Same coordinate. Same string?
+// if (w2.collapseIndex_ < collapseIndex_)
+// {
+// // Take over string index from w2 (and also coordinate but this
+// // was same)
+// operator=(w2);
+// return true;
+// }
+// else
+// {
+// return false;
+// }
+// }
+// else
+// {
+// // Find nearest coordinate
+// if (magSqr(w2.collapsePoint_) < magSqr(collapsePoint_))
+// {
+// operator=(w2);
+// return true;
+// }
+// else
+// {
+// return false;
+// }
+// }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointEdgeCollapse::pointEdgeCollapse()
+:
+ collapsePoint_(GREAT, GREAT, GREAT),
+ collapseIndex_(-2)
+{}
+
+
+// Construct from origin, distance
+inline Foam::pointEdgeCollapse::pointEdgeCollapse
+(
+ const point& collapsePoint,
+ const label collapseIndex
+)
+:
+ collapsePoint_(collapsePoint),
+ collapseIndex_(collapseIndex)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+inline const Foam::point& Foam::pointEdgeCollapse::collapsePoint() const
+{
+ return collapsePoint_;
+}
+
+
+inline Foam::label Foam::pointEdgeCollapse::collapseIndex() const
+{
+ return collapseIndex_;
+}
+
+
+inline bool Foam::pointEdgeCollapse::samePoint(const point& pt) const
+{
+ bool isLegal1 = (cmptMin(collapsePoint_) < 0.5*GREAT);
+ bool isLegal2 = (cmptMin(pt) < 0.5*GREAT);
+
+ if (isLegal1 && isLegal2)
+ {
+ return mag(collapsePoint_ - pt) < 1e-9;//SMALL;
+ }
+ else
+ {
+ return isLegal1 == isLegal2;
+ }
+}
+
+
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::valid(TrackingData& td) const
+{
+ return collapseIndex_ != -2;
+}
+
+
+template<class TrackingData>
+inline void Foam::pointEdgeCollapse::leaveDomain
+(
+ const polyPatch& patch,
+ const label patchPointI,
+ const point& coord,
+ TrackingData& td
+)
+{
+ collapsePoint_ -= coord;
+}
+
+
+template<class TrackingData>
+inline void Foam::pointEdgeCollapse::transform
+(
+ const tensor& rotTensor,
+ TrackingData& td
+)
+{
+ collapsePoint_ = Foam::transform(rotTensor, collapsePoint_);
+}
+
+
+// Update absolute geometric quantities. Note that distance (dist_)
+// is not affected by leaving/entering domain.
+template<class TrackingData>
+inline void Foam::pointEdgeCollapse::enterDomain
+(
+ const polyPatch& patch,
+ const label patchPointI,
+ const point& coord,
+ TrackingData& td
+)
+{
+ // back to absolute form
+ collapsePoint_ += coord;
+}
+
+
+// Update this with information from connected edge
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updatePoint
+(
+ const polyMesh& mesh,
+ const label pointI,
+ const label edgeI,
+ const pointEdgeCollapse& edgeInfo,
+ const scalar tol,
+ TrackingData& td
+)
+{
+ return update(edgeInfo, tol, td);
+}
+
+
+// Update this with new information on same point
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updatePoint
+(
+ const polyMesh& mesh,
+ const label pointI,
+ const pointEdgeCollapse& newPointInfo,
+ const scalar tol,
+ TrackingData& td
+)
+{
+ return update(newPointInfo, tol, td);
+}
+
+
+// Update this with new information on same point. No extra information.
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updatePoint
+(
+ const pointEdgeCollapse& newPointInfo,
+ const scalar tol,
+ TrackingData& td
+)
+{
+ return update(newPointInfo, tol, td);
+}
+
+
+// Update this with information from connected point
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updateEdge
+(
+ const polyMesh& mesh,
+ const label edgeI,
+ const label pointI,
+ const pointEdgeCollapse& pointInfo,
+ const scalar tol,
+ TrackingData& td
+)
+{
+ return update(pointInfo, tol, td);
+}
+
+
+template <class TrackingData>
+inline bool Foam::pointEdgeCollapse::equal
+(
+ const pointEdgeCollapse& rhs,
+ TrackingData& td
+) const
+{
+ return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeCollapse::operator==
+(
+ const Foam::pointEdgeCollapse& rhs
+) const
+{
+ return
+ collapseIndex_ == rhs.collapseIndex_
+ && samePoint(rhs.collapsePoint_);
+}
+
+
+inline bool Foam::pointEdgeCollapse::operator!=
+(
+ const Foam::pointEdgeCollapse& rhs
+) const
+{
+ return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index 415c6facdca40204dfe3e097e10452ec756a9933..22872abf47d6b1c5f87150831ef4ce1d22a82b3a 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -27,6 +27,10 @@ License
#include "polyMesh.H"
#include "polyTopoChange.H"
#include "ListOps.H"
+#include "globalMeshData.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "syncTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -121,6 +125,7 @@ void Foam::edgeCollapser::filterFace(const label faceI, face& f) const
}
}
+
// Check for pinched face. Tries to correct
// - consecutive duplicate vertex. Removes duplicate vertex.
// - duplicate vertex with one other vertex in between (spike).
@@ -190,15 +195,15 @@ void Foam::edgeCollapser::printRegions() const
if (master != -1)
{
- Info<< "Region:" << regionI << nl
+ Pout<< "Region:" << regionI << nl
<< " master:" << master
- << ' ' << mesh_.points()[master] << nl;
+ << ' ' << pointRegionMasterLocation_[regionI] << nl;
forAll(pointRegion_, pointI)
{
if (pointRegion_[pointI] == regionI && pointI != master)
{
- Info<< " slave:" << pointI
+ Pout<< " slave:" << pointI
<< ' ' << mesh_.points()[pointI] << nl;
}
}
@@ -272,6 +277,7 @@ Foam::edgeCollapser::edgeCollapser(const polyMesh& mesh)
:
mesh_(mesh),
pointRegion_(mesh.nPoints(), -1),
+ pointRegionMasterLocation_(mesh.nPoints() / 100),
pointRegionMaster_(mesh.nPoints() / 100),
freeRegions_()
{}
@@ -289,6 +295,8 @@ bool Foam::edgeCollapser::unaffectedEdge(const label edgeI) const
bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
{
+ const pointField& points = mesh_.points();
+
const edge& e = mesh_.edges()[edgeI];
label pointRegion0 = pointRegion_[e[0]];
@@ -310,7 +318,7 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
{
FatalErrorIn
("edgeCollapser::collapseEdge(const label, const label)")
- << "Problem : freeed region :" << freeRegion
+ << "Problem : freed region :" << freeRegion
<< " has already master "
<< pointRegionMaster_[freeRegion]
<< abort(FatalError);
@@ -327,13 +335,22 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
pointRegion_[e[1]] = freeRegion;
pointRegionMaster_(freeRegion) = master;
+ pointRegionMasterLocation_(freeRegion) = points[master];
}
else
{
// e[1] is part of collapse network, e[0] not. Add e0 to e1 region.
pointRegion_[e[0]] = pointRegion1;
- pointRegionMaster_[pointRegion1] = master;
+ if
+ (
+ pointRegionMaster_[pointRegion1] == e[0]
+ || pointRegionMaster_[pointRegion1] == e[1]
+ )
+ {
+ pointRegionMaster_[pointRegion1] = master;
+ pointRegionMasterLocation_[pointRegion1] = points[master];
+ }
}
}
else
@@ -343,7 +360,15 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
// e[0] is part of collapse network. Add e1 to e0 region
pointRegion_[e[1]] = pointRegion0;
- pointRegionMaster_[pointRegion0] = master;
+ if
+ (
+ pointRegionMaster_[pointRegion0] == e[0]
+ || pointRegionMaster_[pointRegion0] == e[1]
+ )
+ {
+ pointRegionMaster_[pointRegion0] = master;
+ pointRegionMasterLocation_[pointRegion0] = points[master];
+ }
}
else if (pointRegion0 != pointRegion1)
{
@@ -356,6 +381,9 @@ bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
// Use minRegion as region for combined net, free maxRegion.
pointRegionMaster_[minRegion] = master;
pointRegionMaster_[maxRegion] = -1;
+ pointRegionMasterLocation_[minRegion] = points[master];
+ pointRegionMasterLocation_[maxRegion] = point(0, 0, 0);
+
freeRegions_.insert(maxRegion);
if (minRegion != pointRegion0)
@@ -380,12 +408,61 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
const labelList& faceNeighbour = mesh_.faceNeighbour();
const labelListList& pointFaces = mesh_.pointFaces();
const labelListList& cellEdges = mesh_.cellEdges();
+ const pointZoneMesh& pointZones = mesh_.pointZones();
+
- // Print regions:
- //printRegions()
bool meshChanged = false;
+ // Synchronise pointRegionMasterLocation_
+ const globalMeshData& globalData = mesh_.globalData();
+ const mapDistribute& map = globalData.globalPointSlavesMap();
+ const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
+ const labelList& meshPoints = coupledPatch.meshPoints();
+ const Map<label>& meshPointMap = coupledPatch.meshPointMap();
+
+
+ List<point> newPoints = coupledPatch.localPoints();
+
+ for (label pI = 0; pI < coupledPatch.nPoints(); ++pI)
+ {
+ const label pointRegionMaster = pointRegion_[meshPoints[pI]];
+
+ if (pointRegionMaster != -1)
+ {
+ newPoints[pI]
+ = pointRegionMasterLocation_[pointRegionMaster];
+ }
+ }
+
+ globalData.syncData
+ (
+ newPoints,
+ globalData.globalPointSlaves(),
+ globalData.globalPointTransformedSlaves(),
+ map,
+ minMagSqrEqOp<point>()
+ );
+
+ OFstream str1("newPoints_" + name(Pstream::myProcNo()) + ".obj");
+ forAll(pointRegion_, pI)
+ {
+ if (meshPointMap.found(pI))
+ {
+ meshTools::writeOBJ(str1, newPoints[meshPointMap[pI]]);
+ }
+ }
+
+ for (label pI = 0; pI < coupledPatch.nPoints(); ++pI)
+ {
+ const label pointRegionMaster = pointRegion_[meshPoints[pI]];
+
+ if (pointRegionMaster != -1)
+ {
+ pointRegionMasterLocation_[pointRegionMaster]
+ = newPoints[pI];
+ }
+ }
// Current faces (is also collapseStatus: f.size() < 3)
faceList newFaces(mesh_.faces());
@@ -393,7 +470,6 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
// Current cellCollapse status
boolList cellRemoved(mesh_.nCells(), false);
-
do
{
// Update face collapse from edge collapses
@@ -521,12 +597,49 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
}
}
+ // Modify the point location of the remaining points
+ forAll(pointRegion_, pointI)
+ {
+ const label pointRegion = pointRegion_[pointI];
+
+ if
+ (
+ !pointRemoved(pointI)
+ && meshPointMap.found(pointI)
+ )
+ {
+ meshMod.modifyPoint
+ (
+ pointI,
+ newPoints[meshPointMap[pointI]],
+ pointZones.whichZone(pointI),
+ false
+ );
+ }
+ else if
+ (
+ pointRegion != -1
+ && !pointRemoved(pointI)
+ && !meshPointMap.found(pointI)
+ )
+ {
+ const point& collapsePoint
+ = pointRegionMasterLocation_[pointRegion];
+
+ meshMod.modifyPoint
+ (
+ pointI,
+ collapsePoint,
+ pointZones.whichZone(pointI),
+ false
+ );
+ }
+ }
const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
const faceZoneMesh& faceZones = mesh_.faceZones();
-
// Renumber faces that use points
forAll(pointRegion_, pointI)
{
@@ -585,6 +698,9 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
}
}
+ // Print regions:
+// printRegions();
+
return meshChanged;
}
@@ -593,8 +709,10 @@ void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map)
{
pointRegion_.setSize(mesh_.nPoints());
pointRegion_ = -1;
+
// Reset count, do not remove underlying storage
pointRegionMaster_.clear();
+ pointRegionMasterLocation_.clear();
freeRegions_.clear();
}
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
index 98d8ef30d5bce148647499c65cbcbde299df9fc1..3a10a9319749081af5960967cb6eaf09ec3e17d4 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
@@ -39,6 +39,7 @@ SourceFiles
#include "labelList.H"
#include "DynamicList.H"
+#include "point.H"
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,6 +67,10 @@ class edgeCollapser
//- For every point -1 or region number
labelList pointRegion_;
+ //- Actual location of the point to collapse to for every region master
+ // point. This will be forced to be consistent across processors
+ DynamicList<point> pointRegionMasterLocation_;
+
//- -1 or master vertex for region number
DynamicList<label> pointRegionMaster_;