diff --git a/applications/test/Circulator/Make/files b/applications/test/Circulator/Make/files
new file mode 100644
index 0000000000..b6a4e63dc2
--- /dev/null
+++ b/applications/test/Circulator/Make/files
@@ -0,0 +1,3 @@
+Test-Circulator.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-Circulator
diff --git a/applications/test/Circulator/Make/options b/applications/test/Circulator/Make/options
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/applications/test/Circulator/Test-Circulator.C b/applications/test/Circulator/Test-Circulator.C
new file mode 100644
index 0000000000..19a21f6b79
--- /dev/null
+++ b/applications/test/Circulator/Test-Circulator.C
@@ -0,0 +1,268 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 aCirc(a);
+ const_circulator bCirc(b);
+
+ // Rotate face b until its element matches the starting element of face a.
+ do
+ {
+ if (aCirc() == bCirc())
+ {
+ // Set bCirc fulcrum to its iterator and increment the iterators
+ bCirc.setFulcrumToIterator();
+ ++aCirc;
+ ++bCirc;
+
+ break;
+ }
+ } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
+
+ // Look forwards around the faces for a match
+ do
+ {
+ if (aCirc() != bCirc())
+ {
+ break;
+ }
+ }
+ while
+ (
+ aCirc.circulate(CirculatorBase::CLOCKWISE),
+ bCirc.circulate(CirculatorBase::CLOCKWISE)
+ );
+
+ // If the circulator has stopped then faces a and b matched.
+ if (!aCirc.circulate())
+ {
+ return 1;
+ }
+ else
+ {
+ // Reset the circulators back to their fulcrum
+ aCirc.setIteratorToFulcrum();
+ bCirc.setIteratorToFulcrum();
+ ++aCirc;
+ --bCirc;
+ }
+
+ // Look backwards around the faces for a match
+ do
+ {
+ if (aCirc() != bCirc())
+ {
+ break;
+ }
+ }
+ while
+ (
+ aCirc.circulate(CirculatorBase::CLOCKWISE),
+ bCirc.circulate(CirculatorBase::ANTICLOCKWISE)
+ );
+
+ // If the circulator has stopped then faces a and b matched.
+ if (!aCirc.circulate())
+ {
+ return -1;
+ }
+
+ return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+ Info<< "Test the implementation of a circular iterator" << nl << endl;
+
+ Info<< "Test const circulator. First go forwards, then backwards."
+ << nl << endl;
+
+ face f(identity(4));
+
+ const_circulator cStart(f);
+
+ if (cStart.size()) do
+ {
+ Info<< "Iterate forwards over face (prev/curr/next) : "
+ << cStart.prev() << " / " << cStart() << " / " << cStart.next()
+ << endl;
+
+ } while (cStart.circulate(CirculatorBase::CLOCKWISE));
+
+ if (cStart.size()) do
+ {
+ Info<< "Iterate backwards over face : " << cStart() << endl;
+
+ } while (cStart.circulate(CirculatorBase::ANTICLOCKWISE));
+
+
+ Info<< nl << nl << "Test non-const circulator" << nl << endl;
+
+ circulator cStart2(f);
+
+ Info<< "Face before : " << f << endl;
+
+ if (cStart2.size()) do
+ {
+ Info<< "Iterate forwards over face (prev/curr/next) : "
+ << cStart2.prev() << " / " << cStart2() << " / " << cStart2.next()
+ << endl;
+
+ } while (cStart2.circulate(CirculatorBase::CLOCKWISE));
+
+ if (cStart2.size()) do
+ {
+ Info<< "Iterate forwards over face, adding 1 to each element : "
+ << cStart2();
+
+ cStart2() += 1;
+
+ Info<< " -> " << cStart2() << endl;
+ } while (cStart2.circulate(CirculatorBase::CLOCKWISE));
+
+ Info<< "Face after : " << f << endl;
+
+
+ Info<< nl << nl << "Compare two faces: " << endl;
+ face a(identity(5));
+ Info<< "Compare " << a << " and " << a << " Match = " << compare(a, a)
+ << endl;
+
+ face b(reverseList(a));
+ Info<< "Compare " << a << " and " << b << " Match = " << compare(a, b)
+ << endl;
+
+ face c(a);
+ c[4] = 3;
+ Info<< "Compare " << a << " and " << c << " Match = " << compare(a, c)
+ << endl;
+
+ face d(rotateList(a, 2));
+ Info<< "Compare " << a << " and " << d << " Match = " << compare(a, d)
+ << endl;
+
+ face g(labelList(5, 1));
+ face h(g);
+ Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+ << endl;
+
+ g[0] = 2;
+ h[3] = 2;
+ Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+ << endl;
+
+ g[4] = 3;
+ h[4] = 3;
+ Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+ << endl;
+
+ face face1(identity(1));
+ Info<< "Compare " << face1 << " and " << face1
+ << " Match = " << compare(face1, face1) << endl;
+
+ Info<< nl << nl << "Zero face" << nl << endl;
+
+ face fZero;
+ circulator cZero(fZero);
+
+ if (cZero.size()) do
+ {
+ Info<< "Iterate forwards over face : " << cZero() << endl;
+
+ } while (cZero.circulate(CirculatorBase::CLOCKWISE));
+
+ fZero = face(identity(5));
+
+ // circulator was invalidated so reset
+ cZero = circulator(fZero);
+
+ do
+ {
+ Info<< "Iterate forwards over face : " << cZero() << endl;
+
+ } while (cZero.circulate(CirculatorBase::CLOCKWISE));
+
+
+ Info<< nl << nl << "Simultaneously go forwards/backwards over face " << f
+ << nl << endl;
+
+ const_circulator circForward(f);
+ const_circulator circBackward(f);
+
+ if (circForward.size() && circBackward.size()) do
+ {
+ Info<< "Iterate over face forwards : " << circForward()
+ << ", backwards : " << circBackward() << endl;
+ }
+ while
+ (
+ circForward.circulate(CirculatorBase::CLOCKWISE),
+ circBackward.circulate(CirculatorBase::ANTICLOCKWISE)
+ );
+
+ Info<< "\nEnd\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/ListOps/Make/files b/applications/test/ListOps/Make/files
new file mode 100644
index 0000000000..645600a709
--- /dev/null
+++ b/applications/test/ListOps/Make/files
@@ -0,0 +1,2 @@
+Test-ListOps.C
+EXE = $(FOAM_USER_APPBIN)/Test-ListOps
diff --git a/applications/test/ListOps/Make/options b/applications/test/ListOps/Make/options
new file mode 100644
index 0000000000..9e79e9d733
--- /dev/null
+++ b/applications/test/ListOps/Make/options
@@ -0,0 +1,3 @@
+EXE_INC = /*-DFULLDEBUG -O0 -g*/ \
+
+EXE_LIBS = \
diff --git a/applications/test/ListOps/Test-ListOps.C b/applications/test/ListOps/Test-ListOps.C
new file mode 100644
index 0000000000..f6e5695164
--- /dev/null
+++ b/applications/test/ListOps/Test-ListOps.C
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 forwardRotate(identity(5));
+ face testFace(identity(4));
+
+ for (label i = 0; i < 8; ++i)
+ {
+ Info<< "Rotate forward by " << i << " : "
+ << rotateList(forwardRotate, i) << endl;
+ }
+
+ for (label i = 0; i < 8; ++i)
+ {
+ Info<< "Rotate backward by " << i << " : "
+ << rotateList(forwardRotate, -i) << endl;
+ }
+
+ Info<< nl << "Face : " << testFace << endl;
+ Info<< "Rotate by 2 : " << rotateList(testFace, 2) << endl;
+ inplaceRotateList(testFace, -6);
+ Info<< "Rotate inplace by -6 : " << testFace << nl << endl;
+
+ Info<< "Test inplace rotate : " << forwardRotate << endl;
+ inplaceRotateList(forwardRotate, 2);
+ Info<< "Rotate to the right by 2 : " << forwardRotate << endl;
+ inplaceRotateList(forwardRotate, -2);
+ Info<< "Rotate to the left by 2 : " << forwardRotate << endl;
+
+ List subRotate(identity(10));
+ SubList subL(subRotate, 5, 3);
+
+ Info<< "Test inplace rotate on sublist : " << subRotate << endl;
+ inplaceRotateList(subL, 3);
+ Info<< "Rotate to the right by 3 : " << subRotate << endl;
+ inplaceRotateList(subL, -8);
+ Info<< "Rotate to the left by 3 : " << subRotate << endl;
+
+ Info<< nl << nl << "Test Reversing:" << nl << endl;
+
+ Info<< "List : " << identity(5) << endl;
+ Info<< "Reverse : " << reverseList(identity(5)) << endl;
+ Info<< "List : " << identity(6) << endl;
+ Info<< "Reverse : " << reverseList(identity(6)) << nl << endl;
+
+ List test1(identity(5));
+ Info<< "List : " << test1 << endl;
+ inplaceReverseList(test1);
+ Info<< "Inplace Reverse : " << test1 << nl << endl;
+
+ List test2(identity(6));
+ Info<< "List : " << test2 << endl;
+ inplaceReverseList(test2);
+ Info<< "Inplace Reverse : " << test2 << nl << endl;
+
+ face test3(identity(6));
+ Info<< "Face : " << test3 << endl;
+ inplaceReverseList(test3);
+ Info<< "Inplace Reverse : " << test3 << nl << endl;
+
+ FixedList test4(identity(6));
+ Info<< "FixedList : " << test4 << endl;
+ inplaceReverseList(test4);
+ Info<< "Inplace Reverse : " << test4 << nl << endl;
+
+ List test5(identity(9));
+ SubList test5SubList(test5, 4, 3);
+ Info<< "List : " << test5 << endl;
+ inplaceReverseList(test5SubList);
+ Info<< "Reverse Sublist between 3 and 6 : " << test5 << endl;
+
+ Info<< "\nEnd\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/dynamicIndexedOctree/Test-dynamicIndexedOctree.C b/applications/test/dynamicIndexedOctree/Test-dynamicIndexedOctree.C
index b84b39c33b..0ccb8e377f 100644
--- a/applications/test/dynamicIndexedOctree/Test-dynamicIndexedOctree.C
+++ b/applications/test/dynamicIndexedOctree/Test-dynamicIndexedOctree.C
@@ -85,8 +85,8 @@ int main(int argc, char *argv[])
(
dynamicTreeDataPoint(pointList),
overallBb, // overall search domain
- 20, // max levels; n/a
- 100, // maximum ratio of cubes v.s. cells
+ 20, // max levels
+ 100, // maximum ratio of cubes v.s. cells
100.0 // max. duplicity; n/a since no bounding boxes.
);
diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/files b/applications/utilities/mesh/advanced/collapseEdges/Make/files
index a15838abe8..d89ca6e737 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/files
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/files
@@ -1,4 +1,3 @@
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 d1efa61fd5..987eae5ed7 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/options
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/options
@@ -1,8 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
- -I$(LIB_SRC)/finiteVolume/lnInclude \
- -IpointEdgeCollapse
+ -I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-ldynamicMesh \
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
new file mode 100644
index 0000000000..78905475e3
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
@@ -0,0 +1,85 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: http://www.openfoam.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+ version 2.0;
+ format ascii;
+
+ root "";
+ case "";
+ instance "";
+ local "";
+
+ class dictionary;
+ object collapseDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+collapseEdgesCoeffs
+{
+ // Edges shorter than this absolute value will be merged
+ minimumEdgeLength 1e-8;
+
+ // The maximum angle between two edges that share a point attached to
+ // no other edges
+ maximumMergeAngle 30;
+
+ // The amount that minimumEdgeLength will be reduced by for each
+ // edge if that edge's collapse generates a poor quality face
+ reductionFactor 0.5;
+}
+
+
+collapseFacesCoeffs
+{
+ // The initial face length factor
+ initialFaceLengthFactor 0.5;
+
+ // The amount that initialFaceLengthFactor will be reduced by for each
+ // face if its collapse generates a poor quality face
+ reductionFactor $initialFaceLengthFactor;
+
+ // If the face can't be collapsed to an edge, and it has a span less than
+ // the target face length multiplied by this coefficient, collapse it
+ // to a point.
+ maxCollapseFaceToPointSideLengthCoeff 0.3;
+
+ // Allow early collapse of edges to a point
+ allowEarlyCollapseToPoint on;
+
+ // Fraction to premultiply maxCollapseFaceToPointSideLengthCoeff by if
+ // allowEarlyCollapseToPoint is enabled
+ allowEarlyCollapseCoeff 0.2;
+
+ // Defining how close to the midpoint (M) of the projected
+ // vertices line a projected vertex (X) can be before making this
+ // an invalid edge collapse
+ //
+ // X---X-g----------------M----X-----------g----X--X
+ //
+ // Only allow a collapse if all projected vertices are outwith
+ // guardFraction (g) of the distance form the face centre to the
+ // furthest vertex in the considered direction
+ guardFraction 0.1;
+}
+
+
+meshQualityCoeffs
+{
+ // Name of the dictionary that has the mesh quality coefficients used
+ // by motionSmoother::checkMesh
+ meshQualityCoeffDict meshQualityDict;
+
+ // Maximum number of outer iterations is mesh quality checking is enabled
+ maximumIterations 30;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index c0523bf43b..132a602b1a 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
@@ -27,1425 +27,96 @@ Description
- collapse short edges. Length of edges to collapse provided as argument.
- merge two edges if they are in line. Maximum angle provided as argument.
- remove unused points.
+ - collapse faces:
+ - with small areas to a single point
+ - that have a high aspect ratio (i.e. sliver face) to a single edge
- Optionally removes cells. Can remove faces and points but does not check
- for nonsense resulting topology.
+ Optionally checks the resulting mesh for bad faces and reduces the desired
+ face length factor for those faces attached to the bad faces.
When collapsing an edge with one point on the boundary it will leave
the boundary point intact. When both points inside it chooses random. When
both points on boundary random again.
Usage
- - collapseEdges [OPTION]
-
- \param -allowCellCollapse \n
- Allow collapsing of cells to a single point
-
- \param -meshQuality \n
- Only collapse if not exceeding user defined (in \a system/meshQualityDict)
- quality settings
+ - collapseEdges [OPTION]
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
+#include "timeSelector.H"
#include "polyTopoChange.H"
#include "fvMesh.H"
-#include "mapPolyMesh.H"
-#include "mathematicalConstants.H"
-#include "PackedBoolList.H"
-#include "unitConversion.H"
-#include "globalMeshData.H"
-#include "globalIndex.H"
-#include "PointEdgeWave.H"
-#include "pointEdgeCollapse.H"
-#include "motionSmoother.H"
-#include "removePoints.H"
+#include "polyMeshFilter.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-void filterFace
-(
- const label faceI,
- face& f,
- const List& allPointInfo,
- const Map >& collapseStrings
-)
-{
- label newFp = 0;
-
- face oldFace = f;
-
- forAll(f, fp)
- {
- label pointI = f[fp];
-
- label collapsePoint = allPointInfo[pointI].collapseIndex();
-
- if (collapseStrings.found(collapsePoint))
- {
- collapsePoint = collapseStrings[collapsePoint][0];
- }
-
- if (collapsePoint == -1)
- {
- WarningIn
- (
- "filterFace"
- "(const label, face&, const List&)"
- )
- << "Point " << pointI << " was not visited by PointEdgeWave"
- << endl;
- }
- else if (collapsePoint == -2)
- {
- f[newFp++] = pointI;
- }
- else
- {
- if (findIndex(SubList(f, newFp), collapsePoint) == -1)
- {
- f[newFp++] = collapsePoint;
- }
- }
- }
-
-
- // 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(SubList(f, 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);
-}
-
-
-bool setRefinement
-(
- const polyMesh& mesh,
- polyTopoChange& meshMod,
- const List& allPointInfo
-)
-{
- 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 from collapseIndex(=global master point) to set of points
- Map > collapseStrings;
- {
- // 1. Count elements per collapseIndex
- Map nPerIndex(mesh.nPoints()/10);
- forAll(allPointInfo, pointI)
- {
- const label collapseIndex = allPointInfo[pointI].collapseIndex();
-
- if (collapseIndex != -1 && collapseIndex != -2)
- {
- Map::iterator fnd = nPerIndex.find(collapseIndex);
- if (fnd != nPerIndex.end())
- {
- fnd()++;
- }
- else
- {
- nPerIndex.insert(collapseIndex, 1);
- }
- }
- }
-
- // 2. Size
- collapseStrings.resize(2*nPerIndex.size());
- forAllConstIter(Map, nPerIndex, iter)
- {
- collapseStrings.insert(iter.key(), DynamicList(iter()));
- }
-
- // 3. Fill
- 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(allPointInfo, pI)
- {
- const pointEdgeCollapse& pec = allPointInfo[pI];
-
- if (pec.collapseIndex() == -1)
- {
- nUnvisited++;
- }
- else if (pec.collapseIndex() == -2)
- {
- nUncollapsed++;
- }
- else if (pec.collapseIndex() >= 0)
- {
- nCollapsed++;
- }
- }
-
- label nPoints = allPointInfo.size();
-
- reduce(nPoints, sumOp());
- reduce(nUnvisited, sumOp());
- reduce(nUncollapsed, sumOp());
- reduce(nCollapsed, sumOp());
-
- 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;
-
- 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])
- {
- 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)
- {
- Pout<< "Cell:" << cellI
- << " uses faces:" << cFaces
- << " of which too many are marked for removal:"
- << endl
- << " ";
- forAll(cFaces, j)
- {
- if (newFaces[cFaces[j]].size() < 3)
- {
- Pout<< ' '<< cFaces[j];
- }
- }
- Pout<< 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);
-
- forAll(cellRemoved, cellI)
- {
- if (cellRemoved[cellI])
- {
- meshMod.removeCell(cellI, -1);
- }
- }
-
- // Remove faces
- forAll(newFaces, faceI)
- {
- const face& f = newFaces[faceI];
-
- if (f.size() < 3)
- {
- meshMod.removeFace(faceI, -1);
- meshChanged = true;
-
- // Mark face as been done.
- doneFace[faceI] = true;
- }
- else
- {
- // Kept face. Mark vertices
- forAll(f, fp)
- {
- usedPoint[f[fp]] = true;
- }
- }
- }
-
- // Remove unused vertices that have not been marked for removal already
- forAll(usedPoint, pointI)
- {
- if (!usedPoint[pointI])
- {
- 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];
-
- forAll(changedFaces, changedFaceI)
- {
- label faceI = changedFaces[changedFaceI];
-
- if (!doneFace[faceI])
- {
- doneFace[faceI] = true;
-
- // Get current zone info
- label zoneID = faceZones.whichZone(faceI);
-
- bool zoneFlip = false;
-
- if (zoneID >= 0)
- {
- 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 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;
-}
-
-
-// Return master point edge needs to be collapsed to (or -1)
-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[e0] == 0)
- {
- if (boundaryPoint[e1] == 0)
- {
- // Both points on boundary. Choose one to collapse to.
- // Note: should look at feature edges/points!
- masterPoint = e0;
- }
- else
- {
- masterPoint = e0;
- }
- }
- else
- {
- if (boundaryPoint[e1] == 0)
- {
- masterPoint = e1;
- }
- else
- {
- // None on boundary. Choose arbitrary.
- // Note: should look at geometry?
- masterPoint = e0;
- }
- }
-
- return masterPoint;
-}
-
-
-// Create consistent set of collapses.
-// collapseEdge : per edge:
-// -1 : do not collapse
-// 0 : collapse to start
-// 1 : collapse to end
-// Note: collapseEdge has to be parallel consistent (in orientation)
-label syncCollapse
-(
- const polyMesh& mesh,
- const globalIndex& globalStrings,
- const labelList& collapseEdge,
- List& allPointInfo
-)
-{
- const edgeList& edges = mesh.edges();
- const pointField& points = mesh.points();
-
- label nCollapsed = 0;
-
- DynamicList initPoints(mesh.nPoints());
- DynamicList initPointInfo(mesh.nPoints());
- allPointInfo.resize(mesh.nPoints());
-
- // Initialise edges to no collapse
- List allEdgeInfo
- (
- mesh.nEdges(),
- pointEdgeCollapse(vector::zero, -1)
- );
-
- // Mark selected edges for collapse
- forAll(edges, edgeI)
- {
- const edge& e = edges[edgeI];
-
- if (collapseEdge[edgeI] != -1)
- {
- label masterPointI = e[collapseEdge[edgeI]];
-
- const pointEdgeCollapse pec
- (
- points[masterPointI],
- globalStrings.toGlobal(masterPointI)
- );
-
- // Mark as collapsable but with nonsense master so it gets
- // overwritten and starts an update wave
- allEdgeInfo[edgeI] = pointEdgeCollapse
- (
- points[masterPointI],
- labelMax
- );
-
- initPointInfo.append(pec);
- initPoints.append(e.start());
-
- initPointInfo.append(pec);
- initPoints.append(e.end());
-
- nCollapsed++;
- }
- }
-
- PointEdgeWave collapsePropagator
- (
- mesh,
- initPoints,
- initPointInfo,
- allPointInfo,
- allEdgeInfo,
- mesh.globalData().nTotalPoints() // Maximum number of iterations
- );
-
- return nCollapsed;
-}
-
-
-void syncCollapseEdge(const polyMesh& mesh, labelList& collapseEdge)
-{
- // 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 = mesh.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(map.constructSize());
-
- forAll(coupledMeshEdges, eI)
- {
- const label meshEdgeIndex = coupledMeshEdges[eI];
-
- cppEdgeData[eI] = collapseEdge[meshEdgeIndex];
-
- if
- (
- (collapseEdge[meshEdgeIndex] != -1)
- && (meshToPatchSameOrientation[eI] != cppOrientation[eI])
- )
- {
- cppEdgeData[eI] = 1-cppEdgeData[eI];
- }
- }
-
-
- // 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()
- );
-
-
- forAll(coupledMeshEdges, eI)
- {
- const label meshEdgeIndex = coupledMeshEdges[eI];
-
- collapseEdge[meshEdgeIndex] = cppEdgeData[eI];
-
- if
- (
- (cppEdgeData[eI] != -1)
- && (meshToPatchSameOrientation[eI] != cppOrientation[eI])
- )
- {
- collapseEdge[meshEdgeIndex] = 1-collapseEdge[meshEdgeIndex];
- }
- }
-}
-
-
-// Mark (in collapseEdge) any edges to collapse
-label collapseSmallEdges
-(
- const polyMesh& mesh,
- const labelList& boundaryPoint,
- const scalarField& minEdgeLen,
- labelList& collapseEdge
-)
-{
- // Work out which edges to collapse
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- label nCollapsed = 0;
-
- forAll(mesh.edges(), edgeI)
- {
- if (collapseEdge[edgeI] == -1)
- {
- const edge& e = mesh.edges()[edgeI];
-
- if (e.mag(mesh.points()) < minEdgeLen[edgeI])
- {
- label masterPointI = edgeMaster(boundaryPoint, false, e);
- if (masterPointI == e[0])
- {
- collapseEdge[edgeI] = 0;
- }
- else
- {
- collapseEdge[edgeI] = 1;
- }
- nCollapsed++;
- }
- }
- }
- return nCollapsed;
-}
-
-
-// Mark (in collapseEdge) any edges to merge
-label mergeEdges
-(
- const polyMesh& mesh,
- const scalar maxCos,
- const labelList& boundaryPoint,
- const scalarField& minEdgeLen,
- labelList& collapseEdge
-)
-{
- const edgeList& edges = mesh.edges();
- const labelListList& pointEdges = mesh.pointEdges();
-
- // Point removal engine
- removePoints pointRemover(mesh, false);
-
- // Find out points that can be deleted
- boolList pointCanBeDeleted;
- label nTotRemove = pointRemover.countPointUsage(maxCos, pointCanBeDeleted);
-
-
- // Rework point-to-remove into edge-to-collapse.
-
- label nCollapsed = 0;
-
- if (nTotRemove > 0)
- {
- forAll(pointEdges, pointI)
- {
- if (pointCanBeDeleted[pointI])
- {
- const labelList& pEdges = pointEdges[pointI];
-
- if (pEdges.size() == 2)
- {
- // Always the case?
-
- label e0 = pEdges[0];
- label e1 = pEdges[1];
-
- if
- (
- collapseEdge[e0] == -1
- && minEdgeLen[e0] >= 0
- && collapseEdge[e1] == -1
- && minEdgeLen[e1] >= 0
- )
- {
- // Get the two vertices on both sides of the point
- label leftV = edges[e0].otherVertex(pointI);
- label rightV = edges[e1].otherVertex(pointI);
-
- // Can collapse pointI onto either leftV or rightV.
- // Preferentially choose an internal point to hopefully
- // give less distortion
-
- if (boundaryPoint[leftV] == -1)
- {
- collapseEdge[e0] = findIndex(edges[e0], leftV);
- }
- else
- {
- collapseEdge[e1] = findIndex(edges[e1], rightV);
- }
- }
- }
- }
- }
- }
- return nCollapsed;
-}
-
-
-// Make consistent set of collapses that does not collapse any cells
-label consistentCollapse
-(
- const bool allowCellCollapse,
- const polyMesh& mesh,
- const globalIndex& globalStrings,
- labelList& collapseEdge,
- List& allPointInfo
-)
-{
- // Make sure we don't collapse cells
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- while (true)
- {
- // Sync collapseEdge
- syncCollapseEdge(mesh, collapseEdge);
-
-
- // Get collapsed faces
-
- label nAdditionalCollapsed = 0;
-
- PackedBoolList isCollapsedFace(mesh.nFaces());
- forAll(mesh.faceEdges(), faceI)
- {
- const labelList& fEdges = mesh.faceEdges()[faceI];
-
- // Count number of remaining edges
- label nEdges = fEdges.size();
- forAll(fEdges, fEdgeI)
- {
- label edgeI = fEdges[fEdgeI];
- if (collapseEdge[edgeI] != -1)
- {
- nEdges--;
- }
- }
-
- if (nEdges < 3)
- {
- // Face is collapsed.
- isCollapsedFace[faceI] = 1;
-
- if (nEdges == 1)
- {
- // Cannot collapse face down to single edge.
-
- //- Option 1: collapse remaining edge as well. However
- // if this edge is on the coupled processor patch this
- // logic clashes with that of syncCollapseEdge
- // (do not collapse if any not collapse)
- //forAll(fEdges, fEdgeI)
- //{
- // label edgeI = fEdges[fEdgeI];
- // if (collapseEdge[edgeI] == -1)
- // {
- // collapseEdge[edgeI] = 0;
- // nAdditionalCollapsed++;
- // }
- //}
-
- //- Option 2: uncollapse this face.
- forAll(fEdges, fEdgeI)
- {
- label edgeI = fEdges[fEdgeI];
- if (collapseEdge[edgeI] != -1)
- {
- collapseEdge[edgeI] = -1;
- nAdditionalCollapsed++;
- }
- }
- }
- }
- }
-
- //Pout<< "nAdditionalCollapsed : " << nAdditionalCollapsed << endl;
-
-
- label nUncollapsed = 0;
-
- if (!allowCellCollapse)
- {
- // Check collapsed cells
-
- forAll(mesh.cells(), cellI)
- {
- const cell& cFaces = mesh.cells()[cellI];
- label nFaces = cFaces.size();
- forAll(cFaces, i)
- {
- label faceI = cFaces[i];
- if (isCollapsedFace[faceI])
- {
- nFaces--;
- if (nFaces < 4)
- {
- // Unmark this face for collapse
- const labelList& fEdges = mesh.faceEdges()[faceI];
-
- forAll(fEdges, fEdgeI)
- {
- label edgeI = fEdges[fEdgeI];
- if (collapseEdge[edgeI] != -1)
- {
- collapseEdge[edgeI] = -1;
- nUncollapsed++;
- }
- }
-
- // Uncollapsed this face.
- isCollapsedFace[faceI] = 0;
- nFaces++;
- }
- }
- }
- }
- //Pout<< "** disallowing cells : nUncollapsed : "
- // << nUncollapsed << endl;
- }
-
-
- if
- (
- returnReduce(nUncollapsed+nAdditionalCollapsed, sumOp())
- == 0
- )
- {
- break;
- }
- }
-
-
- // Create consistent set of collapses
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // Note: requires collapseEdge to be synchronised. (above loop makes sure
- // of that)
-
- return syncCollapse(mesh, globalStrings, collapseEdge, allPointInfo);
-}
-
-
-// Check mesh and mark points on faces in error
-// Returns boolList with points in error set
-PackedBoolList checkMeshQuality(const polyMesh& mesh)
-{
- //mesh.checkMesh(true);
-
- IOdictionary meshQualityDict
- (
- IOobject
- (
- "meshQualityDict",
- mesh.time().system(),
- mesh.time(),
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- );
-
- labelHashSet badFaces(mesh.nFaces()/100);
- DynamicList checkFaces(mesh.nFaces());
-
- const vectorField& fAreas = mesh.faceAreas();
-
- scalar faceAreaLimit = SMALL;
-
- forAll(fAreas, fI)
- {
- if (mag(fAreas[fI]) > faceAreaLimit)
- {
- checkFaces.append(fI);
- }
- }
-
- motionSmoother::checkMesh
- (
- false,
- mesh,
- meshQualityDict,
- checkFaces,
- badFaces
- );
-
- label nBadFaces = badFaces.size();
- reduce(nBadFaces, sumOp());
-
- Info<< nl << "Number of bad faces : " << nBadFaces << endl;
-
- PackedBoolList isErrorPoint(mesh.nPoints());
- forAllConstIter(labelHashSet, badFaces, iter)
- {
- const face& f = mesh.faces()[iter.key()];
-
- forAll(f, pI)
- {
- isErrorPoint[f[pI]] = 1;
- }
- }
-
- syncTools::syncPointList
- (
- mesh,
- isErrorPoint,
- orEqOp(),
- 0
- );
-
- return isErrorPoint;
-}
-
-
-// Check mesh with collapses (newMesh), updates minEdgeLen, nFrozenEdges
-void checkMeshAndFreezeEdges
-(
- const polyMesh& newMesh,
- const labelList& oldToNewMesh,
- const polyMesh& oldMesh,
- scalarField& minEdgeLen,
- label& nFrozenEdges
-)
-{
- PackedBoolList isErrorPoint = checkMeshQuality(newMesh);
-
- forAll(oldMesh.edges(), edgeI)
- {
- const edge& e = oldMesh.edges()[edgeI];
- label newStart = oldToNewMesh[e[0]];
- label newEnd = oldToNewMesh[e[1]];
-
- if
- (
- (newStart >= 0 && isErrorPoint[newStart])
- || (newEnd >= 0 && isErrorPoint[newEnd])
- )
- {
- // Gradual decrease? For now just hard disable
- if (minEdgeLen[edgeI] > -GREAT/2)
- {
- minEdgeLen[edgeI] = -GREAT;
- nFrozenEdges++;
- }
- }
- }
-}
-
-
-// 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();
-
-
- labelList boundaryPoint(mesh.nPoints(), -1);
-
- // Get all points on a boundary
- label nIntFaces = mesh.nInternalFaces();
- for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
- {
- const face& f = faces[faceI];
-
- forAll(f, fp)
- {
- boundaryPoint[f[fp]] = 0;
- }
- }
-
- // Get all processor boundary points and the processor patch label
- // that they are on.
- forAll(bMesh, patchI)
- {
- const polyPatch& patch = bMesh[patchI];
-
- if (isA(patch))
- {
- if (patchI == 0)
- {
- // We mark 'normal' boundary points with 0 so make sure this
- // coupled patch is not 0.
- FatalErrorIn("findBoundaryPoints(const polyMesh&)")
- << "Your patches should have non-coupled ones before any"
- << " coupled ones. Current patches " << bMesh.names()
- << exit(FatalError);
- }
-
- forAll(patch, fI)
- {
- const face& f = patch[fI];
-
- forAll(f, fp)
- {
- boundaryPoint[f[fp]] = patchI;
- }
- }
- }
- }
-
- return boundaryPoint;
-}
-
-
// Main program:
int main(int argc, char *argv[])
{
+ timeSelector::addOptions(true, false);
argList::addNote
(
- "Merges small and in-line edges.\n"
- "Collapses faces and optionally cells to a point."
+ "Collapses small edges to a point.\n"
+ "Optionally collapse small faces to a point and thin faces to an edge."
+ );
+
+ argList::addBoolOption
+ (
+ "collapseFaces",
+ "Collapse small and sliver faces as well as small edges"
);
# include "addOverwriteOption.H"
-
- argList::validArgs.append("edge length [m]");
- argList::validArgs.append("merge angle (degrees)");
- argList::addBoolOption
- (
- "allowCellCollapse",
- "Allow collapsing of cells to a single point"
- );
- argList::addBoolOption
- (
- "meshQuality",
- "Only collapse if not exceeding given meshQualityDict limits"
- );
-
-
-
# include "setRootCase.H"
# include "createTime.H"
+
runTime.functionObjects().off();
-# include "createPolyMesh.H"
+ instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
+
+# include "createMesh.H"
+
const word oldInstance = mesh.pointsInstance();
- scalar minLen = args.argRead(1);
- const scalar angle = args.argRead(2);
-
- const bool allowCellCollapse = args.optionFound("allowCellCollapse");
const bool overwrite = args.optionFound("overwrite");
- const bool checkQuality = args.optionFound("meshQuality");
- scalar maxCos = Foam::cos(degToRad(angle));
+ const bool collapseFaces = args.optionFound("collapseFaces");
- Info<< "Merging:" << nl
- << " edges with length less than " << minLen << " meters" << nl
- << " edges split by a point with edges in line to within " << angle
- << " degrees" << nl
- << endl;
-
- if (allowCellCollapse)
+ forAll(timeDirs, timeI)
{
- Info<< "Allowing collapse of cells down to a point." << nl
- << endl;
- }
- else
- {
- Info<< "Disallowing collapse of cells down to a point." << nl
- << endl;
- }
+ runTime.setTime(timeDirs[timeI], timeI);
+ Info<< "Time = " << runTime.timeName() << endl;
- if (checkQuality)
- {
- Info<< "Selectively disabling wanted collapses until resulting quality"
- << " satisfies constraints in system/meshQualityDict" << nl
- << endl;
- }
+ polyMeshFilter meshFilter(mesh);
+ // newMesh will be empty until it is filtered
+ const autoPtr& newMesh = meshFilter.filteredMesh();
-
- // To mark master of collapes
- globalIndex globalStrings(mesh.nPoints());
-
-
- // Local collapse length. Any edge below this length gets (attempted)
- // collapsed. Use either aa gradually decreasing value
- // (so from minLen to e.g. 0.5*minLen) or a hard limit (GREAT)
- scalarField minEdgeLen(mesh.nEdges(), minLen);
- label nFrozenEdges = 0;
-
-
-
- // Initial mesh check
- // ~~~~~~~~~~~~~~~~~~
- // Do not allow collapses in regions of error.
- // Updates minEdgeLen, nFrozenEdges
- if (checkQuality)
- {
- checkMeshAndFreezeEdges
- (
- mesh,
- identity(mesh.nPoints()),
- mesh,
- minEdgeLen,
- nFrozenEdges
- );
- Info<< "Initial frozen edges "
- << returnReduce(nFrozenEdges, sumOp())
- << " out of " << returnReduce(mesh.nEdges(), sumOp())
- << endl;
- }
-
-
- // Mark points on boundary
- const labelList boundaryPoint = findBoundaryPoints(mesh);
-
- // Copy of current set of topology changes. Used to generate final mesh.
- polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
-
- // Keep track of whether mesh has changed at all
- bool meshChanged = false;
-
-
-
- // Main loop
- // ~~~~~~~~~
- // It tries and do some collapses, checks the resulting mesh and
- // 'freezes' some edges (by marking in minEdgeLen) and tries again.
- // This will iterate ultimately to the situation where every edge is
- // frozen and nothing gets collapsed.
-
- do
- {
- // Per edge collapse status:
- // -1 : not collapsed
- // 0 : collapse to start
- // 1 : collapse to end
- labelList collapseEdge(mesh.nEdges(), -1);
-
-
- // Work out which edges to collapse
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // This is by looking at minEdgeLen (to avoid frozen edges)
- // and marking in collapseEdge.
-
- // Small edges
- label nSmallCollapsed = collapseSmallEdges
- (
- mesh,
- boundaryPoint,
- minEdgeLen,
- collapseEdge
- );
- reduce(nSmallCollapsed, sumOp());
-
- Info<< indent << "Collapsing " << nSmallCollapsed
- << " small edges" << endl;
-
-
- // Inline edges
- label nMerged = mergeEdges
- (
- mesh,
- maxCos,
- boundaryPoint,
- minEdgeLen,
- collapseEdge
- );
-
- reduce(nMerged, sumOp());
- Info<< indent << "Collapsing " << nMerged << " in line edges"
- << endl;
-
-
- // Merge edge collapses into consistent collapse-network. Make sure
- // no cells get collapsed.
- List allPointInfo;
- label nLocalCollapse = consistentCollapse
- (
- allowCellCollapse,
- mesh,
- globalStrings,
- collapseEdge,
- allPointInfo
- );
-
- reduce(nLocalCollapse, sumOp());
- Info<< "nLocalCollapse = " << nLocalCollapse << endl;
-
- if (nLocalCollapse == 0)
+ // Filter small edges only. This reduces the number of faces so that
+ // the face filtering is sped up.
+ label nBadFaces = meshFilter.filterEdges(0);
{
- break;
+ polyTopoChange meshMod(newMesh);
+
+ meshMod.changeMesh(mesh, false);
}
-
- // There are collapses so mesh will get changed
- meshChanged = true;
-
-
- // Apply collapses to current mesh
- polyTopoChange meshMod(mesh);
-
- // Insert mesh refinement into polyTopoChange.
- setRefinement(mesh, meshMod, allPointInfo);
-
- // Do all changes
- Info<< indent << "Applying changes to the mesh" << nl
- //<< decrIndent
- << endl;
-
- savedMeshMod = meshMod;
-
- autoPtr newMeshPtr;
- autoPtr mapPtr = meshMod.makeMesh
- (
- newMeshPtr,
- IOobject
- (
- mesh.name(),
- mesh.instance(),
- mesh.time(), // register with runTime
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- mesh,
- true // parallel sync
- );
-
- fvMesh& newMesh = newMeshPtr();
- const mapPolyMesh& map = mapPtr();
-
- // Update fields
- newMesh.updateMesh(map);
- if (map.hasMotionPoints())
+ if (collapseFaces)
{
- newMesh.movePoints(map.preMotionPoints());
+ // Filter faces. Pass in the number of bad faces that are present
+ // from the previous edge filtering to use as a stopping criterion.
+ meshFilter.filter(nBadFaces);
+ {
+ polyTopoChange meshMod(newMesh);
+
+ meshMod.changeMesh(mesh, false);
+ }
}
-
-
-
- // If no checks needed exit.
- if (!checkQuality)
- {
- break;
- }
-
- // Mesh check
- // ~~~~~~~~~~~~~~~~~~
- // Do not allow collapses in regions of error.
- // Updates minEdgeLen, nFrozenEdges
- label nOldFrozenEdges = returnReduce(nFrozenEdges, sumOp());
- checkMeshAndFreezeEdges
- (
- newMesh,
- map.reversePointMap(),
- mesh,
- minEdgeLen,
- nFrozenEdges
- );
- label nNewFrozenEdges = returnReduce(nFrozenEdges, sumOp());
-
- Info<< "Frozen edges "
- << returnReduce(nFrozenEdges, sumOp())
- << " out of " << returnReduce(mesh.nEdges(), sumOp())
- << endl;
-
- if (nNewFrozenEdges == nOldFrozenEdges)
- {
- break;
- }
-
- } while (true);
-
-
- if (meshChanged)
- {
- // Apply changes to current mesh
- autoPtr mapPtr = savedMeshMod.changeMesh(mesh, false);
- const mapPolyMesh& map = mapPtr();
-
- // Update fields
- mesh.updateMesh(map);
- if (map.hasMotionPoints())
- {
- mesh.movePoints(map.preMotionPoints());
- }
-
-
// Write resulting mesh
if (!overwrite)
{
@@ -1462,7 +133,6 @@ int main(int argc, char *argv[])
mesh.write();
}
-
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
index adb8bdcd5b..a7db08ba71 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
+++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
@@ -55,6 +55,7 @@ Description
#include "edgeCollapser.H"
#include "meshTools.H"
#include "Pair.H"
+#include "globalIndex.H"
using namespace Foam;
@@ -569,26 +570,47 @@ int main(int argc, char *argv[])
// Mesh change engine
edgeCollapser cutter(mesh);
- pointField newPoints(mesh.points());
+ const edgeList& edges = mesh.edges();
+ const pointField& points = mesh.points();
+
+ pointField newPoints(points);
+
+ PackedBoolList collapseEdge(mesh.nEdges());
+ Map collapsePointToLocation(mesh.nPoints());
// Get new positions and construct collapse network
forAllConstIter(Map, edgeToPos, iter)
{
label edgeI = iter.key();
- const edge& e = mesh.edges()[edgeI];
+ const edge& e = edges[edgeI];
+
+ collapseEdge[edgeI] = true;
+ collapsePointToLocation.set(e[1], points[e[0]]);
- cutter.collapseEdge(edgeI, e[0]);
newPoints[e[0]] = iter();
}
// Move master point to destination.
mesh.movePoints(newPoints);
+ List allPointInfo;
+ const globalIndex globalPoints(mesh.nPoints());
+ labelList pointPriority(mesh.nPoints(), 0);
+
+ cutter.consistentCollapse
+ (
+ globalPoints,
+ pointPriority,
+ collapsePointToLocation,
+ collapseEdge,
+ allPointInfo
+ );
+
// Topo change container
polyTopoChange meshMod(mesh);
// Insert
- cutter.setRefinement(meshMod);
+ cutter.setRefinement(allPointInfo, meshMod);
// Do changes
autoPtr morphMap = meshMod.changeMesh(mesh, false);
diff --git a/applications/utilities/mesh/generation/cvMesh/Allwmake b/applications/utilities/mesh/generation/cvMesh/Allwmake
index eb453138ce..d88e8cee4a 100755
--- a/applications/utilities/mesh/generation/cvMesh/Allwmake
+++ b/applications/utilities/mesh/generation/cvMesh/Allwmake
@@ -4,7 +4,8 @@ set -x
wmake libso conformalVoronoiMesh
wmake
-wmake cvMeshBackgroundMesh
+#wmake cvMeshBackgroundMesh
(cd cvMeshSurfaceSimplify && ./Allwmake)
+wmake cellSizeAndAlignmentGrid
# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
index 27241b2965..9461fa3725 100644
--- a/applications/utilities/mesh/generation/cvMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/Make/options
@@ -2,7 +2,7 @@ EXE_DEBUG = -DFULLDEBUG -g -O0
EXE_FROUNDING_MATH = -frounding-math
EXE_NDEBUG = -DNDEBUG
-CGAL_EXACT =
+CGAL_EXACT = /*-DCGAL_DONT_USE_LAZY_KERNEL*/
CGAL_INEXACT = -DCGAL_INEXACT
include $(GENERAL_RULES)/CGAL
@@ -19,7 +19,9 @@ EXE_INC = \
-I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
- -I$(LIB_SRC)/triSurface/lnInclude
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -IvectorTools
EXE_LIBS = \
$(CGAL_LIBS) \
@@ -32,4 +34,5 @@ EXE_LIBS = \
-ledgeMesh \
-lfileFormats \
-ltriSurface \
- -ldynamicMesh
+ -ldynamicMesh \
+ -lsampling
diff --git a/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/files b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/files
new file mode 100644
index 0000000000..83b77fdc77
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/files
@@ -0,0 +1,2 @@
+cellSizeAndAlignmentGrid.C
+EXE = $(FOAM_USER_APPBIN)/cellSizeAndAlignmentGrid
diff --git a/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/options b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/options
new file mode 100644
index 0000000000..31d0d80858
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/options
@@ -0,0 +1,40 @@
+EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+
+CGAL_EXACT = /*-DCGAL_DONT_USE_LAZY_KERNEL*/
+CGAL_INEXACT = -DCGAL_INEXACT
+
+include $(GENERAL_RULES)/CGAL
+
+
+EXE_INC = \
+ ${EXE_FROUNDING_MATH} \
+ ${EXE_NDEBUG} \
+ ${CGAL_INEXACT} \
+ ${CGAL_INC} \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+ -I$(LIB_SRC)/edgeMesh/lnInclude \
+ -I$(HOME)/OpenFOAM/OpenFOAM-dev/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
+ -I$(HOME)/OpenFOAM/OpenFOAM-dev/applications/utilities/mesh/generation/cvMesh/vectorTools
+
+EXE_LIBS = \
+ $(CGAL_LIBS) \
+ -lmpfr \
+ -lboost_thread \
+ -lconformalVoronoiMesh \
+ -lfiniteVolume \
+ -lmeshTools \
+ -ldecompositionMethods \
+ -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
+ -ledgeMesh \
+ -ltriSurface \
+ -ldynamicMesh \
+ -lsampling \
+ -lfileFormats
diff --git a/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C
new file mode 100644
index 0000000000..27cc47f823
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C
@@ -0,0 +1,711 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+Foam::tmp > filterFarPoints
+(
+ const Triangulation& mesh,
+ const Field& field
+)
+{
+ tmp > tNewField(new Field(field.size()));
+ Field& newField = tNewField();
+
+ label added = 0;
+ label count = 0;
+
+ for
+ (
+ typename Triangulation::Finite_vertices_iterator vit =
+ mesh.finite_vertices_begin();
+ vit != mesh.finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->real())
+ {
+ newField[added++] = field[count];
+ }
+
+ count++;
+ }
+
+ newField.resize(added);
+
+ return tNewField;
+}
+
+
+template
+autoPtr buildMap
+(
+ const T& mesh,
+ labelListList& pointPoints
+)
+{
+ pointPoints.setSize(mesh.vertexCount());
+
+ globalIndex globalIndexing(mesh.vertexCount());
+
+ for
+ (
+ typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
+ vit != mesh.finite_vertices_end();
+ ++vit
+ )
+ {
+ if (!vit->real())
+ {
+ continue;
+ }
+
+ std::list adjVerts;
+ mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
+
+ DynamicList indices(adjVerts.size());
+
+ for
+ (
+ typename std::list::const_iterator
+ adjVertI = adjVerts.begin();
+ adjVertI != adjVerts.end();
+ ++adjVertI
+ )
+ {
+ typename T::Vertex_handle vh = *adjVertI;
+
+ if (!vh->farPoint())
+ {
+ indices.append
+ (
+ globalIndexing.toGlobal(vh->procIndex(), vh->index())
+ );
+ }
+ }
+
+ pointPoints[vit->index()].transfer(indices);
+ }
+
+ List > compactMap;
+
+ return autoPtr
+ (
+ new mapDistribute
+ (
+ globalIndexing,
+ pointPoints,
+ compactMap
+ )
+ );
+}
+
+
+template
+Foam::tmp buildAlignmentField(const T& mesh)
+{
+ tmp tAlignments
+ (
+ new triadField(mesh.vertexCount(), triad::unset)
+ );
+ triadField& alignments = tAlignments();
+
+ for
+ (
+ typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
+ vit != mesh.finite_vertices_end();
+ ++vit
+ )
+ {
+ if (!vit->real())
+ {
+ continue;
+ }
+
+ alignments[vit->index()] = triad
+ (
+ vit->alignment().x(),
+ vit->alignment().y(),
+ vit->alignment().z()
+ );
+ }
+
+ return tAlignments;
+}
+
+
+template
+Foam::tmp buildPointField(const T& mesh)
+{
+ tmp tPoints
+ (
+ new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
+ );
+ pointField& points = tPoints();
+
+ for
+ (
+ typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
+ vit != mesh.finite_vertices_end();
+ ++vit
+ )
+ {
+ if (!vit->real())
+ {
+ continue;
+ }
+
+ points[vit->index()] = topoint(vit->point());
+ }
+
+ return tPoints;
+}
+
+
+int main(int argc, char *argv[])
+{
+ #include "setRootCase.H"
+ #include "createTime.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ label maxRefinementIterations = 0;
+ label maxSmoothingIterations = 200;
+ scalar minResidual = 0;
+ scalar defaultCellSize = 0.0004;
+ scalar nearFeatDistSqrCoeff = 1e-8;
+
+
+ // Need to decouple vertex and cell type from this class?
+ // Vertex must have:
+ // + index
+ // + procIndex
+ // - type should be optional
+ cellShapeControlMesh mesh(runTime);
+
+ IOdictionary cvMeshDict
+ (
+ IOobject
+ (
+ "cvMeshDict",
+ runTime.system(),
+ runTime,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+ );
+
+ Random rndGen(64293*Pstream::myProcNo());
+
+ searchableSurfaces allGeometry
+ (
+ IOobject
+ (
+ "cvSearchableSurfaces",
+ runTime.constant(),
+ "triSurface",
+ runTime,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ ),
+ cvMeshDict.subDict("geometry")
+ );
+
+ conformationSurfaces geometryToConformTo
+ (
+ runTime,
+ rndGen,
+ allGeometry,
+ cvMeshDict.subDict("surfaceConformation")
+ );
+
+ autoPtr bMesh;
+ if (Pstream::parRun())
+ {
+ bMesh.set
+ (
+ new backgroundMeshDecomposition
+ (
+ runTime,
+ rndGen,
+ geometryToConformTo,
+ cvMeshDict.subDict("backgroundMeshDecomposition")
+ )
+ );
+ }
+
+ // Nice to have IO for the delaunay mesh
+ // IO depend on vertex type.
+ //
+ // Define a delaunay mesh as:
+ // + list of points of the triangulation
+ // + optionally a list of cells
+
+ Info<< nl << "Loop over surfaces" << endl;
+
+ forAll(geometryToConformTo.surfaces(), sI)
+ {
+ const label surfI = geometryToConformTo.surfaces()[sI];
+
+ const searchableSurface& surface =
+ geometryToConformTo.geometry()[surfI];
+
+ Info<< nl << "Inserting points from surface " << surface.name()
+ << " (" << surface.type() << ")" << endl;
+
+ const tmp tpoints = surface.points();
+ const pointField& points = tpoints();
+
+ Info<< " Number of points = " << points.size() << endl;
+
+ forAll(points, pI)
+ {
+ // Is the point in the extendedFeatureEdgeMesh? If so get the
+ // point normal, otherwise get the surface normal from
+ // searchableSurface
+
+ pointIndexHit info;
+ label infoFeature;
+ geometryToConformTo.findFeaturePointNearest
+ (
+ points[pI],
+ nearFeatDistSqrCoeff,
+ info,
+ infoFeature
+ );
+
+
+ autoPtr pointAlignment;
+
+ if (info.hit())
+ {
+ const extendedFeatureEdgeMesh& features =
+ geometryToConformTo.features()[infoFeature];
+
+ vectorField norms = features.featurePointNormals(info.index());
+
+ // Create a triad from these norms.
+ pointAlignment.set(new triad());
+ forAll(norms, nI)
+ {
+ pointAlignment() += norms[nI];
+ }
+
+ pointAlignment().normalize();
+ pointAlignment().orthogonalize();
+ }
+ else
+ {
+ geometryToConformTo.findEdgeNearest
+ (
+ points[pI],
+ nearFeatDistSqrCoeff,
+ info,
+ infoFeature
+ );
+
+ if (info.hit())
+ {
+ const extendedFeatureEdgeMesh& features =
+ geometryToConformTo.features()[infoFeature];
+
+ vectorField norms = features.edgeNormals(info.index());
+
+ // Create a triad from these norms.
+ pointAlignment.set(new triad());
+ forAll(norms, nI)
+ {
+ pointAlignment() += norms[nI];
+ }
+
+ pointAlignment().normalize();
+ pointAlignment().orthogonalize();
+ }
+ else
+ {
+ pointField ptField(1, points[pI]);
+ scalarField distField(1, nearFeatDistSqrCoeff);
+ List infoList(1, pointIndexHit());
+
+ surface.findNearest(ptField, distField, infoList);
+
+ vectorField normals(1);
+ surface.getNormal(infoList, normals);
+
+ pointAlignment.set(new triad(normals[0]));
+ }
+ }
+
+ if (Pstream::parRun())
+ {
+ if (bMesh().positionOnThisProcessor(points[pI]))
+ {
+ CellSizeDelaunay::Vertex_handle vh = mesh.insert
+ (
+ points[pI],
+ defaultCellSize,
+ pointAlignment()
+ );
+ }
+ }
+ else
+ {
+ CellSizeDelaunay::Vertex_handle vh = mesh.insert
+ (
+ points[pI],
+ defaultCellSize,
+ pointAlignment()
+ );
+ }
+ }
+ }
+
+
+ for (label iter = 0; iter < maxRefinementIterations; ++iter)
+ {
+ DynamicList ptsToInsert;
+
+ for
+ (
+ CellSizeDelaunay::Finite_cells_iterator cit =
+ mesh.finite_cells_begin();
+ cit != mesh.finite_cells_end();
+ ++cit
+ )
+ {
+ const point newPoint =
+ topoint
+ (
+ CGAL::centroid
+ (
+ cit->vertex(0)->point(),
+ cit->vertex(1)->point(),
+ cit->vertex(2)->point(),
+ cit->vertex(3)->point()
+ )
+ );
+
+ if (geometryToConformTo.inside(newPoint))
+ {
+ ptsToInsert.append(newPoint);
+ }
+ }
+
+ Info<< " Adding " << returnReduce(ptsToInsert.size(), sumOp())
+ << endl;
+
+ forAll(ptsToInsert, ptI)
+ {
+ mesh.insert
+ (
+ ptsToInsert[ptI],
+ defaultCellSize,
+ triad::unset
+ );
+ }
+ }
+
+
+ if (Pstream::parRun())
+ {
+ mesh.distribute(bMesh);
+ }
+
+ labelListList pointPoints;
+ autoPtr meshDistributor = buildMap(mesh, pointPoints);
+
+ triadField alignments = buildAlignmentField(mesh);
+ pointField points = buildPointField(mesh);
+
+ mesh.printInfo(Info);
+
+
+ // Setup the sizes and alignments on each point
+ triadField fixedAlignments(mesh.vertexCount(), triad::unset);
+
+ for
+ (
+ CellSizeDelaunay::Finite_vertices_iterator vit =
+ mesh.finite_vertices_begin();
+ vit != mesh.finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->real())
+ {
+ const tensor& alignment = vit->alignment();
+
+ fixedAlignments[vit->index()] = triad
+ (
+ alignment.x(),
+ alignment.y(),
+ alignment.z()
+ );
+ }
+ }
+
+
+ Info<< nl << "Smoothing alignments" << endl;
+
+ for (label iter = 0; iter < maxSmoothingIterations; iter++)
+ {
+ Info<< "Iteration " << iter;
+
+ meshDistributor().distribute(points);
+ meshDistributor().distribute(alignments);
+
+ scalar residual = 0;
+
+ triadField triadAv(alignments.size(), triad::unset);
+
+ forAll(pointPoints, pI)
+ {
+ const labelList& pPoints = pointPoints[pI];
+
+ if (pPoints.empty())
+ {
+ continue;
+ }
+
+ const triad& oldTriad = alignments[pI];
+ triad& newTriad = triadAv[pI];
+
+ forAll(pPoints, adjPointI)
+ {
+ const label adjPointIndex = pPoints[adjPointI];
+
+ scalar dist = mag(points[pI] - points[adjPointIndex]);
+
+ dist = max(dist, SMALL);
+
+ triad tmpTriad = alignments[adjPointIndex];
+
+ for (direction dir = 0; dir < 3; dir++)
+ {
+ if (tmpTriad.set(dir))
+ {
+ tmpTriad[dir] *= (1.0/dist);
+ }
+ }
+
+ newTriad += tmpTriad;
+ }
+
+ newTriad.normalize();
+ newTriad.orthogonalize();
+ newTriad = newTriad.sortxyz();
+
+ // Enforce the boundary conditions
+ const triad& fixedAlignment = fixedAlignments[pI];
+
+ label nFixed = 0;
+
+ forAll(fixedAlignment, dirI)
+ {
+ if (fixedAlignment[dirI] != triad::unset[dirI])
+ {
+ nFixed++;
+ }
+ }
+
+ if (nFixed == 1)
+ {
+ forAll(fixedAlignment, dirI)
+ {
+ if (fixedAlignment.set(dirI))
+ {
+ newTriad.align(fixedAlignment[dirI]);
+ }
+ }
+ }
+ else if (nFixed == 2)
+ {
+ forAll(fixedAlignment, dirI)
+ {
+ if (fixedAlignment.set(dirI))
+ {
+ newTriad[dirI] = fixedAlignment[dirI];
+ }
+ else
+ {
+ newTriad[dirI] = triad::unset[dirI];
+ }
+ }
+
+ newTriad.orthogonalize();
+ }
+ else if (nFixed == 3)
+ {
+ forAll(fixedAlignment, dirI)
+ {
+ if (fixedAlignment.set(dirI))
+ {
+ newTriad[dirI] = fixedAlignment[dirI];
+ }
+ }
+ }
+
+ if (newTriad.set(vector::X) && oldTriad.set(vector::X))
+ {
+ scalar dotProd = (oldTriad.x() & newTriad.x());
+
+ scalar diff = mag(dotProd) - 1.0;
+ residual += mag(diff);
+ }
+ if (newTriad.set(vector::Y) && oldTriad.set(vector::Y))
+ {
+ scalar dotProd = (oldTriad.y() & newTriad.y());
+
+ scalar diff = mag(dotProd) - 1.0;
+ residual += mag(diff);
+ }
+ if (newTriad.set(vector::Z) && oldTriad.set(vector::Z))
+ {
+ scalar dotProd = (oldTriad.z() & newTriad.z());
+
+ scalar diff = mag(dotProd) - 1.0;
+ residual += mag(diff);
+ }
+ }
+
+ forAll(alignments, pI)
+ {
+ alignments[pI] = triadAv[pI].sortxyz();
+ }
+
+ reduce(residual, sumOp());
+
+ Info<< ", Residual = " << residual << endl;
+
+ if (residual <= minResidual)
+ {
+ break;
+ }
+ }
+
+
+ // Write alignments to a .obj file
+ OFstream str(runTime.path()/"alignments.obj");
+
+ forAll(alignments, pI)
+ {
+ const triad& tri = alignments[pI];
+
+ if (tri.set())
+ {
+ forAll(tri, dirI)
+ {
+ meshTools::writeOBJ(str, points[pI], tri[dirI] + points[pI]);
+ }
+ }
+ }
+
+
+ // Remove the far points
+ pointIOField pointsIO
+ (
+ IOobject
+ (
+ "points",
+ runTime.constant(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ filterFarPoints(mesh, points)
+ );
+
+ scalarField sizes(points.size(), defaultCellSize);
+ scalarIOField sizesIO
+ (
+ IOobject
+ (
+ "sizes",
+ runTime.constant(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ filterFarPoints(mesh, sizes)
+ );
+
+ triadIOField alignmentsIO
+ (
+ IOobject
+ (
+ "alignments",
+ runTime.constant(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ filterFarPoints(mesh, alignments)
+ );
+
+ pointsIO.write();
+ sizesIO.write();
+ alignmentsIO.write();
+
+ Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+
+ Info<< "\nEnd\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/changesLaurence b/applications/utilities/mesh/generation/cvMesh/changesLaurence
deleted file mode 100644
index 5f55eff462..0000000000
--- a/applications/utilities/mesh/generation/cvMesh/changesLaurence
+++ /dev/null
@@ -1,5 +0,0 @@
-wmake/rules/General/CGAL
-
- -lboost_thread
- -lboost_thread-mt
-
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/meshQualityDict b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/meshQualityDict
new file mode 100644
index 0000000000..fa5319e087
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/meshQualityDict
@@ -0,0 +1,73 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+ version 2.0;
+ format ascii;
+
+ root "";
+ case "";
+ instance "";
+ local "";
+
+ class dictionary;
+ object meshQualityDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Maximum non-orthogonality allowed. Set to 180 to disable.
+maxNonOrtho 65;
+
+//- Max skewness allowed. Set to <0 to disable.
+maxBoundarySkewness 50;
+maxInternalSkewness 10;
+
+//- Max concaveness allowed. Is angle (in degrees) below which concavity
+// is allowed. 0 is straight face, <0 would be convex face.
+// Set to 180 to disable.
+maxConcave 80;
+
+//- Minimum quality of the tet formed by the face-centre
+// and variable base point minimum decomposition triangles and
+// the cell centre. This has to be a positive number for tracking
+// to work. Set to very negative number (e.g. -1E30) to
+// disable.
+// <0 = inside out tet,
+// 0 = flat tet
+// 1 = regular tet
+minTetQuality 1e-30;
+
+//- Minimum pyramid volume. Is absolute volume of cell pyramid.
+// Set to a sensible fraction of the smallest cell volume expected.
+// Set to very negative number (e.g. -1E30) to disable.
+minVol 1e-20;
+
+//- Minimum face area. Set to <0 to disable.
+minArea -1;
+
+//- Minimum face twist. Set to <-1 to disable. dot product of face normal
+//- and face centre triangles normal
+minTwist 0.001;
+
+//- minimum normalised cell determinant
+//- 1 = hex, <= 0 = folded or flattened illegal cell
+minDeterminant 0.001;
+
+//- minFaceWeight (0 -> 0.5)
+minFaceWeight 0.02;
+
+//- minVolRatio (0 -> 1)
+minVolRatio 0.01;
+
+//must be >0 for Fluent compatibility
+minTriangleTwist -1;
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
new file mode 100644
index 0000000000..a0327ce04f
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+Foam::DelaunayMesh::DelaunayMesh()
+:
+ Triangulation(),
+ vertexCount_(0),
+ cellCount_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+template
+Foam::DelaunayMesh::~DelaunayMesh()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+template
+void Foam::DelaunayMesh::reset()
+{
+ Info<< "Clearing triangulation" << endl;
+
+ this->clear();
+
+ resetVertexCount();
+ resetCellCount();
+}
+
+
+template
+void Foam::DelaunayMesh::insertPoints(const List& vertices)
+{
+ rangeInsertWithInfo
+ (
+ vertices.begin(),
+ vertices.end(),
+ true
+ );
+}
+
+
+template
+bool Foam::DelaunayMesh::Traits_for_spatial_sort::Less_x_3::
+operator()
+(
+ const Point_3& p,
+ const Point_3& q
+) const
+{
+ return typename Gt::Less_x_3()(*(p.first), *(q.first));
+}
+
+template
+bool Foam::DelaunayMesh::Traits_for_spatial_sort::Less_y_3::
+operator()
+(
+ const Point_3& p,
+ const Point_3& q
+) const
+{
+ return typename Gt::Less_y_3()(*(p.first), *(q.first));
+}
+
+template
+bool Foam::DelaunayMesh::Traits_for_spatial_sort::Less_z_3::
+operator()
+(
+ const Point_3& p,
+ const Point_3& q
+) const
+{
+ return typename Gt::Less_z_3()(*(p.first), *(q.first));
+}
+
+template
+typename Foam::DelaunayMesh::Traits_for_spatial_sort::Less_x_3
+Foam::DelaunayMesh::Traits_for_spatial_sort::less_x_3_object()
+const
+{
+ return Less_x_3();
+}
+
+template
+typename Foam::DelaunayMesh::Traits_for_spatial_sort::Less_y_3
+Foam::DelaunayMesh::Traits_for_spatial_sort::less_y_3_object()
+const
+{
+ return Less_y_3();
+}
+
+template
+typename Foam::DelaunayMesh::Traits_for_spatial_sort::Less_z_3
+Foam::DelaunayMesh::Traits_for_spatial_sort::less_z_3_object()
+const
+{
+ return Less_z_3();
+}
+
+
+template
+template
+void Foam::DelaunayMesh::rangeInsertWithInfo
+(
+ PointIterator begin,
+ PointIterator end,
+ bool printErrors
+)
+{
+ typedef DynamicList
+ <
+ std::pair
+ <
+ const typename Triangulation::Point*,
+ label
+ >
+ > vectorPairPointIndex;
+
+ vectorPairPointIndex points;
+
+ label count = 0;
+ for (PointIterator it = begin; it != end; ++it)
+ {
+ points.append
+ (
+ std::make_pair(&(it->point()), count++)
+ );
+ }
+
+ std::random_shuffle(points.begin(), points.end());
+
+ spatial_sort
+ (
+ points.begin(),
+ points.end(),
+ Traits_for_spatial_sort()
+ );
+
+ Vertex_handle hint;
+
+ for
+ (
+ typename vectorPairPointIndex::const_iterator p = points.begin();
+ p != points.end();
+ ++p
+ )
+ {
+ const size_t checkInsertion = Triangulation::number_of_vertices();
+
+ hint = this->insert(*(p->first), hint);
+
+ const Vb& vert = *(begin + p->second);
+
+ if (checkInsertion != Triangulation::number_of_vertices() - 1)
+ {
+ if (printErrors)
+ {
+ Vertex_handle nearV =
+ Triangulation::nearest_vertex(*(p->first));
+
+ Pout<< "Failed insertion : " << vert.info()
+ << " nearest : " << nearV->info();
+ }
+ }
+ else
+ {
+ hint->index() = getNewVertexIndex();
+ hint->type() = vert.type();
+ hint->procIndex() = vert.procIndex();
+ hint->targetCellSize() = vert.targetCellSize();
+ hint->alignment() = vert.alignment();
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "DelaunayMeshIO.C"
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
new file mode 100644
index 0000000000..bdeee880e7
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+class DelaunayMesh
+:
+ public Triangulation
+{
+public:
+
+ typedef typename Triangulation::Cell_handle Cell_handle;
+ typedef typename Triangulation::Vertex_handle Vertex_handle;
+ typedef typename Triangulation::Point Point;
+ typedef typename Triangulation::Facet Facet;
+
+ typedef typename Triangulation::Finite_vertices_iterator
+ Finite_vertices_iterator;
+ typedef typename Triangulation::Finite_cells_iterator
+ Finite_cells_iterator;
+ typedef typename Triangulation::Finite_facets_iterator
+ Finite_facets_iterator;
+
+ typedef HashSet
+ <
+ Pair,
+ FixedList::Hash<>
+ > labelPairHashSet;
+
+
+private:
+
+ // Private data
+
+ //- Keep track of the number of vertices that have been added.
+ // This allows a unique index to be assigned to each vertex.
+ mutable label vertexCount_;
+
+ //- Keep track of the number of cells that have been added.
+ // This allows a unique index to be assigned to each cell.
+ mutable label cellCount_;
+
+ //- Spatial sort traits to use with a pair of point pointers and an int.
+ // Taken from a post on the CGAL lists: 2010-01/msg00004.html by
+ // Sebastien Loriot (Geometry Factory).
+ struct Traits_for_spatial_sort
+ :
+ public Triangulation::Geom_traits
+ {
+ typedef typename Triangulation::Geom_traits Gt;
+
+ typedef std::pair
+ Point_3;
+
+ struct Less_x_3
+ {
+ bool operator()(const Point_3& p, const Point_3& q) const;
+ };
+
+ struct Less_y_3
+ {
+ bool operator()(const Point_3& p, const Point_3& q) const;
+ };
+
+ struct Less_z_3
+ {
+ bool operator()(const Point_3& p, const Point_3& q) const;
+ };
+
+ Less_x_3 less_x_3_object() const;
+ Less_y_3 less_y_3_object() const;
+ Less_z_3 less_z_3_object() const;
+ };
+
+
+ // Private Member Functions
+
+ void sortFaces
+ (
+ faceList& faces,
+ labelList& owner,
+ labelList& neighbour
+ ) const;
+
+ void addPatches
+ (
+ const label nInternalFaces,
+ faceList& faces,
+ labelList& owner,
+ labelList& patchSizes,
+ labelList& patchStarts,
+ const List >& patchFaces,
+ const List >& patchOwners
+ ) const;
+
+ //- Disallow default bitwise copy construct
+ DelaunayMesh(const DelaunayMesh&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const DelaunayMesh&);
+
+
+public:
+
+ // Constructors
+
+ //- Construct from components
+ DelaunayMesh();
+
+
+ //- Destructor
+ ~DelaunayMesh();
+
+
+ // Member Functions
+
+ inline label getNewVertexIndex() const;
+
+ inline label getNewCellIndex() const;
+
+ inline label cellCount() const;
+
+ inline void resetCellCount();
+
+ inline label vertexCount() const;
+
+ inline void resetVertexCount();
+
+
+ //- Remove the entire triangulation
+ void reset();
+
+ void insertPoints(const List& vertices);
+
+ //- Function inserting points into a triangulation and setting the
+ // index and type data of the point in the correct order. This is
+ // faster than inserting points individually.
+ //
+ // Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
+ // Sebastien Loriot (Geometry Factory).
+ template
+ void rangeInsertWithInfo
+ (
+ PointIterator begin,
+ PointIterator end,
+ bool printErrors = true
+ );
+
+
+ // Queries
+
+ void printInfo(Ostream& os) const;
+
+ //- Create an fvMesh from the triangulation.
+ // The mesh is not parallel consistent - only used for viewing
+ autoPtr createMesh
+ (
+ const fileName& name,
+ const Time& runTime,
+ labelList& vertexMap,
+ labelList& cellMap
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "DelaunayMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "DelaunayMesh.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshI.H
new file mode 100644
index 0000000000..841e5c9024
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshI.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+inline Foam::label Foam::DelaunayMesh::getNewVertexIndex() const
+{
+ label id = vertexCount_++;
+
+ if (id == labelMax)
+ {
+ WarningIn
+ (
+ "Foam::DelaunayMesh::getNewVertexIndex() const"
+ ) << "Vertex counter has overflowed." << endl;
+ }
+
+ return id;
+}
+
+
+template
+inline Foam::label Foam::DelaunayMesh::getNewCellIndex() const
+{
+ label id = cellCount_++;
+
+ if (id == labelMax)
+ {
+ WarningIn
+ (
+ "Foam::DelaunayMesh::getNewCellIndex() const"
+ ) << "Cell counter has overflowed." << endl;
+ }
+
+ return id;
+}
+
+
+template
+Foam::label Foam::DelaunayMesh::cellCount() const
+{
+ return cellCount_;
+}
+
+
+template
+void Foam::DelaunayMesh::resetCellCount()
+{
+ cellCount_ = 0;
+}
+
+
+template
+Foam::label Foam::DelaunayMesh::vertexCount() const
+{
+ return vertexCount_;
+}
+
+
+template
+void Foam::DelaunayMesh::resetVertexCount()
+{
+ vertexCount_ = 0;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshIO.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshIO.C
new file mode 100644
index 0000000000..63a9889308
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshIO.C
@@ -0,0 +1,391 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+void Foam::DelaunayMesh::sortFaces
+(
+ faceList& faces,
+ labelList& owner,
+ labelList& neighbour
+) const
+{
+ // Upper triangular order:
+ // + owner is sorted in ascending cell order
+ // + within each block of equal value for owner, neighbour is sorted in
+ // ascending cell order.
+ // + faces sorted to correspond
+ // e.g.
+ // owner | neighbour
+ // 0 | 2
+ // 0 | 23
+ // 0 | 71
+ // 1 | 23
+ // 1 | 24
+ // 1 | 91
+
+ List ownerNeighbourPair(owner.size());
+
+ forAll(ownerNeighbourPair, oNI)
+ {
+ ownerNeighbourPair[oNI] = labelPair(owner[oNI], neighbour[oNI]);
+ }
+
+ Info<< nl
+ << "Sorting faces, owner and neighbour into upper triangular order"
+ << endl;
+
+ labelList oldToNew;
+
+ sortedOrder(ownerNeighbourPair, oldToNew);
+
+ oldToNew = invert(oldToNew.size(), oldToNew);
+
+ inplaceReorder(oldToNew, faces);
+ inplaceReorder(oldToNew, owner);
+ inplaceReorder(oldToNew, neighbour);
+}
+
+
+template
+void Foam::DelaunayMesh::addPatches
+(
+ const label nInternalFaces,
+ faceList& faces,
+ labelList& owner,
+ labelList& patchSizes,
+ labelList& patchStarts,
+ const List >& patchFaces,
+ const List >& patchOwners
+) const
+{
+ label nPatches = patchFaces.size();
+
+ patchSizes.setSize(nPatches, -1);
+ patchStarts.setSize(nPatches, -1);
+
+ label nBoundaryFaces = 0;
+
+ forAll(patchFaces, p)
+ {
+ patchSizes[p] = patchFaces[p].size();
+ patchStarts[p] = nInternalFaces + nBoundaryFaces;
+
+ nBoundaryFaces += patchSizes[p];
+ }
+
+ faces.setSize(nInternalFaces + nBoundaryFaces);
+ owner.setSize(nInternalFaces + nBoundaryFaces);
+
+ label faceI = nInternalFaces;
+
+ forAll(patchFaces, p)
+ {
+ forAll(patchFaces[p], f)
+ {
+ faces[faceI] = patchFaces[p][f];
+ owner[faceI] = patchOwners[p][f];
+
+ faceI++;
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
+
+template
+void Foam::DelaunayMesh::printInfo(Ostream& os) const
+{
+ PrintTable triInfoTable("Mesh Statistics");
+
+ triInfoTable.add("Points", Triangulation::number_of_vertices());
+ triInfoTable.add("Edges", Triangulation::number_of_finite_edges());
+ triInfoTable.add("Faces", Triangulation::number_of_finite_facets());
+ triInfoTable.add("Cells", Triangulation::number_of_finite_cells());
+
+ scalar minSize = GREAT;
+ scalar maxSize = 0;
+
+ for
+ (
+ Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
+ vit != Triangulation::finite_vertices_end();
+ ++vit
+ )
+ {
+ if (!vit->farPoint())
+ {
+ minSize = min(vit->targetCellSize(), minSize);
+ maxSize = max(vit->targetCellSize(), maxSize);
+ }
+ }
+
+ Info<< incrIndent;
+ triInfoTable.print(Info, true, true);
+
+ Info<< "Size (Min/Max) = "
+ << returnReduce(minSize, minOp()) << " "
+ << returnReduce(maxSize, maxOp()) << endl;
+
+ Info<< decrIndent;
+}
+
+
+template
+Foam::autoPtr
+Foam::DelaunayMesh::createMesh
+(
+ const fileName& name,
+ const Time& runTime,
+ labelList& vertexMap,
+ labelList& cellMap
+) const
+{
+ pointField points(Triangulation::number_of_vertices());
+ faceList faces(Triangulation::number_of_finite_facets());
+ labelList owner(Triangulation::number_of_finite_facets());
+ labelList neighbour(Triangulation::number_of_finite_facets());
+
+ wordList patchNames(1, "cvMesh_defaultPatch");
+ wordList patchTypes(1, wallPolyPatch::typeName);
+
+ labelList patchSizes(1, 0);
+ labelList patchStarts(1, 0);
+
+ List > patchFaces(1, DynamicList());
+ List > patchOwners(1, DynamicList());
+
+ vertexMap.setSize(Triangulation::number_of_vertices());
+ cellMap.setSize(Triangulation::number_of_finite_cells());
+
+ // Calculate pts and a map of point index to location in pts.
+ label vertI = 0;
+
+ for
+ (
+ Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
+ vit != Triangulation::finite_vertices_end();
+ ++vit
+ )
+ {
+ if (!vit->farPoint())
+ {
+ vertexMap[vit->index()] = vertI;
+ points[vertI] = topoint(vit->point());
+ vertI++;
+ }
+ }
+
+ points.setSize(vertI);
+
+ // Index the cells
+ label cellI = 0;
+
+ for
+ (
+ Finite_cells_iterator cit = Triangulation::finite_cells_begin();
+ cit != Triangulation::finite_cells_end();
+ ++cit
+ )
+ {
+ if
+ (
+ !cit->hasFarPoint()
+ && !Triangulation::is_infinite(cit)
+ )
+ {
+ cellMap[cit->cellIndex()] = cellI++;
+ }
+ }
+
+ label faceI = 0;
+ labelList verticesOnTriFace(3, -1);
+ face newFace(verticesOnTriFace);
+
+ for
+ (
+ Finite_facets_iterator fit = Triangulation::finite_facets_begin();
+ fit != Triangulation::finite_facets_end();
+ ++fit
+ )
+ {
+ const Cell_handle c1(fit->first);
+ const int oppositeVertex = fit->second;
+ const Cell_handle c2(c1->neighbor(oppositeVertex));
+
+ label c1I = Cb::ctFar;
+ bool c1Real = false;
+ if (!c1->hasFarPoint() && !Triangulation::is_infinite(c1))
+ {
+ c1I = cellMap[c1->cellIndex()];
+ c1Real = true;
+ }
+
+ label c2I = Cb::ctFar;
+ bool c2Real = false;
+ if (!c2->hasFarPoint() && !Triangulation::is_infinite(c2))
+ {
+ c2I = cellMap[c2->cellIndex()];
+ c2Real = true;
+ }
+
+ if (!c1Real && !c2Real)
+ {
+ // Both tets are outside, skip
+ continue;
+ }
+
+ label ownerCell = -1;
+ label neighbourCell = -1;
+
+ for (label i = 0; i < 3; i++)
+ {
+ verticesOnTriFace[i] = vertexMap
+ [
+ c1->vertex
+ (
+ Triangulation::vertex_triple_index(oppositeVertex, i)
+ )->index()
+ ];
+ }
+
+ newFace = face(verticesOnTriFace);
+
+ if (!c1Real || !c2Real)
+ {
+ // Boundary face...
+ if (!c1Real)
+ {
+ //... with c1 outside
+ ownerCell = c2I;
+ }
+ else
+ {
+ // ... with c2 outside
+ ownerCell = c1I;
+
+ reverse(newFace);
+ }
+
+ patchFaces[0].append(newFace);
+ patchOwners[0].append(ownerCell);
+ }
+ else
+ {
+ // Internal face...
+ if (c1I < c2I)
+ {
+ // ...with c1 as the ownerCell
+ ownerCell = c1I;
+ neighbourCell = c2I;
+
+ reverse(newFace);
+ }
+ else
+ {
+ // ...with c2 as the ownerCell
+ ownerCell = c2I;
+ neighbourCell = c1I;
+ }
+
+ faces[faceI] = newFace;
+ owner[faceI] = ownerCell;
+ neighbour[faceI] = neighbourCell;
+ faceI++;
+ }
+ }
+
+ faces.setSize(faceI);
+ owner.setSize(faceI);
+ neighbour.setSize(faceI);
+
+ sortFaces(faces, owner, neighbour);
+
+ addPatches
+ (
+ faceI,
+ faces,
+ owner,
+ patchSizes,
+ patchStarts,
+ patchFaces,
+ patchOwners
+ );
+
+ autoPtr meshPtr
+ (
+ new fvMesh
+ (
+ IOobject
+ (
+ name,
+ runTime.timeName(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ xferMove(points),
+ xferMove(faces),
+ xferMove(owner),
+ xferMove(neighbour)
+ )
+ );
+
+ List patches(patchStarts.size());
+
+ label nValidPatches = 0;
+
+ forAll(patches, p)
+ {
+ patches[nValidPatches] = polyPatch::New
+ (
+ patchTypes[p],
+ patchNames[p],
+ patchSizes[p],
+ patchStarts[p],
+ nValidPatches,
+ meshPtr().boundaryMesh()
+ ).ptr();
+
+ nValidPatches++;
+ }
+
+ patches.setSize(nValidPatches);
+
+ meshPtr().addFvPatches(patches);
+
+ return meshPtr;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
new file mode 100644
index 0000000000..9a5351a1f4
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
@@ -0,0 +1,936 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+Foam::autoPtr
+Foam::DistributedDelaunayMesh::buildMap
+(
+ const List& toProc
+)
+{
+ // Determine send map
+ // ~~~~~~~~~~~~~~~~~~
+
+ // 1. Count
+ labelList nSend(Pstream::nProcs(), 0);
+
+ forAll(toProc, i)
+ {
+ label procI = toProc[i];
+
+ nSend[procI]++;
+ }
+
+ // Send over how many I need to receive
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ labelListList sendSizes(Pstream::nProcs());
+
+ sendSizes[Pstream::myProcNo()] = nSend;
+
+ combineReduce(sendSizes, UPstream::listEq());
+
+ // 2. Size sendMap
+ labelListList sendMap(Pstream::nProcs());
+
+ forAll(nSend, procI)
+ {
+ sendMap[procI].setSize(nSend[procI]);
+
+ nSend[procI] = 0;
+ }
+
+ // 3. Fill sendMap
+ forAll(toProc, i)
+ {
+ label procI = toProc[i];
+
+ sendMap[procI][nSend[procI]++] = i;
+ }
+
+ // Determine receive map
+ // ~~~~~~~~~~~~~~~~~~~~~
+
+ labelListList constructMap(Pstream::nProcs());
+
+ // Local transfers first
+ constructMap[Pstream::myProcNo()] = identity
+ (
+ sendMap[Pstream::myProcNo()].size()
+ );
+
+ label constructSize = constructMap[Pstream::myProcNo()].size();
+
+ forAll(constructMap, procI)
+ {
+ if (procI != Pstream::myProcNo())
+ {
+ label nRecv = sendSizes[procI][Pstream::myProcNo()];
+
+ constructMap[procI].setSize(nRecv);
+
+ for (label i = 0; i < nRecv; i++)
+ {
+ constructMap[procI][i] = constructSize++;
+ }
+ }
+ }
+
+ return autoPtr
+ (
+ new mapDistribute
+ (
+ constructSize,
+ sendMap.xfer(),
+ constructMap.xfer()
+ )
+ );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template
+Foam::DistributedDelaunayMesh::DistributedDelaunayMesh()
+:
+ DelaunayMesh(),
+ allBackgroundMeshBounds_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+template
+Foam::DistributedDelaunayMesh::~DistributedDelaunayMesh()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template
+bool Foam::DistributedDelaunayMesh::distributeBoundBoxes
+(
+ const boundBox& bb
+)
+{
+ allBackgroundMeshBounds_.reset(new List(Pstream::nProcs()));
+
+ // Give the bounds of every processor to every other processor
+ allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
+
+ Pstream::gatherList(allBackgroundMeshBounds_());
+ Pstream::scatterList(allBackgroundMeshBounds_());
+
+ return true;
+}
+
+
+template
+bool Foam::DistributedDelaunayMesh::isLocal
+(
+ const Vertex_handle& v
+) const
+{
+ return isLocal(v->procIndex());
+}
+
+
+template
+bool Foam::DistributedDelaunayMesh::isLocal
+(
+ const label localProcIndex
+) const
+{
+ return localProcIndex == Pstream::myProcNo();
+}
+
+
+template
+Foam::labelList Foam::DistributedDelaunayMesh::overlapProcessors
+(
+ const point& centre,
+ const scalar radiusSqr
+) const
+{
+ DynamicList toProc(Pstream::nProcs());
+
+ forAll(allBackgroundMeshBounds_(), procI)
+ {
+ // Test against the bounding box of the processor
+ if
+ (
+ !isLocal(procI)
+ && allBackgroundMeshBounds_()[procI].overlaps(centre, radiusSqr)
+ )
+ {
+ toProc.append(procI);
+ }
+ }
+
+ return toProc;
+}
+
+
+template
+bool Foam::DistributedDelaunayMesh::checkProcBoundaryCell
+(
+ const Cell_handle& cit,
+ Map& circumsphereOverlaps
+) const
+{
+ const Foam::point& cc = cit->dual();
+
+ const scalar crSqr = magSqr
+ (
+ cc - topoint(cit->vertex(0)->point())
+ );
+
+ labelList circumsphereOverlap = overlapProcessors
+ (
+ cc,
+ sqr(1.01)*crSqr
+ );
+
+ cit->cellIndex() = this->getNewCellIndex();
+
+ if (!circumsphereOverlap.empty())
+ {
+ circumsphereOverlaps.insert(cit->cellIndex(), circumsphereOverlap);
+
+ return true;
+ }
+
+ return false;
+}
+
+
+template
+void Foam::DistributedDelaunayMesh::findProcessorBoundaryCells
+(
+ Map& circumsphereOverlaps
+) const
+{
+ // Start by assuming that all the cells have no index
+ // If they do, they have already been visited so ignore them
+
+ labelHashSet cellToCheck
+ (
+ Triangulation::number_of_finite_cells()
+ /Pstream::nProcs()
+ );
+
+ for
+ (
+ All_cells_iterator cit = Triangulation::all_cells_begin();
+ cit != Triangulation::all_cells_end();
+ ++cit
+ )
+ {
+ if (Triangulation::is_infinite(cit))
+ {
+ // Index of infinite vertex in this cell.
+ int i = cit->index(Triangulation::infinite_vertex());
+
+ Cell_handle c = cit->neighbor(i);
+
+ if (c->unassigned())
+ {
+ c->cellIndex() = this->getNewCellIndex();
+
+ if (checkProcBoundaryCell(c, circumsphereOverlaps))
+ {
+ cellToCheck.insert(c->cellIndex());
+ }
+ }
+ }
+ else if (cit->parallelDualVertex())
+ {
+ if (cit->unassigned())
+ {
+ if (checkProcBoundaryCell(cit, circumsphereOverlaps))
+ {
+ cellToCheck.insert(cit->cellIndex());
+ }
+ }
+ }
+ }
+
+ for
+ (
+ Finite_cells_iterator cit = Triangulation::finite_cells_begin();
+ cit != Triangulation::finite_cells_end();
+ ++cit
+ )
+ {
+ if (cellToCheck.found(cit->cellIndex()))
+ {
+ // Get the neighbours and check them
+ for (label adjCellI = 0; adjCellI < 4; ++adjCellI)
+ {
+ Cell_handle citNeighbor = cit->neighbor(adjCellI);
+
+ // Ignore if has far point or previously visited
+ if
+ (
+ !citNeighbor->unassigned()
+ || !citNeighbor->internalOrBoundaryDualVertex()
+ || Triangulation::is_infinite(citNeighbor)
+ )
+ {
+ continue;
+ }
+
+ checkProcBoundaryCell
+ (
+ citNeighbor,
+ circumsphereOverlaps
+ );
+ }
+ }
+ }
+}
+
+
+template
+void Foam::DistributedDelaunayMesh::markVerticesToRefer
+(
+ const Map& circumsphereOverlaps,
+ PtrList& referralVertices,
+ DynamicList& targetProcessor,
+ DynamicList& parallelInfluenceVertices
+)
+{
+ // Relying on the order of iteration of cells being the same as before
+ for
+ (
+ Finite_cells_iterator cit = Triangulation::finite_cells_begin();
+ cit != Triangulation::finite_cells_end();
+ ++cit
+ )
+ {
+ if (Triangulation::is_infinite(cit))
+ {
+ continue;
+ }
+
+ Map::const_iterator iter =
+ circumsphereOverlaps.find(cit->cellIndex());
+
+ // Pre-tested circumsphere potential influence
+ if (iter != circumsphereOverlaps.cend())
+ {
+ const labelList& citOverlaps = iter();
+
+ forAll(citOverlaps, cOI)
+ {
+ label procI = citOverlaps[cOI];
+
+ for (int i = 0; i < 4; i++)
+ {
+ Vertex_handle v = cit->vertex(i);
+
+ if (v->farPoint())
+ {
+ continue;
+ }
+
+ label vProcIndex = v->procIndex();
+ label vIndex = v->index();
+
+ const labelPair procIndexPair(vProcIndex, vIndex);
+
+ // Using the hashSet to ensure that each vertex is only
+ // referred once to each processor.
+ // Do not refer a vertex to its own processor.
+ if (vProcIndex != procI)
+ {
+ if (referralVertices[procI].insert(procIndexPair))
+ {
+ targetProcessor.append(procI);
+
+ parallelInfluenceVertices.append
+ (
+ Vb
+ (
+ v->point(),
+ v->index(),
+ v->type(),
+ v->procIndex()
+ )
+ );
+
+ parallelInfluenceVertices.last().targetCellSize() =
+ v->targetCellSize();
+ parallelInfluenceVertices.last().alignment() =
+ v->alignment();
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+
+template
+Foam::label Foam::DistributedDelaunayMesh::referVertices
+(
+ const DynamicList& targetProcessor,
+ DynamicList& parallelVertices,
+ PtrList& referralVertices,
+ labelPairHashSet& receivedVertices
+)
+{
+ DynamicList referredVertices(targetProcessor.size());
+
+ const label preDistributionSize = parallelVertices.size();
+
+ mapDistribute pointMap = buildMap(targetProcessor);
+
+ // Make a copy of the original list.
+ DynamicList originalParallelVertices(parallelVertices);
+
+ pointMap.distribute(parallelVertices);
+
+ for (label procI = 0; procI < Pstream::nProcs(); procI++)
+ {
+ const labelList& constructMap = pointMap.constructMap()[procI];
+
+ if (constructMap.size())
+ {
+ forAll(constructMap, i)
+ {
+ const Vb& v = parallelVertices[constructMap[i]];
+
+ if
+ (
+ v.procIndex() != Pstream::myProcNo()
+ && !receivedVertices.found(labelPair(v.procIndex(), v.index()))
+ )
+ {
+ referredVertices.append(v);
+
+ receivedVertices.insert
+ (
+ labelPair(v.procIndex(), v.index())
+ );
+ }
+ }
+ }
+ }
+
+ label preInsertionSize = Triangulation::number_of_vertices();
+
+ labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
+ (
+ referredVertices.begin(),
+ referredVertices.end()
+ );
+
+ if (!pointsNotInserted.empty())
+ {
+ for
+ (
+ typename labelPairHashSet::const_iterator iter
+ = pointsNotInserted.begin();
+ iter != pointsNotInserted.end();
+ ++iter
+ )
+ {
+ if (receivedVertices.found(iter.key()))
+ {
+ receivedVertices.erase(iter.key());
+ }
+ }
+ }
+
+ boolList pointInserted(parallelVertices.size(), true);
+
+ forAll(parallelVertices, vI)
+ {
+ const labelPair procIndexI
+ (
+ parallelVertices[vI].procIndex(),
+ parallelVertices[vI].index()
+ );
+
+ if (pointsNotInserted.found(procIndexI))
+ {
+ pointInserted[vI] = false;
+ }
+ }
+
+ pointMap.reverseDistribute(preDistributionSize, pointInserted);
+
+ forAll(originalParallelVertices, vI)
+ {
+ const label procIndex = targetProcessor[vI];
+
+ if (!pointInserted[vI])
+ {
+ if (referralVertices[procIndex].size())
+ {
+ if
+ (
+ !referralVertices[procIndex].unset
+ (
+ labelPair
+ (
+ originalParallelVertices[vI].procIndex(),
+ originalParallelVertices[vI].index()
+ )
+ )
+ )
+ {
+ Pout<< "*** not found "
+ << originalParallelVertices[vI].procIndex()
+ << " " << originalParallelVertices[vI].index() << endl;
+ }
+
+ }
+ }
+ }
+
+ label postInsertionSize = Triangulation::number_of_vertices();
+
+ reduce(preInsertionSize, sumOp());
+ reduce(postInsertionSize, sumOp());
+
+ label nTotalToInsert = referredVertices.size();
+
+ reduce(nTotalToInsert, sumOp());
+
+ if (preInsertionSize + nTotalToInsert != postInsertionSize)
+ {
+ label nNotInserted =
+ returnReduce(pointsNotInserted.size(), sumOp());
+
+ Info<< " Inserted = "
+ << setw(name(label(Triangulation::number_of_finite_cells())).size())
+ << nTotalToInsert - nNotInserted
+ << " / " << nTotalToInsert << endl;
+
+ nTotalToInsert -= nNotInserted;
+ }
+ else
+ {
+ Info<< " Inserted = " << nTotalToInsert << endl;
+ }
+
+ return nTotalToInsert;
+}
+
+
+template
+void Foam::DistributedDelaunayMesh::sync
+(
+ const boundBox& bb,
+ PtrList& referralVertices,
+ labelPairHashSet& receivedVertices,
+ bool iterateReferral
+)
+{
+ if (!Pstream::parRun())
+ {
+ return;
+ }
+
+ if (allBackgroundMeshBounds_.empty())
+ {
+ distributeBoundBoxes(bb);
+ }
+
+ label nVerts = Triangulation::number_of_vertices();
+ label nCells = Triangulation::number_of_finite_cells();
+
+ DynamicList parallelInfluenceVertices(0.1*nVerts);
+ DynamicList targetProcessor(0.1*nVerts);
+
+ // Some of these values will not be used, i.e. for non-real cells
+ DynamicList circumcentre(0.1*nVerts);
+ DynamicList circumradiusSqr(0.1*nVerts);
+
+ Map circumsphereOverlaps(nCells);
+
+ findProcessorBoundaryCells(circumsphereOverlaps);
+
+ Info<< " Influences = "
+ << setw(name(nCells).size())
+ << returnReduce(circumsphereOverlaps.size(), sumOp()) << " / "
+ << returnReduce(nCells, sumOp());
+
+ markVerticesToRefer
+ (
+ circumsphereOverlaps,
+ referralVertices,
+ targetProcessor,
+ parallelInfluenceVertices
+ );
+
+ referVertices
+ (
+ targetProcessor,
+ parallelInfluenceVertices,
+ referralVertices,
+ receivedVertices
+ );
+
+ if (iterateReferral)
+ {
+ label oldNReferred = 0;
+ label nIterations = 1;
+
+ Info<< incrIndent << indent
+ << "Iteratively referring referred vertices..."
+ << endl;
+ do
+ {
+ Info<< indent << "Iteration " << nIterations++ << ":";
+
+ circumsphereOverlaps.clear();
+ targetProcessor.clear();
+ parallelInfluenceVertices.clear();
+
+ findProcessorBoundaryCells(circumsphereOverlaps);
+
+ nCells = Triangulation::number_of_finite_cells();
+
+ Info<< " Influences = "
+ << setw(name(nCells).size())
+ << returnReduce(circumsphereOverlaps.size(), sumOp())
+ << " / "
+ << returnReduce(nCells, sumOp());
+
+ markVerticesToRefer
+ (
+ circumsphereOverlaps,
+ referralVertices,
+ targetProcessor,
+ parallelInfluenceVertices
+ );
+
+ label nReferred = referVertices
+ (
+ targetProcessor,
+ parallelInfluenceVertices,
+ referralVertices,
+ receivedVertices
+ );
+
+ if (nReferred == 0 || nReferred == oldNReferred)
+ {
+ break;
+ }
+
+ oldNReferred = nReferred;
+
+ } while (true);
+
+ Info<< decrIndent;
+ }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+template
+bool Foam::DistributedDelaunayMesh::distribute
+(
+ const boundBox& bb
+)
+{
+ notImplemented
+ (
+ "Foam::DistributedDelaunayMesh::distribute"
+ "("
+ " const boundBox& bb"
+ ")"
+ );
+
+ if (!Pstream::parRun())
+ {
+ return false;
+ }
+
+ distributeBoundBoxes(bb);
+
+ return true;
+}
+
+
+template
+Foam::autoPtr
+Foam::DistributedDelaunayMesh::distribute
+(
+ const backgroundMeshDecomposition& decomposition
+)
+{
+ if (!Pstream::parRun())
+ {
+ return autoPtr();
+ }
+
+ distributeBoundBoxes(decomposition.procBounds());
+
+ DynamicList points(Triangulation::number_of_vertices());
+
+ for
+ (
+ Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
+ vit != Triangulation::finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->real())
+ {
+ points.append(topoint(vit->point()));
+ }
+ }
+
+ autoPtr mapDist = decomposition.distributePoints(points);
+
+ return mapDist;
+}
+
+
+template
+void Foam::DistributedDelaunayMesh::sync(const boundBox& bb)
+{
+ if (!Pstream::parRun())
+ {
+ return;
+ }
+
+ if (allBackgroundMeshBounds_.empty())
+ {
+ distributeBoundBoxes(bb);
+ }
+
+ const label nApproxReferred =
+ Triangulation::number_of_vertices()
+ /Pstream::nProcs();
+
+ PtrList referralVertices(Pstream::nProcs());
+ forAll(referralVertices, procI)
+ {
+ if (!isLocal(procI))
+ {
+ referralVertices.set(procI, new labelPairHashSet(nApproxReferred));
+ }
+ }
+
+ labelPairHashSet receivedVertices(nApproxReferred);
+
+ sync
+ (
+ bb,
+ referralVertices,
+ receivedVertices,
+ true
+ );
+}
+
+
+template
+template
+typename Foam::DistributedDelaunayMesh::labelPairHashSet
+Foam::DistributedDelaunayMesh::rangeInsertReferredWithInfo
+(
+ PointIterator begin,
+ PointIterator end,
+ bool printErrors
+)
+{
+ const boundBox& bb = allBackgroundMeshBounds_()[Pstream::myProcNo()];
+
+ typedef DynamicList
+ <
+ std::pair
+ > vectorPairPointIndex;
+
+ vectorPairPointIndex pointsBbDistSqr;
+
+ label count = 0;
+ for (PointIterator it = begin; it != end; ++it)
+ {
+ const pointFromPoint samplePoint = topoint(it->point());
+
+ if (!bb.contains(samplePoint))
+ {
+ const Foam::point nearestPoint = bb.nearest(samplePoint);
+
+ const scalar distFromBbSqr = magSqr(nearestPoint - samplePoint);
+
+ pointsBbDistSqr.append
+ (
+ std::make_pair(distFromBbSqr, count++)
+ );
+ }
+ }
+
+ std::random_shuffle(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
+
+ // Sort in ascending order by the distance of the point from the centre
+ // of the processor bounding box
+ sort(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
+
+ typename Triangulation::Vertex_handle hint;
+
+ typename Triangulation::Locate_type lt;
+ int li, lj;
+
+ label nNotInserted = 0;
+
+ labelPairHashSet uninserted
+ (
+ Triangulation::number_of_vertices()
+ /Pstream::nProcs()
+ );
+
+ for
+ (
+ typename vectorPairPointIndex::const_iterator p =
+ pointsBbDistSqr.begin();
+ p != pointsBbDistSqr.end();
+ ++p
+ )
+ {
+ const size_t checkInsertion = Triangulation::number_of_vertices();
+
+ const Vb& vert = *(begin + p->second);
+ const Point& pointToInsert = vert.point();
+
+ // Locate the point
+ Cell_handle c = Triangulation::locate(pointToInsert, lt, li, lj, hint);
+
+ if (lt == Triangulation::VERTEX)
+ {
+ if (printErrors)
+ {
+ Vertex_handle nearV =
+ Triangulation::nearest_vertex(pointToInsert);
+
+ Pout<< "Failed insertion, point already exists" << nl
+ << "Failed insertion : " << vert.info()
+ << " nearest : " << nearV->info();
+ }
+
+ uninserted.insert(labelPair(vert.procIndex(), vert.index()));
+ nNotInserted++;
+
+ continue;
+ }
+
+ // Get the cells that conflict with p in a vector V,
+ // and a facet on the boundary of this hole in f.
+ std::vector V;
+ typename Triangulation::Facet f;
+
+ Triangulation::find_conflicts
+ (
+ pointToInsert,
+ c,
+ CGAL::Oneset_iterator(f),
+ std::back_inserter(V)
+ );
+
+ bool insert = false;
+ for (size_t i = 0; i < V.size(); ++i)
+ {
+ if (V[i]->real() || V[i]->hasFarPoint())
+ {
+ insert = true;
+ break;
+ }
+ }
+
+ if (insert)
+ {
+ hint = Triangulation::insert_in_hole
+ (
+ pointToInsert,
+ V.begin(),
+ V.end(),
+ f.first,
+ f.second
+ );
+
+ if (checkInsertion != Triangulation::number_of_vertices() - 1)
+ {
+ if (printErrors)
+ {
+ Vertex_handle nearV =
+ Triangulation::nearest_vertex(pointToInsert);
+
+ Pout<< "Failed insertion : " << vert.info()
+ << " nearest : " << nearV->info();
+ }
+ }
+ else
+ {
+ hint->index() = vert.index();
+ hint->type() = vert.type();
+ hint->procIndex() = vert.procIndex();
+ hint->targetCellSize() = vert.targetCellSize();
+ hint->alignment() = vert.alignment();
+ }
+ }
+ else
+ {
+ uninserted.insert(labelPair(vert.procIndex(), vert.index()));
+ nNotInserted++;
+ }
+ }
+
+ return uninserted;
+}
+
+
+// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.H
new file mode 100644
index 0000000000..124f4b39dc
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.H
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+class DistributedDelaunayMesh
+:
+ public DelaunayMesh
+{
+public:
+
+ typedef typename Triangulation::Vertex_handle Vertex_handle;
+ typedef typename Triangulation::Cell_handle Cell_handle;
+ typedef typename Triangulation::Point Point;
+
+ typedef typename Triangulation::Finite_vertices_iterator
+ Finite_vertices_iterator;
+ typedef typename Triangulation::Finite_cells_iterator
+ Finite_cells_iterator;
+ typedef typename Triangulation::All_cells_iterator
+ All_cells_iterator;
+
+ typedef typename DelaunayMesh::labelPairHashSet
+ labelPairHashSet;
+
+
+private:
+
+ autoPtr > allBackgroundMeshBounds_;
+
+
+ // Private Member Functions
+
+ //-
+ bool distributeBoundBoxes(const boundBox& bb);
+
+ //-
+ bool isLocal(const Vertex_handle& v) const;
+
+ bool isLocal(const label localProcIndex) const;
+
+ labelList overlapProcessors
+ (
+ const point& centre,
+ const scalar radiusSqr
+ ) const;
+
+ bool checkProcBoundaryCell
+ (
+ const Cell_handle& cit,
+ Map& circumsphereOverlaps
+ ) const;
+
+ void findProcessorBoundaryCells
+ (
+ Map& circumsphereOverlaps
+ ) const;
+
+ void markVerticesToRefer
+ (
+ const Map& circumsphereOverlaps,
+ PtrList& referralVertices,
+ DynamicList& targetProcessor,
+ DynamicList& parallelInfluenceVertices
+ );
+
+ label referVertices
+ (
+ const DynamicList& targetProcessor,
+ DynamicList& parallelVertices,
+ PtrList& referralVertices,
+ labelPairHashSet& receivedVertices
+ );
+
+ //- Disallow default bitwise copy construct
+ DistributedDelaunayMesh(const DistributedDelaunayMesh