diff --git a/applications/test/vectorTools/Make/files b/applications/test/vectorTools/Make/files
new file mode 100644
index 0000000000..0b30b98f8f
--- /dev/null
+++ b/applications/test/vectorTools/Make/files
@@ -0,0 +1,3 @@
+Test-vectorTools.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-vectorTools
diff --git a/applications/test/vectorTools/Make/options b/applications/test/vectorTools/Make/options
new file mode 100644
index 0000000000..9e015e6078
--- /dev/null
+++ b/applications/test/vectorTools/Make/options
@@ -0,0 +1 @@
+EXE_INC = -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/vectorTools
diff --git a/applications/test/vectorTools/Test-vectorTools.C b/applications/test/vectorTools/Test-vectorTools.C
new file mode 100644
index 0000000000..85d18ed989
--- /dev/null
+++ b/applications/test/vectorTools/Test-vectorTools.C
@@ -0,0 +1,71 @@
+#include "vector.H"
+#include "IOstreams.H"
+#include "vectorTools.H"
+#include "unitConversion.H"
+
+using namespace Foam;
+
+
+void test(const vector& a, const vector& b, const scalar tolerance)
+{
+ Info<< "Vectors " << a << " and " << b
+ << " are (to tolerance of " << tolerance << "): ";
+
+ if (vectorTools::areParallel(a, b, tolerance))
+ Info<< " parallel ";
+
+ if (vectorTools::areOrthogonal(a, b, tolerance))
+ Info<< " orthogonal ";
+
+ if (vectorTools::areAcute(a, b))
+ Info<< " acute ";
+
+ if (vectorTools::areObtuse(a, b))
+ Info<< " obtuse ";
+
+ Info<< ", angle = " << vectorTools::degAngleBetween(a, b);
+
+ Info<< endl;
+}
+
+
+int main()
+{
+ vector a(1.0, 1.0, 1.0);
+ vector b(2.0, 2.0, 2.0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ a = vector(1,0,0);
+ b = vector(0,2,0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ a = vector(1,0,0);
+ b = vector(-1,0,0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ a = vector(1,0,0);
+ b = vector(-1,2,0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ return 0;
+}
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
index 170a0a890d..75f04740ef 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseDict
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
@@ -2,16 +2,24 @@
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
-| \\ / A nd | Web: www.OpenFOAM.org |
+| \\ / A nd | Web: http://www.openfoam.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
+
FoamFile
{
- version 2.0;
- format ascii;
- class dictionary;
- object collapseDict;
+ version 2.0;
+ format ascii;
+
+ root "";
+ case "";
+ instance "";
+ local "";
+
+ class dictionary;
+ object collapseDict;
}
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// If on, after collapsing check the quality of the mesh. If bad faces are
diff --git a/applications/utilities/mesh/generation/Allwmake b/applications/utilities/mesh/generation/Allwmake
index d557fa34f8..48ca530c96 100755
--- a/applications/utilities/mesh/generation/Allwmake
+++ b/applications/utilities/mesh/generation/Allwmake
@@ -9,5 +9,10 @@ extrude2DMesh/Allwmake
wmake snappyHexMesh
+if [ -d "$CGAL_ARCH_PATH" ]
+then
+ cvMesh/Allwmake
+ cv2DMesh/Allwmake
+fi
# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwclean b/applications/utilities/mesh/generation/cv2DMesh/Allwclean
new file mode 100755
index 0000000000..d0ae53e415
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/Allwclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+set -x
+
+wclean libso conformalVoronoi2DMesh
+wclean
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwmake b/applications/utilities/mesh/generation/cv2DMesh/Allwmake
new file mode 100755
index 0000000000..5486849957
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/Allwmake
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+set -x
+
+wmake libso conformalVoronoi2DMesh
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2DKernel.H b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2DKernel.H
new file mode 100644
index 0000000000..28002f962f
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2DKernel.H
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 Vb;
+typedef CGAL::indexedFace Fb;
+
+#ifdef CGAL_HIERARCHY
+
+ // Data structures for hierarchical Delaunay triangulation which is more
+ // efficient but also uses more storage
+ #include "CGAL/Triangulation_hierarchy_2.h"
+
+ typedef CGAL::Triangulation_hierarchy_vertex_base_2 Vbh;
+ typedef CGAL::Triangulation_data_structure_2 Tds;
+ typedef CGAL::Delaunay_triangulation_2 DT;
+ typedef CGAL::Triangulation_hierarchy_2 Delaunay;
+
+#else
+
+ // Data structures for standard Delaunay triangulation
+ typedef CGAL::Triangulation_data_structure_2 Tds;
+ typedef CGAL::Delaunay_triangulation_2 Delaunay;
+
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C
new file mode 100644
index 0000000000..8a5931df61
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C
@@ -0,0 +1,997 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 index() = nVert++;
+ }
+ }
+ }
+
+ Info<< nVert << " vertices inserted" << endl;
+
+ if (meshControls().objOutput())
+ {
+ // Checking validity of triangulation
+ assert(is_valid());
+
+ writeTriangles("initial_triangles.obj", true);
+ writeFaces("initial_faces.obj", true);
+ }
+}
+
+
+void Foam::CV2D::insertSurfacePointPairs()
+{
+ startOfSurfacePointPairs_ = number_of_vertices();
+
+ if (meshControls().insertSurfaceNearestPointPairs())
+ {
+ insertSurfaceNearestPointPairs();
+ }
+
+ write("nearest");
+
+ // Insertion of point-pairs for near-points may cause protrusions
+ // so insertBoundaryConformPointPairs must be executed last
+ if (meshControls().insertSurfaceNearPointPairs())
+ {
+ insertSurfaceNearPointPairs();
+ }
+
+ startOfBoundaryConformPointPairs_ = number_of_vertices();
+}
+
+
+void Foam::CV2D::boundaryConform()
+{
+ if (!meshControls().insertSurfaceNearestPointPairs())
+ {
+ markNearBoundaryPoints();
+ }
+
+ // Mark all the faces as SAVE_CHANGED
+ for
+ (
+ Triangulation::Finite_faces_iterator fit = finite_faces_begin();
+ fit != finite_faces_end();
+ fit++
+ )
+ {
+ fit->faceIndex() = Fb::SAVE_CHANGED;
+ }
+
+ for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++)
+ {
+ label nIntersections = insertBoundaryConformPointPairs
+ (
+ "surfaceIntersections_" + Foam::name(iter) + ".obj"
+ );
+
+ if (nIntersections == 0)
+ {
+ break;
+ }
+ else
+ {
+ Info<< "BC iteration " << iter << ": "
+ << nIntersections << " point-pairs inserted" << endl;
+ }
+
+ // Any faces changed by insertBoundaryConformPointPairs will now
+ // be marked CHANGED, mark those as SAVE_CHANGED and those that
+ // remained SAVE_CHANGED as UNCHANGED
+ for
+ (
+ Triangulation::Finite_faces_iterator fit = finite_faces_begin();
+ fit != finite_faces_end();
+ fit++
+ )
+ {
+ if (fit->faceIndex() == Fb::SAVE_CHANGED)
+ {
+ fit->faceIndex() = Fb::UNCHANGED;
+ }
+ else if (fit->faceIndex() == Fb::CHANGED)
+ {
+ fit->faceIndex() = Fb::SAVE_CHANGED;
+ }
+ }
+ }
+
+ Info<< nl;
+
+ write("boundary");
+}
+
+
+void Foam::CV2D::removeSurfacePointPairs()
+{
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->index() >= startOfSurfacePointPairs_)
+ {
+ remove(vit);
+ }
+ }
+}
+
+
+void Foam::CV2D::newPoints()
+{
+ const scalar relaxation = relaxationModel_->relaxation();
+
+ Info<< "Relaxation = " << relaxation << endl;
+
+ Field dualVertices(number_of_faces());
+
+ label dualVerti = 0;
+
+ // Find the dual point of each tetrahedron and assign it an index.
+ for
+ (
+ Triangulation::Finite_faces_iterator fit = finite_faces_begin();
+ fit != finite_faces_end();
+ ++fit
+ )
+ {
+ fit->faceIndex() = -1;
+
+ if
+ (
+ fit->vertex(0)->internalOrBoundaryPoint()
+ || fit->vertex(1)->internalOrBoundaryPoint()
+ || fit->vertex(2)->internalOrBoundaryPoint()
+ )
+ {
+ fit->faceIndex() = dualVerti;
+
+ dualVertices[dualVerti] = toPoint2D(circumcenter(fit));
+
+ dualVerti++;
+ }
+ }
+
+ dualVertices.setSize(dualVerti);
+
+ Field displacementAccumulator
+ (
+ startOfSurfacePointPairs_,
+ vector2D::zero
+ );
+
+ // Calculate target size and alignment for vertices
+ scalarField sizes
+ (
+ number_of_vertices(),
+ meshControls().minCellSize()
+ );
+
+ Field alignments
+ (
+ number_of_vertices(),
+ vector2D(1, 0)
+ );
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalOrBoundaryPoint())
+ {
+ point2D vert = toPoint2D(vit->point());
+
+ // alignment and size determination
+ pointIndexHit pHit;
+ label hitSurface = -1;
+
+ qSurf_.findSurfaceNearest
+ (
+ toPoint3D(vert),
+ meshControls().span2(),
+ pHit,
+ hitSurface
+ );
+
+ if (pHit.hit())
+ {
+ vectorField norm(1);
+ allGeometry_[hitSurface].getNormal
+ (
+ List(1, pHit),
+ norm
+ );
+
+ alignments[vit->index()] = toPoint2D(norm[0]);
+
+ sizes[vit->index()] =
+ cellSizeControl_.cellSize(toPoint3D(vit->point()));
+ }
+ }
+ }
+
+ // Info<< "Calculated alignments" << endl;
+
+ scalar cosAlignmentAcceptanceAngle = 0.68;
+
+ // Upper and lower edge length ratios for weight
+ scalar u = 1.0;
+ scalar l = 0.7;
+
+ PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true);
+
+ std::list pointsToInsert;
+
+ for
+ (
+ Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+ eit != finite_edges_end();
+ eit++
+ )
+ {
+ Vertex_handle vA = eit->first->vertex(cw(eit->second));
+ Vertex_handle vB = eit->first->vertex(ccw(eit->second));
+
+ if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint())
+ {
+ continue;
+ }
+
+ const point2D& dualV1 = dualVertices[eit->first->faceIndex()];
+ const point2D& dualV2 =
+ dualVertices[eit->first->neighbor(eit->second)->faceIndex()];
+
+ scalar dualEdgeLength = mag(dualV1 - dualV2);
+
+ point2D dVA = toPoint2D(vA->point());
+ point2D dVB = toPoint2D(vB->point());
+
+ Field alignmentDirsA(2);
+
+ alignmentDirsA[0] = alignments[vA->index()];
+ alignmentDirsA[1] = vector2D
+ (
+ -alignmentDirsA[0].y(),
+ alignmentDirsA[0].x()
+ );
+
+ Field alignmentDirsB(2);
+
+ alignmentDirsB[0] = alignments[vB->index()];
+ alignmentDirsB[1] = vector2D
+ (
+ -alignmentDirsB[0].y(),
+ alignmentDirsB[0].x()
+ );
+
+ Field alignmentDirs(2);
+
+ forAll(alignmentDirsA, aA)
+ {
+ const vector2D& a(alignmentDirsA[aA]);
+
+ scalar maxDotProduct = 0.0;
+
+ forAll(alignmentDirsB, aB)
+ {
+ const vector2D& b(alignmentDirsB[aB]);
+
+ scalar dotProduct = a & b;
+
+ if (mag(dotProduct) > maxDotProduct)
+ {
+ maxDotProduct = mag(dotProduct);
+
+ alignmentDirs[aA] = a + sign(dotProduct)*b;
+
+ alignmentDirs[aA] /= mag(alignmentDirs[aA]);
+ }
+ }
+ }
+
+ vector2D rAB = dVA - dVB;
+
+ scalar rABMag = mag(rAB);
+
+ forAll(alignmentDirs, aD)
+ {
+ vector2D& alignmentDir = alignmentDirs[aD];
+
+ if ((rAB & alignmentDir) < 0)
+ {
+ // swap the direction of the alignment so that has the
+ // same sense as rAB
+ alignmentDir *= -1;
+ }
+
+ scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir);
+
+ if (alignmentDotProd > cosAlignmentAcceptanceAngle)
+ {
+ scalar targetFaceSize =
+ 0.5*(sizes[vA->index()] + sizes[vB->index()]);
+
+ // Test for changing aspect ratio on second alignment (first
+ // alignment is neartest surface normal)
+ // if (aD == 1)
+ // {
+ // targetFaceSize *= 2.0;
+ // }
+
+ alignmentDir *= 0.5*targetFaceSize;
+
+ vector2D delta = alignmentDir - 0.5*rAB;
+
+ if (dualEdgeLength < 0.7*targetFaceSize)
+ {
+ delta *= 0;
+ }
+ else if (dualEdgeLength < targetFaceSize)
+ {
+ delta *=
+ (
+ dualEdgeLength
+ /(targetFaceSize*(u - l))
+ - 1/((u/l) - 1)
+ );
+ }
+
+ if
+ (
+ vA->internalPoint()
+ && vB->internalPoint()
+ && rABMag > 1.75*targetFaceSize
+ && dualEdgeLength > 0.05*targetFaceSize
+ && alignmentDotProd > 0.93
+ )
+ {
+ // Point insertion
+ pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
+ }
+ else if
+ (
+ (vA->internalPoint() || vB->internalPoint())
+ && rABMag < 0.65*targetFaceSize
+ )
+ {
+ // Point removal
+
+ // Only insert a point at the midpoint of the short edge
+ // if neither attached point has already been identified
+ // to be removed.
+ if
+ (
+ pointToBeRetained[vA->index()] == true
+ && pointToBeRetained[vB->index()] == true
+ )
+ {
+ pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
+ }
+
+ if (vA->internalPoint())
+ {
+ pointToBeRetained[vA->index()] = false;
+ }
+
+ if (vB->internalPoint())
+ {
+ pointToBeRetained[vB->index()] = false;
+ }
+ }
+ else
+ {
+ if (vA->internalPoint())
+ {
+ displacementAccumulator[vA->index()] += delta;
+ }
+
+ if (vB->internalPoint())
+ {
+ displacementAccumulator[vB->index()] += -delta;
+ }
+ }
+ }
+ }
+ }
+
+ vector2D totalDisp = sum(displacementAccumulator);
+ scalar totalDist = sum(mag(displacementAccumulator));
+
+ // Relax the calculated displacement
+ displacementAccumulator *= relaxation;
+
+ label numberOfNewPoints = pointsToInsert.size();
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalPoint())
+ {
+ if (pointToBeRetained[vit->index()])
+ {
+ pointsToInsert.push_front
+ (
+ toPoint
+ (
+ toPoint2D(vit->point())
+ + displacementAccumulator[vit->index()]
+ )
+ );
+ }
+ }
+ }
+
+ // Clear the triangulation and reinsert the bounding box and feature points.
+ // This is faster than removing and moving points.
+ this->clear();
+
+ insertBoundingBox();
+
+ reinsertFeaturePoints();
+
+ startOfInternalPoints_ = number_of_vertices();
+
+ label nVert = startOfInternalPoints_;
+
+ Info<< "Inserting " << numberOfNewPoints << " new points" << endl;
+
+ // Use the range insert as it is faster than individually inserting points.
+ insert(pointsToInsert.begin(), pointsToInsert.end());
+
+ for
+ (
+ Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if
+ (
+ vit->type() == Vb::INTERNAL_POINT
+ && vit->index() == Vb::INTERNAL_POINT
+ )
+ {
+ vit->index() = nVert++;
+ }
+ }
+
+ Info<< " Total displacement = " << totalDisp << nl
+ << " Total distance = " << totalDist << nl
+ << " Points added = " << pointsToInsert.size()
+ << endl;
+
+ write("internal");
+
+ insertSurfacePointPairs();
+
+ boundaryConform();
+
+
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// Old Method
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+// for
+// (
+// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// ++vit
+// )
+// {
+// if (vit->internalPoint())
+// {
+// // Current dual-cell defining vertex ("centre")
+// point2DFromPoint defVert0 = toPoint2D(vit->point());
+
+// Triangulation::Edge_circulator ec = incident_edges(vit);
+// Triangulation::Edge_circulator ecStart = ec;
+
+// // Circulate around the edges to find the first which is not
+// // infinite
+// do
+// {
+// if (!is_infinite(ec)) break;
+// } while (++ec != ecStart);
+
+// // Store the start-end of the first non-infinte edge
+// point2D de0 = toPoint2D(circumcenter(ec->first));
+
+// // Keep track of the maximum edge length^2
+// scalar maxEdgeLen2 = 0.0;
+
+// // Keep track of the index of the longest edge
+// label edgecd0i = -1;
+
+// // Edge counter
+// label edgei = 0;
+
+// do
+// {
+// if (!is_infinite(ec))
+// {
+// // Get the end of the current edge
+// point2D de1 = toPoint2D
+// (
+// circumcenter(ec->first->neighbor(ec->second))
+// );
+
+// // Store the current edge vector
+// edges[edgei] = de1 - de0;
+
+// // Store the edge mid-point in the vertices array
+// vertices[edgei] = 0.5*(de1 + de0);
+
+// // Move the current edge end into the edge start for the
+// // next iteration
+// de0 = de1;
+
+// // Keep track of the longest edge
+
+// scalar edgeLen2 = magSqr(edges[edgei]);
+
+// if (edgeLen2 > maxEdgeLen2)
+// {
+// maxEdgeLen2 = edgeLen2;
+// edgecd0i = edgei;
+// }
+
+// edgei++;
+// }
+// } while (++ec != ecStart);
+
+// // Initialise cd0 such that the mesh will align
+// // in in the x-y directions
+// vector2D cd0(1, 0);
+
+// if (meshControls().relaxOrientation())
+// {
+// // Get the longest edge from the array and use as the primary
+// // direction of the coordinate system of the "square" cell
+// cd0 = edges[edgecd0i];
+// }
+
+// if (meshControls().nearWallAlignedDist() > 0)
+// {
+// pointIndexHit pHit = qSurf_.tree().findNearest
+// (
+// toPoint3D(defVert0),
+// meshControls().nearWallAlignedDist2()
+// );
+
+// if (pHit.hit())
+// {
+// cd0 = toPoint2D(faceNormals[pHit.index()]);
+// }
+// }
+
+// // Rotate by 45deg needed to create an averaging procedure which
+// // encourages the cells to be square
+// cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x());
+
+// // Normalise the primary coordinate direction
+// cd0 /= mag(cd0);
+
+// // Calculate the orthogonal coordinate direction
+// vector2D cd1(-cd0.y(), cd0.x());
+
+
+// // Restart the circulator
+// ec = ecStart;
+
+// // ... and the counter
+// edgei = 0;
+
+// // Initialise the displacement for the centre and sum-weights
+// vector2D disp = vector2D::zero;
+// scalar sumw = 0;
+
+// do
+// {
+// if (!is_infinite(ec))
+// {
+// // Pick up the current edge
+// const vector2D& ei = edges[edgei];
+
+// // Calculate the centre to edge-centre vector
+// vector2D deltai = vertices[edgei] - defVert0;
+
+// // Set the weight for this edge contribution
+// scalar w = 1;
+
+// if (meshControls().squares())
+// {
+// w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x());
+// // alternative weights
+// //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x());
+// //w = magSqr(ei)*mag(deltai);
+
+// // Use the following for an ~square mesh
+// // Find the coordinate contributions for this edge delta
+// scalar cd0deltai = cd0 & deltai;
+// scalar cd1deltai = cd1 & deltai;
+
+// // Create a "square" displacement
+// if (mag(cd0deltai) > mag(cd1deltai))
+// {
+// disp += (w*cd0deltai)*cd0;
+// }
+// else
+// {
+// disp += (w*cd1deltai)*cd1;
+// }
+// }
+// else
+// {
+// // Use this for a hexagon/pentagon mesh
+// disp += w*deltai;
+// }
+
+// // Sum the weights
+// sumw += w;
+// }
+// else
+// {
+// FatalErrorIn("CV2D::newPoints() const")
+// << "Infinite triangle found in internal mesh"
+// << exit(FatalError);
+// }
+
+// edgei++;
+
+// } while (++ec != ecStart);
+
+// // Calculate the average displacement
+// disp /= sumw;
+// totalDisp += disp;
+// totalDist += mag(disp);
+
+// // Move the point by a fraction of the average displacement
+// movePoint(vit, defVert0 + relaxation*disp);
+// }
+// }
+
+// Info << "\nTotal displacement = " << totalDisp
+// << " total distance = " << totalDist << endl;
+}
+
+
+//void Foam::CV2D::moveInternalPoints(const point2DField& newPoints)
+//{
+// label pointI = 0;
+
+// for
+// (
+// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// ++vit
+// )
+// {
+// if (vit->internalPoint())
+// {
+// movePoint(vit, newPoints[pointI++]);
+// }
+// }
+//}
+
+
+void Foam::CV2D::write() const
+{
+ if (meshControls().objOutput())
+ {
+ writeFaces("allFaces.obj", false);
+ writeFaces("faces.obj", true);
+ writeTriangles("allTriangles.obj", false);
+ writeTriangles("triangles.obj", true);
+ writePatch("patch.pch");
+ }
+}
+
+
+void Foam::CV2D::write(const word& stage) const
+{
+ if (meshControls().objOutput())
+ {
+ Foam::mkDir(stage + "Faces");
+ Foam::mkDir(stage + "Triangles");
+
+ writeFaces
+ (
+ stage
+ + "Faces/allFaces_"
+ + runTime_.timeName()
+ + ".obj",
+ false
+ );
+
+ writeTriangles
+ (
+ stage
+ + "Triangles/allTriangles_"
+ + runTime_.timeName()
+ + ".obj",
+ false
+ );
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H
new file mode 100644
index 0000000000..4948196fc2
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H
@@ -0,0 +1,474 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 relaxationModel_;
+
+ //- z-level
+ scalar z_;
+
+ //- Keep track of the start of the internal points
+ label startOfInternalPoints_;
+
+ //- Keep track of the start of the surface point-pairs
+ label startOfSurfacePointPairs_;
+
+ //- Keep track of the boundary conform point-pairs
+ // stored after the insertion of the surface point-pairs in case
+ // the boundary conform function is called more than once without
+ // removing and insertin the surface point-pairs
+ label startOfBoundaryConformPointPairs_;
+
+ //- Store the feature points
+ std::list featurePoints_;
+
+ //- Temporary storage for a dual-cell
+ static const label maxNvert = 20;
+ mutable point2D vertices[maxNvert+1];
+ mutable vector2D edges[maxNvert+1];
+
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ CV2D(const CV2D&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const CV2D&);
+
+
+ //- Insert point and return it's index
+ inline label insertPoint
+ (
+ const point2D& pt,
+ const label type
+ );
+
+ //- Insert point and return it's index
+ inline label insertPoint
+ (
+ const point2D& pt,
+ const label index,
+ const label type
+ );
+
+ inline label insertPoint
+ (
+ const Point& p,
+ const label index,
+ const label type
+ );
+
+ inline bool insertMirrorPoint
+ (
+ const point2D& nearSurfPt,
+ const point2D& surfPt
+ );
+
+ //- Insert a point-pair at a distance ppDist either side of
+ // surface point point surfPt in the direction n
+ inline void insertPointPair
+ (
+ const scalar mirrorDist,
+ const point2D& surfPt,
+ const vector2D& n
+ );
+
+ //- Create the initial mesh from the bounding-box
+ void insertBoundingBox();
+
+ //- Check if a point is within a line.
+ bool on2DLine(const point2D& p, const linePointRef& line);
+
+ //- Insert point groups at the feature points.
+ void insertFeaturePoints();
+
+ //- Re-insert point groups at the feature points.
+ void reinsertFeaturePoints();
+
+ //- Insert point-pairs at the given set of points using the surface
+ // normals corresponding to the given set of surface triangles
+ // and write the inserted point locations to the given file.
+ void insertPointPairs
+ (
+ const DynamicList& nearSurfacePoints,
+ const DynamicList& surfacePoints,
+ const DynamicList& surfaceTris,
+ const DynamicList& surfaceHits,
+ const fileName fName
+ );
+
+ //- Check to see if dual cell specified by given vertex iterator
+ // intersects the boundary and hence reqires a point-pair.
+ bool dualCellSurfaceIntersection
+ (
+ const Triangulation::Finite_vertices_iterator& vit
+ ) const;
+
+ //- Insert point-pairs at the nearest points on the surface to the
+ // control vertex of dual-cells which intersect the boundary in order
+ // to provide a boundary-layer mesh.
+ // NB: This is not guaranteed to close the boundary
+ void insertSurfaceNearestPointPairs();
+
+ //- Insert point-pairs at small dual-cell edges on the surface in order
+ // to improve the boundary-layer mesh generated by
+ // insertSurfaceNearestPointPairs.
+ void insertSurfaceNearPointPairs();
+
+ //- Insert point-pair and correcting the Finite_vertices_iterator
+ // to account for the additional vertices
+ void insertPointPair
+ (
+ Triangulation::Finite_vertices_iterator& vit,
+ const point2D& p,
+ const label trii,
+ const label hitSurface
+ );
+
+ //- Insert point-pair at the best intersection point between the lines
+ // from the dual-cell real centroid and it's vertices and the surface.
+ bool insertPointPairAtIntersection
+ (
+ Triangulation::Finite_vertices_iterator& vit,
+ const point2D& defVert,
+ const point2D vertices[],
+ const scalar maxProtSize
+ );
+
+ //- Insert point-pairs corresponding to dual-cells which intersect
+ // the boundary surface
+ label insertBoundaryConformPointPairs(const fileName& fName);
+
+ void markNearBoundaryPoints();
+
+ //- Restore the Delaunay contraint
+ void fast_restore_Delaunay(Vertex_handle vh);
+
+ // Flip operations used by fast_restore_Delaunay
+ void external_flip(Face_handle& f, int i);
+ bool internal_flip(Face_handle& f, int i);
+
+ //- Write all the faces and all the triangles at a particular stage.
+ void write(const word& stage) const;
+
+
+public:
+
+ //- Runtime type information
+ ClassName("CV2D");
+
+
+ // Constructors
+
+ //- Construct for given surface
+ CV2D(const Time& runTime, const dictionary& controlDict);
+
+
+ //- Destructor
+ ~CV2D();
+
+
+ // Member Functions
+
+ // Access
+
+ inline const cv2DControls& meshControls() const;
+
+
+ // Conversion functions between point2D, point and Point
+
+ inline const point2D& toPoint2D(const point&) const;
+ inline const point2DField toPoint2D(const pointField&) const;
+ inline point toPoint3D(const point2D&) const;
+
+ #ifdef CGAL_INEXACT
+ typedef const point2D& point2DFromPoint;
+ typedef const Point& PointFromPoint2D;
+ #else
+ typedef point2D point2DFromPoint;
+ typedef Point PointFromPoint2D;
+ #endif
+
+ inline point2DFromPoint toPoint2D(const Point&) const;
+ inline PointFromPoint2D toPoint(const point2D&) const;
+ inline point toPoint3D(const Point&) const;
+
+
+ // Point insertion
+
+ //- Create the initial mesh from the given internal points.
+ // Points must be inside the boundary by at least nearness
+ // otherwise they are ignored.
+ void insertPoints
+ (
+ const point2DField& points,
+ const scalar nearness
+ );
+
+ //- Create the initial mesh from the internal points in the given
+ // file. Points outside the geometry are ignored.
+ void insertPoints(const fileName& pointFileName);
+
+ //- Create the initial mesh as a regular grid of points.
+ // Points outside the geometry are ignored.
+ void insertGrid();
+
+ //- Insert all surface point-pairs from
+ // insertSurfaceNearestPointPairs and
+ // findIntersectionForOutsideCentroid
+ void insertSurfacePointPairs();
+
+ //- Insert point-pairs where there are protrusions into
+ // or out of the surface
+ void boundaryConform();
+
+
+ // Point removal
+
+ //- Remove the point-pairs introduced by insertSurfacePointPairs
+ // and boundaryConform
+ void removeSurfacePointPairs();
+
+
+ // Point motion
+
+ inline void movePoint(const Vertex_handle& vh, const Point& P);
+
+ //- Move the internal points to the given new locations and update
+ // the triangulation to ensure it is Delaunay
+ // void moveInternalPoints(const point2DField& newPoints);
+
+ //- Calculate the displacements to create the new points
+ void newPoints();
+
+ //- Extract patch names and sizes.
+ void extractPatches
+ (
+ wordList& patchNames,
+ labelList& patchSizes,
+ EdgeMap& mapEdgesRegion,
+ EdgeMap& indirectPatchEdge
+ ) const;
+
+
+ // Write
+
+ //- Write internal points to .obj file
+ void writePoints(const fileName& fName, bool internalOnly) const;
+
+ //- Write triangles as .obj file
+ void writeTriangles(const fileName& fName, bool internalOnly) const;
+
+ //- Write dual faces as .obj file
+ void writeFaces(const fileName& fName, bool internalOnly) const;
+
+ //- Calculates dual points (circumcentres of tets) and faces
+ // (point-cell walk of tets).
+ // Returns:
+ // - dualPoints (in triangle ordering)
+ // - dualFaces (compacted)
+ void calcDual
+ (
+ point2DField& dualPoints,
+ faceList& dualFaces,
+ wordList& patchNames,
+ labelList& patchSizes,
+ EdgeMap& mapEdgesRegion,
+ EdgeMap& indirectPatchEdge
+ ) const;
+
+ //- Write patch
+ void writePatch(const fileName& fName) const;
+
+ void write() const;
+};
+
+
+inline bool boundaryTriangle(const CV2D::Face_handle fc);
+inline bool outsideTriangle(const CV2D::Face_handle fc);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CV2DI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H b/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H
new file mode 100644
index 0000000000..42f4757df1
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H
@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 (p);
+}
+
+
+inline const Foam::point2DField Foam::CV2D::toPoint2D(const pointField& p) const
+{
+ point2DField temp(p.size());
+ forAll(temp, pointI)
+ {
+ temp[pointI] = point2D(p[pointI].x(), p[pointI].y());
+ }
+ return temp;
+}
+
+
+inline Foam::point Foam::CV2D::toPoint3D(const point2D& p) const
+{
+ return point(p.x(), p.y(), z_);
+}
+
+
+#ifdef CGAL_INEXACT
+
+inline Foam::CV2D::point2DFromPoint Foam::CV2D::toPoint2D(const Point& P) const
+{
+ return reinterpret_cast(P);
+}
+
+
+inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const
+{
+ return reinterpret_cast(p);
+}
+
+#else
+
+inline Foam::CV2D::point2DFromPoint Foam::CV2D::toPoint2D(const Point& P) const
+{
+ return point2D(CGAL::to_double(P.x()), CGAL::to_double(P.y()));
+}
+
+
+inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const
+{
+ return Point(p.x(), p.y());
+}
+
+#endif
+
+
+inline Foam::point Foam::CV2D::toPoint3D(const Point& P) const
+{
+ return point(CGAL::to_double(P.x()), CGAL::to_double(P.y()), z_);
+}
+
+
+inline void Foam::CV2D::movePoint(const Vertex_handle& vh, const Point& P)
+{
+ int i = vh->index();
+ int t = vh->type();
+
+ remove(vh);
+
+ Vertex_handle newVh = insert(P);
+
+ newVh->index() = i;
+ newVh->type() = t;
+
+ // label i = vh->index();
+ // move(vh, P);
+ // vh->index() = i;
+
+ //vh->set_point(P);
+ //fast_restore_Delaunay(vh);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
+
+inline bool Foam::boundaryTriangle(const CV2D::Face_handle fc)
+{
+ return boundaryTriangle
+ (
+ *fc->vertex(0),
+ *fc->vertex(1),
+ *fc->vertex(2)
+ );
+}
+
+
+inline bool Foam::outsideTriangle(const CV2D::Face_handle fc)
+{
+ return outsideTriangle
+ (
+ *fc->vertex(0),
+ *fc->vertex(1),
+ *fc->vertex(2)
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C b/applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C
new file mode 100644
index 0000000000..0223c28ebe
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C
@@ -0,0 +1,386 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 & mapEdgesRegion,
+ EdgeMap& indirectPatchEdge
+) const
+{
+ label nPatches = qSurf_.patchNames().size() + 1;
+ label defaultPatchIndex = qSurf_.patchNames().size();
+
+ patchNames.setSize(nPatches);
+ patchSizes.setSize(nPatches, 0);
+ mapEdgesRegion.clear();
+
+ const wordList& existingPatches = qSurf_.patchNames();
+
+ forAll(existingPatches, sP)
+ {
+ patchNames[sP] = existingPatches[sP];
+ }
+
+ patchNames[defaultPatchIndex] = "CV2D_default_patch";
+
+ for
+ (
+ Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+ eit != finite_edges_end();
+ ++eit
+ )
+ {
+ Face_handle fOwner = eit->first;
+ Face_handle fNeighbor = fOwner->neighbor(eit->second);
+
+ Vertex_handle vA = fOwner->vertex(cw(eit->second));
+ Vertex_handle vB = fOwner->vertex(ccw(eit->second));
+
+ if
+ (
+ (vA->internalOrBoundaryPoint() && !vB->internalOrBoundaryPoint())
+ || (vB->internalOrBoundaryPoint() && !vA->internalOrBoundaryPoint())
+ )
+ {
+ point ptA = toPoint3D(vA->point());
+ point ptB = toPoint3D(vB->point());
+
+ label patchIndex = qSurf_.findPatch(ptA, ptB);
+
+ if (patchIndex == -1)
+ {
+ patchIndex = defaultPatchIndex;
+
+ WarningIn("Foam::CV2D::extractPatches")
+ << "Dual face found that is not on a surface "
+ << "patch. Adding to CV2D_default_patch."
+ << endl;
+ }
+
+ edge e(fOwner->faceIndex(), fNeighbor->faceIndex());
+ patchSizes[patchIndex]++;
+ mapEdgesRegion.insert(e, patchIndex);
+
+ if (!pointPair(*vA, *vB))
+ {
+ indirectPatchEdge.insert(e, 1);
+ }
+ }
+ }
+}
+
+
+void Foam::CV2D::calcDual
+(
+ point2DField& dualPoints,
+ faceList& dualFaces,
+ wordList& patchNames,
+ labelList& patchSizes,
+ EdgeMap& mapEdgesRegion,
+ EdgeMap& indirectPatchEdge
+) const
+{
+ // Dual points stored in triangle order.
+ dualPoints.setSize(number_of_faces());
+ label dualVerti = 0;
+
+ for
+ (
+ Triangulation::Finite_faces_iterator fit = finite_faces_begin();
+ fit != finite_faces_end();
+ ++fit
+ )
+ {
+ if
+ (
+ fit->vertex(0)->internalOrBoundaryPoint()
+ || fit->vertex(1)->internalOrBoundaryPoint()
+ || fit->vertex(2)->internalOrBoundaryPoint()
+ )
+ {
+ fit->faceIndex() = dualVerti;
+
+ dualPoints[dualVerti++] = toPoint2D(circumcenter(fit));
+ }
+ else
+ {
+ fit->faceIndex() = -1;
+ }
+ }
+
+ dualPoints.setSize(dualVerti);
+
+ extractPatches(patchNames, patchSizes, mapEdgesRegion, indirectPatchEdge);
+
+ forAll(patchNames, patchI)
+ {
+ Info<< "Patch " << patchNames[patchI]
+ << " has size " << patchSizes[patchI] << endl;
+ }
+
+ // Create dual faces
+ // ~~~~~~~~~~~~~~~~~
+
+ dualFaces.setSize(number_of_vertices());
+ label dualFacei = 0;
+ labelList faceVerts(maxNvert);
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalOrBoundaryPoint())
+ {
+ Face_circulator fcStart = incident_faces(vit);
+ Face_circulator fc = fcStart;
+ label verti = 0;
+
+ do
+ {
+ if (!is_infinite(fc))
+ {
+ if (fc->faceIndex() < 0)
+ {
+ FatalErrorIn
+ (
+ "Foam::CV2D::calcDual"
+ "(point2DField& dualPoints, faceList& dualFaces)"
+ )<< "Dual face uses vertex defined by a triangle"
+ " defined by an external point"
+ << exit(FatalError);
+ }
+
+ // Look up the index of the triangle
+ faceVerts[verti++] = fc->faceIndex();
+ }
+ } while (++fc != fcStart);
+
+ if (faceVerts.size() > 2)
+ {
+ dualFaces[dualFacei++] =
+ face(labelList::subList(faceVerts, verti));
+ }
+ else
+ {
+ Info<< "From triangle point:" << vit->index()
+ << " coord:" << toPoint2D(vit->point())
+ << " generated illegal dualFace:" << faceVerts
+ << endl;
+ }
+ }
+ }
+
+ dualFaces.setSize(dualFacei);
+}
+
+
+void Foam::CV2D::writePatch(const fileName& fName) const
+{
+ point2DField dual2DPoints;
+ faceList dualFaces;
+ wordList patchNames;
+ labelList patchSizes;
+ EdgeMap mapEdgesRegion;
+ EdgeMap indirectPatchEdge;
+
+ calcDual
+ (
+ dual2DPoints,
+ dualFaces,
+ patchNames,
+ patchSizes,
+ mapEdgesRegion,
+ indirectPatchEdge
+ );
+
+ pointField dualPoints(dual2DPoints.size());
+ forAll(dualPoints, ip)
+ {
+ dualPoints[ip] = toPoint3D(dual2DPoints[ip]);
+ }
+
+ // Dump as primitive patch to be read by extrudeMesh.
+ OFstream str(fName);
+
+ Info<< "Writing patch to be used with extrudeMesh to file " << fName
+ << endl;
+
+ str << dualPoints << nl << dualFaces << nl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Make/files b/applications/utilities/mesh/generation/cv2DMesh/Make/files
new file mode 100755
index 0000000000..f7f70afa2a
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/Make/files
@@ -0,0 +1,12 @@
+#include CGAL_FILES
+
+CV2D.C
+insertFeaturePoints.C
+insertSurfaceNearestPointPairs.C
+insertSurfaceNearPointPairs.C
+insertBoundaryConformPointPairs.C
+CV2DIO.C
+shortEdgeFilter2D.C
+cv2DMesh.C
+
+EXE = $(FOAM_APPBIN)/cv2DMesh
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Make/options b/applications/utilities/mesh/generation/cv2DMesh/Make/options
new file mode 100755
index 0000000000..37c8ed4aea
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/Make/options
@@ -0,0 +1,43 @@
+EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+
+include $(GENERAL_RULES)/CGAL
+FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"'
+
+EXE_INC = \
+ ${EXE_FROUNDING_MATH} \
+ ${EXE_NDEBUG} \
+ ${CGAL_INC} \
+ -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
+ -I../cvMesh/vectorTools \
+ -IconformalVoronoi2DMesh/lnInclude \
+ -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
+ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/surfMesh/lnInclude \
+ -I$(LIB_SRC)/edgeMesh/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude \
+ -I$(LIB_SRC)/mesh/extrudeModel/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+
+EXE_LIBS = \
+ $(CGAL_LIBS) \
+ -lboost_thread \
+ -lmpfr \
+ -lextrude2DMesh \
+ -lextrudeModel \
+ -lcv2DMesh \
+ -lconformalVoronoiMesh \
+ -lmeshTools \
+ -lsurfMesh \
+ -ledgeMesh \
+ -ltriSurface \
+ -ldynamicMesh \
+ -ldecompositionMethods \
+ -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
+ -lsampling \
+ -lfileFormats
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files
new file mode 100755
index 0000000000..e35fec7b71
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files
@@ -0,0 +1,3 @@
+cv2DControls/cv2DControls.C
+
+LIB = $(FOAM_LIBBIN)/libcv2DMesh
\ No newline at end of file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options
new file mode 100755
index 0000000000..79be6f3a7d
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options
@@ -0,0 +1,3 @@
+EXE_INC =
+
+LIB_LIBS =
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.C b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.C
new file mode 100644
index 0000000000..1ec4442e10
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.C
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 ("objOutput", false)),
+
+ meshedSurfaceOutput_
+ (
+ motionControl_.lookupOrDefault("meshedSurfaceOutput", false)
+ ),
+
+ randomiseInitialGrid_(conformationControl_.lookup("randomiseInitialGrid")),
+ randomPerturbation_
+ (
+ readScalar(conformationControl_.lookup("randomPerturbation"))
+ ),
+
+ maxBoundaryConformingIter_
+ (
+ readLabel(conformationControl_.lookup("maxBoundaryConformingIter"))
+ ),
+
+ span_
+ (
+ max(mag(bb.max().x()), mag(bb.min().x()))
+ + max(mag(bb.max().y()), mag(bb.min().y()))
+ ),
+ span2_(Foam::sqr(span_)),
+
+ minEdgeLen_
+ (
+ readScalar(conformationControl_.lookup("minEdgeLenCoeff"))
+ *minCellSize_
+ ),
+ minEdgeLen2_(Foam::sqr(minEdgeLen_)),
+
+ maxNotchLen_
+ (
+ readScalar(conformationControl_.lookup("maxNotchLenCoeff"))
+ *minCellSize_
+ ),
+ maxNotchLen2_(Foam::sqr(maxNotchLen_)),
+
+ minNearPointDist_
+ (
+ readScalar(conformationControl_.lookup("minNearPointDistCoeff"))
+ *minCellSize_
+ ),
+ minNearPointDist2_(Foam::sqr(minNearPointDist_)),
+
+ ppDist_
+ (
+ readScalar(conformationControl_.lookup("pointPairDistanceCoeff"))
+ *minCellSize_
+ )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::cv2DControls::~cv2DControls()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+void Foam::cv2DControls::write(Ostream& os) const
+{
+ os.indentLevel() = 1;
+ os.precision(2);
+ os.flags(ios_base::scientific);
+
+ os << nl << "Outputting CV2D Mesher controls:" << nl
+ << token::BEGIN_BLOCK << nl
+ << indent << "minCellSize2_ : " << minCellSize2_ << nl
+ << indent << "span_ / span2_ : " << span_ << " / " << span2_ << nl
+ << indent << "maxNotchLen2_ : " << maxNotchLen2_ << nl
+ << indent << "minNearPointDist2_ : " << minNearPointDist2_ << nl
+ << indent << "nearWallAlignedDist2_ : " << nearWallAlignedDist2_ << nl
+ << indent << "ppDist_ : " << ppDist_ << nl
+ << indent << "minEdgeLen2_ : " << minEdgeLen2_ << nl
+ << token::END_BLOCK << endl;
+}
+
+
+// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const cv2DControls& s)
+{
+ s.write(os);
+ return os;
+}
+
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
new file mode 100644
index 0000000000..fb357dd62d
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
@@ -0,0 +1,260 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 patches(poly2DMesh.patchNames().size());
+ label countPatches = 0;
+
+ forAll(patches, patchI)
+ {
+ if (poly2DMesh.patchSizes()[patchI] != 0)
+ {
+ patches[countPatches] = new polyPatch
+ (
+ poly2DMesh.patchNames()[patchI],
+ poly2DMesh.patchSizes()[patchI],
+ poly2DMesh.patchStarts()[patchI],
+ countPatches,
+ pMesh.boundaryMesh(),
+ word::null
+ );
+
+ countPatches++;
+ }
+ }
+ patches.setSize(countPatches);
+ pMesh.addPatches(patches);
+
+ if (extrude)
+ {
+ Info<< "Begin extruding the polyMesh:" << endl;
+
+ {
+ // Point generator
+ autoPtr model(extrudeModel::New(extrusionDict));
+
+ extrude2DMesh extruder(pMesh, extrusionDict, model());
+
+ extruder.addFrontBackPatches();
+
+ polyTopoChange meshMod(pMesh.boundaryMesh().size());
+
+ extruder.setRefinement(meshMod);
+
+ autoPtr morphMap = meshMod.changeMesh(pMesh, false);
+
+ pMesh.updateMesh(morphMap);
+ }
+ }
+
+ if (!overwrite)
+ {
+ runTime++;
+ }
+ else
+ {
+ pMesh.setInstance("constant");
+ }
+
+ pMesh.write();
+
+ Info<< "Finished extruding in = "
+ << runTime.cpuTimeIncrement() << " s." << endl;
+
+ Info<< nl << "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict b/applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict
new file mode 100644
index 0000000000..0ef6d2a0c4
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict
@@ -0,0 +1,208 @@
+/*--------------------------------*- 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 cv2DMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+geometry
+{
+ laurence_clean_preciser.stl
+ {
+ name laurence_clean_preciser;
+ type closedTriSurfaceMesh;
+ //type triSurfaceMesh;
+ }
+// refinementBox
+// {
+// type searchableBox;
+// min (-0.5 0.35 -1000);
+// max (-0.5 0.35 1000);
+// }
+// refinementSphere
+// {
+// type searchableSphere;
+// centre (0.85 0.4 0.0);
+// radius 0.01;
+// }
+}
+
+surfaceConformation
+{
+ locationInMesh (-2.8 0.7 0.5);
+
+ pointPairDistanceCoeff 0.005;
+
+ minEdgeLenCoeff 0.005;
+
+ maxNotchLenCoeff 0.003;
+
+ minNearPointDistCoeff 0.0025;
+
+ maxQuadAngle 125;
+
+ // Insert near-boundary point mirror or point-pairs
+ insertSurfaceNearestPointPairs yes;
+
+ // Mirror near-boundary points rather than insert point-pairs
+ mirrorPoints no;
+
+ // Insert point-pairs vor dual-cell vertices very near the surface
+ insertSurfaceNearPointPairs yes;
+
+ // Maximum number of iterations used in boundaryConform.
+ maxBoundaryConformingIter 5;
+
+ geometryToConformTo
+ {
+ laurence_clean_preciser
+ {
+ featureMethod extendedFeatureEdgeMesh;
+ extendedFeatureEdgeMesh "laurence_clean_preciser.extendedFeatureEdgeMesh";
+ }
+ }
+
+ additionalFeatures
+ {
+ }
+
+ // Choose if to randomise the initial grid created by insertGrid.
+ randomiseInitialGrid yes;
+
+ // Perturbation fraction, 1 = cell-size.
+ randomPerturbation 0.1;
+
+}
+
+
+motionControl
+{
+ defaultCellSize 0.05;
+
+ // Assign a priority to all requests for cell sizes, the highest overrules.
+ defaultPriority 0;
+
+ cellSizeControlGeometry
+ {
+ laurence_clean_preciser
+ {
+ priority 1;
+ mode bothSides;
+ cellSizeFunction linearDistance;
+ linearDistanceCoeffs
+ {
+ distanceCellSize 0.05;
+ surfaceCellSize 0.01;
+ distance 0.5;
+ }
+ uniformCoeffs
+ {
+ cellSize 0.01;
+ }
+ }
+// refinementBox
+// {
+// priority 1;
+// mode outside;
+// cellSizeFunction linearDistance;
+// linearDistanceCoeffs
+// {
+// distanceCellSize 0.04;
+// surfaceCellSize 0.005;
+// distance 0.2;
+// }
+// }
+// refinementSphere
+// {
+// priority 1;
+// mode outside;
+// cellSizeFunction linearDistance;
+// linearDistanceCoeffs
+// {
+// distanceCellSize 0.04;
+// surfaceCellSize 0.005;
+// distance 0.2;
+// }
+// }
+ }
+
+ relaxationModel adaptiveLinear;
+
+ adaptiveLinearCoeffs
+ {
+ relaxationStart 0.5;
+ relaxationEnd 0.0;
+ }
+
+ objOutput no;
+
+ // Near-wall region where cells are aligned with the wall specified as a number
+ // of cell layers
+ nearWallAlignedDist 3;
+
+}
+
+shortEdgeFilter
+{
+ // Factor to multiply the average of a face's edge lengths by.
+ // If an edge of that face is smaller than that value then delete it.
+ shortEdgeFilterFactor 0.2;
+
+ // Weighting for the lengths of edges that are attached to the boundaries.
+ // Used when calculating the length of an edge. Default 2.0.
+ edgeAttachedToBoundaryFactor 2.0;
+}
+
+extrusion
+{
+ extrude on;
+
+ extrudeModel linearDirection;
+ //extrudeModel wedge;
+
+ patchInfo
+ {
+ //type empty;
+ //startFace
+ }
+
+ patchType empty;
+ //patchType wedge;
+
+ nLayers 1;
+
+ expansionRatio 1.0; //0.9;
+
+ linearDirectionCoeffs
+ {
+ direction (0 0 1);
+ thickness 0.1;
+ }
+
+ wedgeCoeffs
+ {
+ axisPt (0 0 0);
+ axis (1 0 0);
+ angle 10;
+ }
+
+ thickness 0.1;
+}
diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedFace.H b/applications/utilities/mesh/generation/cv2DMesh/indexedFace.H
new file mode 100644
index 0000000000..0c30fcc952
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/indexedFace.H
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 used to keep
+ track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedFace_H
+#define indexedFace_H
+
+#include
+#include "indexedVertex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+ Class indexedFace Declaration
+\*---------------------------------------------------------------------------*/
+
+template >
+class indexedFace
+:
+ public Fb
+{
+ // Private data
+
+ //- The index for this triangle face
+ // -1: triangle and changed and associated data needs updating
+ // >=0: index of triangles, set by external update algorithm
+ int index_;
+
+
+public:
+
+ enum faceTypes
+ {
+ UNCHANGED = 0,
+ CHANGED = -1,
+ SAVE_CHANGED = -2
+ };
+
+ typedef typename Fb::Vertex_handle Vertex_handle;
+ typedef typename Fb::Face_handle Face_handle;
+
+ template < typename TDS2 >
+ struct Rebind_TDS
+ {
+ typedef typename Fb::template Rebind_TDS::Other Fb2;
+ typedef indexedFace Other;
+ };
+
+
+ // Constructors
+
+ inline indexedFace();
+
+ inline indexedFace
+ (
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2
+ );
+
+ inline indexedFace
+ (
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2,
+ Face_handle n0,
+ Face_handle n1,
+ Face_handle n2
+ );
+
+
+ // Member Functions
+
+ inline void set_vertex(int i, Vertex_handle v);
+
+ inline void set_vertices();
+
+ inline void set_vertices
+ (
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2
+ );
+
+ inline int& faceIndex();
+
+ inline int faceIndex() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "indexedFaceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedFaceI.H b/applications/utilities/mesh/generation/cv2DMesh/indexedFaceI.H
new file mode 100644
index 0000000000..367cf94433
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/indexedFaceI.H
@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 CGAL::indexedFace::indexedFace()
+:
+ Fb(),
+ index_(CHANGED)
+{}
+
+
+template
+inline CGAL::indexedFace::indexedFace
+(
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2
+)
+:
+ Fb(v0, v1, v2),
+ index_(CHANGED)
+{}
+
+
+template
+inline CGAL::indexedFace::indexedFace
+(
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2,
+ Face_handle n0,
+ Face_handle n1,
+ Face_handle n2
+)
+:
+ Fb(v0, v1, v2, n0, n1, n2),
+ index_(CHANGED)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template
+inline void CGAL::indexedFace::set_vertex(int i, Vertex_handle v)
+{
+ index_ = CHANGED;
+ Fb::set_vertex(i, v);
+}
+
+
+template
+inline void CGAL::indexedFace::set_vertices()
+{
+ index_ = CHANGED;
+ Fb::set_vertices();
+}
+
+
+template
+inline void CGAL::indexedFace::set_vertices
+(
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2
+)
+{
+ index_ = CHANGED;
+ Fb::set_vertices(v0, v1, v2);
+}
+
+
+template
+inline int& CGAL::indexedFace::faceIndex()
+{
+ return index_;
+}
+
+
+template
+inline int CGAL::indexedFace::faceIndex() const
+{
+ return index_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H b/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H
new file mode 100644
index 0000000000..ae58a0bfe5
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 used to keep
+ track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedVertex_H
+#define indexedVertex_H
+
+#include
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+// Forward declaration of friend functions and operators
+
+template
+class indexedVertex;
+
+template
+bool pointPair
+(
+ const indexedVertex& v0,
+ const indexedVertex& v1
+);
+
+template
+bool boundaryTriangle
+(
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+);
+
+template
+bool outsideTriangle
+(
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+);
+
+/*---------------------------------------------------------------------------*\
+ Class indexedVertex Declaration
+\*---------------------------------------------------------------------------*/
+
+template >
+class indexedVertex
+:
+ public Vb
+{
+ // Private data
+
+ //- The index for this triangle vertex
+ int index_;
+
+ //- Index of pair-point :
+ // NEAR_BOUNDARY_POINT : internal near boundary point.
+ // INTERNAL_POINT : internal point.
+ // FAR_POINT : far-point.
+ // >= 0 : part of point-pair. Index of other point.
+ // Lowest numbered is inside one (master).
+ int type_;
+
+
+public:
+
+ enum pointTypes
+ {
+ NEAR_BOUNDARY_POINT = -4,
+ INTERNAL_POINT = -3,
+ MIRROR_POINT = -2,
+ FAR_POINT = -1
+ };
+
+ typedef typename Vb::Vertex_handle Vertex_handle;
+ typedef typename Vb::Face_handle Face_handle;
+ typedef typename Vb::Point Point;
+
+ template
+ struct Rebind_TDS
+ {
+ typedef typename Vb::template Rebind_TDS::Other Vb2;
+ typedef indexedVertex Other;
+ };
+
+
+ // Constructors
+
+ inline indexedVertex();
+
+ inline indexedVertex(const Point& p);
+
+ inline indexedVertex(const Point& p, const int index, const int& type);
+
+ inline indexedVertex(const Point& p, Face_handle f);
+
+ inline indexedVertex(Face_handle f);
+
+
+ // Member Functions
+
+ inline int& index();
+
+ inline int index() const;
+
+ inline int& type();
+
+ inline int type() const;
+
+ //- Is point a far-point
+ inline bool farPoint() const;
+
+ //- Is point internal, i.e. not on boundary
+ inline bool internalPoint() const;
+
+ //- Is point internal and near the boundary
+ inline bool nearBoundary() const;
+
+ //- Set the point to be near the boundary
+ inline void setNearBoundary();
+
+ //- Is point a mirror point
+ inline bool mirrorPoint() const;
+
+ //- Either master or slave of pointPair.
+ inline bool pairPoint() const;
+
+ //- Master of a pointPair is the lowest numbered one.
+ inline bool ppMaster() const;
+
+ //- Slave of a pointPair is the highest numbered one.
+ inline bool ppSlave() const;
+
+ //- Either original internal point or master of pointPair.
+ inline bool internalOrBoundaryPoint() const;
+
+ //- Is point near the boundary or part of the boundary definition
+ inline bool nearOrOnBoundary() const;
+
+
+ // Friend Functions
+
+ //- Do the two given vertices consitute a boundary point-pair
+ friend bool pointPair
+ (
+ const indexedVertex& v0,
+ const indexedVertex& v1
+ );
+
+ //- Do the three given vertices consitute a boundary triangle
+ friend bool boundaryTriangle
+ (
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+ );
+
+ //- Do the three given vertices consitute an outside triangle
+ friend bool outsideTriangle
+ (
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+ );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "indexedVertexI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedVertexI.H b/applications/utilities/mesh/generation/cv2DMesh/indexedVertexI.H
new file mode 100644
index 0000000000..2d6280274f
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/indexedVertexI.H
@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 CGAL::indexedVertex::indexedVertex()
+:
+ Vb(),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT)
+{}
+
+
+template
+inline CGAL::indexedVertex::indexedVertex(const Point& p)
+:
+ Vb(p),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT)
+{}
+
+
+template
+inline CGAL::indexedVertex::indexedVertex
+(
+ const Point& p,
+ const int index,
+ const int& type
+)
+:
+ Vb(p),
+ index_(index),
+ type_(type)
+{}
+
+
+template
+inline CGAL::indexedVertex::indexedVertex(const Point& p, Face_handle f)
+:
+ Vb(f, p),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT)
+{}
+
+
+template
+inline CGAL::indexedVertex::indexedVertex(Face_handle f)
+:
+ Vb(f),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template
+inline int& CGAL::indexedVertex::index()
+{
+ return index_;
+}
+
+
+template
+inline int CGAL::indexedVertex::index() const
+{
+ return index_;
+}
+
+
+template
+inline int& CGAL::indexedVertex::type()
+{
+ return type_;
+}
+
+
+template
+inline int CGAL::indexedVertex::type() const
+{
+ return type_;
+}
+
+
+template
+inline bool CGAL::indexedVertex::farPoint() const
+{
+ return type_ == FAR_POINT;
+}
+
+
+template
+inline bool CGAL::indexedVertex::internalPoint() const
+{
+ return type_ <= INTERNAL_POINT;
+}
+
+
+template
+inline bool CGAL::indexedVertex::nearBoundary() const
+{
+ return type_ == NEAR_BOUNDARY_POINT;
+}
+
+
+template
+inline void CGAL::indexedVertex::setNearBoundary()
+{
+ type_ = NEAR_BOUNDARY_POINT;
+}
+
+
+template
+inline bool CGAL::indexedVertex::mirrorPoint() const
+{
+ return type_ == MIRROR_POINT;
+}
+
+
+template
+inline bool CGAL::indexedVertex::pairPoint() const
+{
+ return type_ >= 0;
+}
+
+
+template
+inline bool CGAL::indexedVertex::ppMaster() const
+{
+ if (type_ > index_)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+
+template
+inline bool CGAL::indexedVertex::ppSlave() const
+{
+ if (type_ >= 0 && type_ < index_)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+
+template
+inline bool CGAL::indexedVertex::internalOrBoundaryPoint() const
+{
+ return internalPoint() || ppMaster();
+}
+
+
+template
+inline bool CGAL::indexedVertex::nearOrOnBoundary() const
+{
+ return pairPoint() || mirrorPoint() || nearBoundary();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
+
+template
+bool CGAL::pointPair
+(
+ const indexedVertex& v0,
+ const indexedVertex& v1
+)
+{
+ return v0.index_ == v1.type_ || v1.index_ == v0.type_;
+}
+
+
+template
+bool CGAL::boundaryTriangle
+(
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+)
+{
+ return (v0.pairPoint() && pointPair(v1, v2))
+ || (v1.pairPoint() && pointPair(v2, v0))
+ || (v2.pairPoint() && pointPair(v0, v1));
+}
+
+
+template
+bool CGAL::outsideTriangle
+(
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+)
+{
+ return (v0.farPoint() || v0.ppSlave())
+ || (v1.farPoint() || v1.ppSlave())
+ || (v2.farPoint() || v2.ppSlave());
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C
new file mode 100644
index 0000000000..7db9f2dae4
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C
@@ -0,0 +1,323 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 (1, pHit),
+ norm
+ );
+
+ insertPointPair
+ (
+ meshControls().ppDist(),
+ p,
+ toPoint2D(norm[0])
+ );
+ }
+
+ vit = Triangulation::Finite_vertices_iterator
+ (
+ CGAL::Filter_iterator
+ <
+ Triangulation::All_vertices_iterator,
+ Triangulation::Infinite_tester
+ >(finite_vertices_end(), vit.predicate(), vit.base())
+ );
+}
+
+
+bool Foam::CV2D::insertPointPairAtIntersection
+(
+ Triangulation::Finite_vertices_iterator& vit,
+ const point2D& defVert,
+ const point2D vertices[],
+ const scalar maxProtSize2
+)
+{
+ bool found = false;
+ point2D interPoint;
+ label interTri = -1;
+ label interHitSurface = -1;
+ scalar interDist2 = 0;
+
+ Face_circulator fcStart = incident_faces(vit);
+ Face_circulator fc = fcStart;
+ label vi = 0;
+
+ do
+ {
+ if (!is_infinite(fc))
+ {
+ pointIndexHit pHit;
+ label hitSurface = -1;
+
+ qSurf_.findSurfaceNearestIntersection
+ (
+ toPoint3D(defVert),
+ toPoint3D(vertices[vi]),
+ pHit,
+ hitSurface
+ );
+
+ if (pHit.hit())
+ {
+ scalar dist2 =
+ magSqr(toPoint2D(pHit.hitPoint()) - vertices[vi]);
+
+ // Check the point is further away than the furthest so far
+ if (dist2 > interDist2)
+ {
+ scalar mps2 = maxProtSize2;
+
+ // If this is a boundary triangle reset the tolerance
+ // to avoid finding a hit point very close to a boundary
+ // vertex
+ if (boundaryTriangle(fc))
+ {
+ mps2 = meshControls().maxNotchLen2();
+ }
+
+ if (dist2 > mps2)
+ {
+ found = true;
+ interPoint = toPoint2D(pHit.hitPoint());
+ interTri = pHit.index();
+ interDist2 = dist2;
+ interHitSurface = hitSurface;
+ }
+ }
+ }
+
+ vi++;
+ }
+ } while (++fc != fcStart);
+
+ if (found)
+ {
+ insertPointPair(vit, interPoint, interTri, interHitSurface);
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+
+Foam::label Foam::CV2D::insertBoundaryConformPointPairs
+(
+ const fileName& fName
+)
+{
+ label nIntersections = 0;
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ vit++
+ )
+ {
+ // Consider only those points part of point-pairs or near boundary
+ if (!vit->nearOrOnBoundary())
+ {
+ continue;
+ }
+
+ // Counter-clockwise circulator
+ Face_circulator fcStart = incident_faces(vit);
+ Face_circulator fc = fcStart;
+
+ bool infinite = false;
+ bool changed = false;
+
+ do
+ {
+ if (is_infinite(fc))
+ {
+ infinite = true;
+ break;
+ }
+ else if (fc->faceIndex() < Fb::UNCHANGED)
+ {
+ changed = true;
+ break;
+ }
+ } while (++fc != fcStart);
+
+ // If the dual-cell is connected to the infinite point or none of the
+ // faces whose circumcentres it uses have changed ignore
+ if (infinite || !changed) continue;
+
+ fc = fcStart;
+ label nVerts = 0;
+
+ do
+ {
+ vertices[nVerts++] = toPoint2D(circumcenter(fc));
+
+ if (nVerts == maxNvert)
+ {
+ break;
+ }
+ } while (++fc != fcStart);
+
+ // Check if dual-cell has a large number of faces in which case
+ // assumed to be in the far-field and reject
+ if (nVerts == maxNvert) continue;
+
+ // Set n+1 vertex to the first vertex for easy circulating
+ vertices[nVerts] = vertices[0];
+
+ // Convert triangle vertex to OpenFOAM point
+ point2DFromPoint defVert = toPoint2D(vit->point());
+
+ scalar maxProtSize2 = meshControls().maxNotchLen2();
+
+ if (vit->internalOrBoundaryPoint())
+ {
+ // Calculate metrics of the dual-cell
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ // The perimeter of the dual-cell
+ scalar perimeter = 0;
+
+ // Twice the area of the dual-cell
+ scalar areaT2 = 0;
+
+ for (int vi=0; vi 2*meshControls().minCellSize()
+ && areaT2 > meshControls().minCellSize2()
+ && cellWidth > 0.5*meshControls().minEdgeLen()
+ )
+ {
+ maxProtSize2 = 0.25*meshControls().maxNotchLen2();
+ }
+ }
+
+ if (insertPointPairAtIntersection(vit, defVert, vertices, maxProtSize2))
+ {
+ nIntersections++;
+ }
+ }
+
+ return nIntersections;
+}
+
+
+void Foam::CV2D::markNearBoundaryPoints()
+{
+ label count = 0;
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ vit++
+ )
+ {
+ if (vit->internalPoint())
+ {
+ point vert(toPoint3D(vit->point()));
+
+ pointIndexHit pHit;
+ label hitSurface = -1;
+
+ qSurf_.findSurfaceNearest
+ (
+ vert,
+ 4*meshControls().minCellSize2(),
+ pHit,
+ hitSurface
+ );
+
+ if (pHit.hit())
+ {
+ vit->setNearBoundary();
+ ++count;
+ }
+ }
+ }
+
+ Info<< count << " points marked as being near a boundary" << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C b/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C
new file mode 100644
index 0000000000..464557078f
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C
@@ -0,0 +1,394 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 & feMeshes
+ (
+ qSurf_.features()
+ );
+
+ if (feMeshes.empty())
+ {
+ WarningIn("CV2D::insertFeaturePoints")
+ << "Extended Feature Edge Mesh is empty so no feature points will "
+ << "be found." << nl
+ << " Use: featureMethod extendedFeatureEdgeMesh;" << nl
+ << endl;
+ }
+
+ forAll(feMeshes, i)
+ {
+ const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+ const edgeList& edges = feMesh.edges();
+ const pointField& points = feMesh.points();
+
+ if (debug)
+ {
+ label nConvex = feMesh.concaveStart() - feMesh.convexStart();
+ label nConcave = feMesh.mixedStart() - feMesh.concaveStart();
+ label nMixed = feMesh.nonFeatureStart() - feMesh.mixedStart();
+ label nExternal = feMesh.internalStart() - feMesh.externalStart();
+ label nInternal = feMesh.flatStart() - feMesh.internalStart();
+ label nFlat = feMesh.openStart() - feMesh.flatStart();
+ label nOpen = feMesh.multipleStart() - feMesh.openStart();
+ label nMultiple = edges.size() - feMesh.multipleStart();
+
+ Info<< "Inserting Feature Points:" << nl
+ << " Convex points: " << nConvex << nl
+ << " Concave points: " << nConcave << nl
+ << " Mixed points: " << nMixed << nl
+ << " External edges: " << nExternal << nl
+ << " Internal edges: " << nInternal << nl
+ << " Flat edges: " << nFlat << nl
+ << " Open edges: " << nOpen << nl
+ << " Multiple edges: " << nMultiple << endl;
+ }
+
+ // Args: (base point, normal)
+ // @todo allow user to input this
+ plane zPlane(vector(0, 0, z_), vector(0, 0, 1));
+
+ if (debug)
+ {
+ Info<< " plane: " << zPlane << " " << z_ << endl;
+ }
+
+ forAll(edges, edgeI)
+ {
+ const edge& e = feMesh.edges()[edgeI];
+
+ const point& ep0 = points[e.start()];
+ const point& ep1 = points[e.end()];
+
+ const linePointRef line(ep0, ep1);
+
+ scalar intersect = zPlane.lineIntersect(line);
+
+ point2D featPoint = toPoint2D(intersect * (ep1 - ep0) + ep0);
+
+ if (on2DLine(featPoint, line))
+ {
+ vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI));
+
+ vector2D cornerNormal = sum(fpn);
+ cornerNormal /= mag(cornerNormal);
+
+ if (debug)
+ {
+ Info<< nl << " line: " << line << nl
+ << " vec: " << line.vec() << nl
+ << " featurePoint: " << featPoint << nl
+ << " line length: " << line.mag() << nl
+ << " intersect: " << intersect << endl;
+ }
+
+ if
+ (
+ feMesh.getEdgeStatus(edgeI)
+ == extendedFeatureEdgeMesh::EXTERNAL
+ )
+ {
+ // Convex Point
+ Foam::point2D internalPt =
+ featPoint - meshControls().ppDist()*cornerNormal;
+
+ if (debug)
+ {
+ Info<< "PREC: " << internalPt << nl
+ << " : " << featPoint << nl
+ << " : " << meshControls().ppDist() << nl
+ << " : " << cornerNormal << endl;
+ }
+
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(internalPt),
+ nVert,
+ nVert + 1
+ )
+ );
+ label masterPtIndex = nVert++;
+
+ forAll(fpn, nI)
+ {
+ const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
+
+ plane planeN = plane(toPoint3D(featPoint), n3D);
+
+ Foam::point2D externalPt =
+ internalPt
+ + (
+ 2.0
+ * planeN.distance(toPoint3D(internalPt))
+ * fpn[nI]
+ );
+
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(externalPt),
+ nVert++,
+ masterPtIndex
+ )
+ );
+
+ if (debug)
+ {
+ Info<< " side point: " << externalPt << endl;
+ }
+ }
+
+ if (debug)
+ {
+ Info<< "Convex Point: " << featPoint << nl
+ << " corner norm: " << cornerNormal << nl
+ << " reference: " << internalPt << endl;
+ }
+ }
+ else if
+ (
+ feMesh.getEdgeStatus(edgeI)
+ == extendedFeatureEdgeMesh::INTERNAL
+ )
+ {
+ // Concave Point
+ Foam::point2D externalPt =
+ featPoint + meshControls().ppDist()*cornerNormal;
+
+ Foam::point2D refPt =
+ featPoint - meshControls().ppDist()*cornerNormal;
+
+ label slavePointIndex = 0;
+
+ scalar totalAngle =
+ radToDeg
+ (
+ constant::mathematical::pi
+ + acos(mag(fpn[0] & fpn[1]))
+ );
+
+ // Number of quadrants the angle should be split into
+ int nQuads =
+ int(totalAngle/meshControls().maxQuadAngle()) + 1;
+
+ // The number of additional master points needed to
+ //obtain the required number of quadrants.
+ int nAddPoints = min(max(nQuads - 2, 0), 2);
+
+ // index of reflMaster
+ label reflectedMaster = nVert + 2 + nAddPoints;
+
+ if (debug)
+ {
+ Info<< "Concave Point: " << featPoint << nl
+ << " corner norm: " << cornerNormal << nl
+ << " external: " << externalPt << nl
+ << " reference: " << refPt << nl
+ << " angle: " << totalAngle << nl
+ << " nQuads: " << nQuads << nl
+ << " nAddPoints: " << nAddPoints << endl;
+ }
+
+ forAll(fpn, nI)
+ {
+ const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
+
+ plane planeN = plane(toPoint3D(featPoint), n3D);
+
+ Foam::point2D internalPt =
+ externalPt
+ - (
+ 2.0
+ * planeN.distance(toPoint3D(externalPt))
+ * fpn[nI]
+ );
+
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(internalPt),
+ nVert,
+ reflectedMaster
+ )
+ );
+ slavePointIndex = nVert++;
+
+ if (debug)
+ {
+ Info<< "Internal Point: " << internalPt << endl;
+ }
+ }
+
+ if (nAddPoints == 1)
+ {
+ // One additional point is the reflection of the slave
+ // point, i.e., the original reference point
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(refPt),
+ nVert++,
+ reflectedMaster
+ )
+ );
+
+ if (debug)
+ {
+ Info<< "ref Point: " << refPt << endl;
+ }
+ }
+ else if (nAddPoints == 2)
+ {
+ point2D reflectedAa =
+ refPt - ((featPoint - externalPt) & fpn[1])*fpn[1];
+
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(reflectedAa),
+ nVert++,
+ reflectedMaster
+ )
+ );
+
+ point2D reflectedBb =
+ refPt - ((featPoint - externalPt) & fpn[0])*fpn[0];
+
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(reflectedBb),
+ nVert++,
+ reflectedMaster
+ )
+ );
+
+ if (debug)
+ {
+ Info<< "refA Point: " << reflectedAa << nl
+ << "refb Point: " << reflectedBb << endl;
+ }
+ }
+
+ featurePoints_.push_back
+ (
+ Vb
+ (
+ toPoint(externalPt),
+ nVert++,
+ slavePointIndex
+ )
+ );
+ }
+ else
+ {
+ WarningIn("void Foam::CV2D::insertFeaturePoints()")
+ << "Feature Edge " << edges[edgeI] << nl
+ << " points(" << points[edges[edgeI].start()]
+ << ", " << points[edges[edgeI].end()] << ")" << nl
+ << " is not labelled as either concave or convex, it"
+ << " is labelled as (#2 = flat): "
+ << feMesh.getEdgeStatus(edgeI) << endl;
+ }
+ }
+ else
+ {
+ WarningIn("void Foam::CV2D::insertFeaturePoints()")
+ << "Point " << featPoint << " is not on the line "
+ << line << endl;
+ }
+ }
+ }
+
+ // Insert the feature points.
+ reinsertFeaturePoints();
+
+ if (meshControls().objOutput())
+ {
+ writePoints("feat_allPoints.obj", false);
+ writeFaces("feat_allFaces.obj", false);
+ writeFaces("feat_faces.obj", true);
+ writeTriangles("feat_triangles.obj", true);
+ }
+}
+
+
+void Foam::CV2D::reinsertFeaturePoints()
+{
+ for
+ (
+ std::list::iterator vit=featurePoints_.begin();
+ vit != featurePoints_.end();
+ ++vit
+ )
+ {
+ insertPoint
+ (
+ toPoint2D(vit->point()),
+ vit->index(),
+ vit->type()
+ );
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
new file mode 100644
index 0000000000..700041f814
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 (1, pHit),
+ norm
+ );
+
+ insertPointPair
+ (
+ meshControls().ppDist(),
+ toPoint2D(pHit.hitPoint()),
+ toPoint2D(norm[0])
+ );
+
+ nNearPoints++;
+
+ // Correct the edge iterator for the change in the
+ // number of edges following the point-pair insertion
+ eit = Finite_edges_iterator
+ (
+ finite_edges_end().base(),
+ eit.predicate(),
+ eit.base()
+ );
+ }
+ }
+ }
+ }
+ }
+
+ Info<< nNearPoints << " point-pairs inserted" << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C
new file mode 100644
index 0000000000..49249241e7
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C
@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 & nearSurfacePoints,
+ const DynamicList& surfacePoints,
+ const DynamicList& surfaceTris,
+ const DynamicList& surfaceHits,
+ const fileName fName
+)
+{
+ if (meshControls().mirrorPoints())
+ {
+ forAll(surfacePoints, ppi)
+ {
+ insertMirrorPoint
+ (
+ nearSurfacePoints[ppi],
+ surfacePoints[ppi]
+ );
+ }
+ }
+ else
+ {
+ forAll(surfacePoints, ppi)
+ {
+ pointIndexHit pHit
+ (
+ true,
+ toPoint3D(surfacePoints[ppi]),
+ surfaceTris[ppi]
+ );
+
+ vectorField norm(1);
+ qSurf_.geometry()[surfaceHits[ppi]].getNormal
+ (
+ List(1, pHit),
+ norm
+ );
+
+ insertPointPair
+ (
+ meshControls().ppDist(),
+ surfacePoints[ppi],
+ toPoint2D(norm[0])
+ );
+ }
+ }
+
+ Info<< surfacePoints.size() << " point-pairs inserted" << endl;
+
+ if (meshControls().objOutput())
+ {
+ OFstream str(fName);
+ label vertI = 0;
+
+ forAll(surfacePoints, ppi)
+ {
+ meshTools::writeOBJ(str, toPoint3D(surfacePoints[ppi]));
+ vertI++;
+ }
+
+ Info<< "insertPointPairs: Written " << surfacePoints.size()
+ << " inserted point-pair locations to file "
+ << str.name() << endl;
+ }
+}
+
+
+void Foam::CV2D::insertSurfaceNearestPointPairs()
+{
+ Info<< "insertSurfaceNearestPointPairs: ";
+
+ label nSurfacePointsEst =
+ min
+ (
+ number_of_vertices(),
+ size_t(10*sqrt(scalar(number_of_vertices())))
+ );
+
+ DynamicList nearSurfacePoints(nSurfacePointsEst);
+ DynamicList surfacePoints(nSurfacePointsEst);
+ DynamicList surfaceTris(nSurfacePointsEst);
+ DynamicList surfaceHits(nSurfacePointsEst);
+
+ // Local references to surface mesh addressing
+// const pointField& localPoints = qSurf_.localPoints();
+// const labelListList& edgeFaces = qSurf_.edgeFaces();
+// const vectorField& faceNormals = qSurf_.faceNormals();
+// const labelListList& faceEdges = qSurf_.faceEdges();
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ vit++
+ )
+ {
+ if (vit->internalPoint())
+ {
+ point2DFromPoint vert(toPoint2D(vit->point()));
+
+ pointIndexHit pHit;
+ label hitSurface = -1;
+
+ qSurf_.findSurfaceNearest
+ (
+ toPoint3D(vert),
+ 4*meshControls().minCellSize2(),
+ pHit,
+ hitSurface
+ );
+
+ if (pHit.hit())
+ {
+ vit->setNearBoundary();
+
+ // Reference to the nearest triangle
+// const labelledTri& f = qSurf_[hitSurface];
+
+// // Find where point is on triangle.
+// // Note tolerance needed is relative one
+// // (used in comparing normalized [0..1] triangle coordinates).
+// label nearType, nearLabel;
+// triPointRef
+// (
+// localPoints[f[0]],
+// localPoints[f[1]],
+// localPoints[f[2]]
+// ).classify(pHit.hitPoint(), nearType, nearLabel);
+
+// // If point is on a edge check if it is an internal feature
+
+// bool internalFeatureEdge = false;
+
+// if (nearType == triPointRef::EDGE)
+// {
+// label edgeI = faceEdges[hitSurface][nearLabel];
+// const labelList& eFaces = edgeFaces[edgeI];
+
+// if
+// (
+// eFaces.size() == 2
+// && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
+// < -0.2
+// )
+// {
+// internalFeatureEdge = true;
+// }
+// }
+
+ if (dualCellSurfaceIntersection(vit)) //&& !internalFeatureEdge)
+ {
+ nearSurfacePoints.append(vert);
+ surfacePoints.append(toPoint2D(pHit.hitPoint()));
+ surfaceTris.append(pHit.index());
+ surfaceHits.append(hitSurface);
+ }
+ }
+ }
+ }
+
+ insertPointPairs
+ (
+ nearSurfacePoints,
+ surfacePoints,
+ surfaceTris,
+ surfaceHits,
+ "surfaceNearestIntersections.obj"
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C
new file mode 100644
index 0000000000..b125351314
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C
@@ -0,0 +1,533 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 ("edgeAttachedToBoundaryFactor", 2.0)
+ ),
+ patchNames_(wordList()),
+ patchSizes_(labelList()),
+ mapEdgesRegion_(),
+ indirectPatchEdge_()
+{
+ point2DField points2D;
+ faceList faces;
+
+ cv2Dmesh.calcDual
+ (
+ points2D,
+ faces,
+ patchNames_,
+ patchSizes_,
+ mapEdgesRegion_,
+ indirectPatchEdge_
+ );
+
+ pointField points(points2D.size());
+ forAll(points, ip)
+ {
+ points[ip] = cv2Dmesh.toPoint3D(points2D[ip]);
+ }
+
+ if (debug)
+ {
+ OFstream str("indirectPatchEdges.obj");
+ label count = 0;
+
+ Info<< "Writing indirectPatchEdges to " << str.name() << endl;
+
+ forAllConstIter(EdgeMap, indirectPatchEdge_, iter)
+ {
+ const edge& e = iter.key();
+
+ meshTools::writeOBJ
+ (
+ str,
+ points[e.start()],
+ points[e.end()],
+ count
+ );
+ }
+ }
+
+ points2D.clear();
+
+ ms_ = MeshedSurface(xferMove(points), xferMove(faces));
+
+ Info<< "Meshed surface stats before edge filtering :" << endl;
+ ms_.writeStats(Info);
+
+ if (debug)
+ {
+ writeInfo(Info);
+
+ ms_.write("MeshedSurface_preFilter.obj");
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::shortEdgeFilter2D::~shortEdgeFilter2D()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+void
+Foam::shortEdgeFilter2D::filter()
+{
+ // These are global indices.
+ const pointField& points = ms_.points();
+ const edgeList& edges = ms_.edges();
+ const faceList& faces = ms_.faces();
+ const labelList& meshPoints = ms_.meshPoints();
+ const labelList& boundaryPoints = ms_.boundaryPoints();
+
+ label maxChain = 0;
+ label nPointsToRemove = 0;
+
+ labelList pointsToRemove(ms_.points().size(), -1);
+
+ // List of number of vertices in a face.
+ labelList newFaceVertexCount(faces.size(), -1);
+ forAll(faces, faceI)
+ {
+ newFaceVertexCount[faceI] = faces[faceI].size();
+ }
+
+ // Check if the point is a boundary point. Flag if it is so that
+ // it will not be deleted.
+ boolList boundaryPointFlags(points.size(), false);
+ // This has been removed, otherwise small edges on the boundary are not
+ // removed.
+ /* forAll(boundaryPointFlags, pointI)
+ {
+ forAll(boundaryPoints, bPoint)
+ {
+ if (meshPoints[boundaryPoints[bPoint]] == pointI)
+ {
+ boundaryPointFlags[pointI] = true;
+ }
+ }
+ }*/
+
+ // Check if an edge has a boundary point. It it does the edge length
+ // will be doubled when working out its length.
+ Info<< " Marking edges attached to boundaries." << endl;
+ boolList edgeAttachedToBoundary(edges.size(), false);
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+ const label startVertex = e.start();
+ const label endVertex = e.end();
+
+ forAll(boundaryPoints, bPoint)
+ {
+ if
+ (
+ boundaryPoints[bPoint] == startVertex
+ || boundaryPoints[bPoint] == endVertex
+ )
+ {
+ edgeAttachedToBoundary[edgeI] = true;
+ }
+ }
+ }
+
+ forAll(edges, edgeI)
+ {
+ const edge& e = edges[edgeI];
+
+ // get the vertices of that edge.
+ const label startVertex = e.start();
+ const label endVertex = e.end();
+
+ scalar edgeLength =
+ mag
+ (
+ points[meshPoints[startVertex]]
+ - points[meshPoints[endVertex]]
+ );
+
+ if (edgeAttachedToBoundary[edgeI])
+ {
+ edgeLength *= edgeAttachedToBoundaryFactor_;
+ }
+
+ scalar shortEdgeFilterValue = 0.0;
+
+ const labelList& psEdges = ms_.pointEdges()[startVertex];
+ const labelList& peEdges = ms_.pointEdges()[endVertex];
+
+ forAll(psEdges, psEdgeI)
+ {
+ const edge& psE = edges[psEdges[psEdgeI]];
+ if (edgeI != psEdges[psEdgeI])
+ {
+ shortEdgeFilterValue +=
+ mag
+ (
+ points[meshPoints[psE.start()]]
+ - points[meshPoints[psE.end()]]
+ );
+ }
+ }
+
+ forAll(peEdges, peEdgeI)
+ {
+ const edge& peE = edges[peEdges[peEdgeI]];
+ if (edgeI != peEdges[peEdgeI])
+ {
+ shortEdgeFilterValue +=
+ mag
+ (
+ points[meshPoints[peE.start()]]
+ - points[meshPoints[peE.end()]]
+ );
+ }
+ }
+
+ shortEdgeFilterValue *=
+ shortEdgeFilterFactor_
+ /(psEdges.size() + peEdges.size() - 2);
+
+
+ edge lookupInPatchEdge
+ (
+ meshPoints[startVertex],
+ meshPoints[endVertex]
+ );
+
+ if
+ (
+ edgeLength < shortEdgeFilterValue
+ || indirectPatchEdge_.found(lookupInPatchEdge)
+ )
+ {
+ bool flagDegenerateFace = false;
+ const labelList& pFaces = ms_.pointFaces()[startVertex];
+
+ forAll(pFaces, pFaceI)
+ {
+ const face& f = ms_.localFaces()[pFaces[pFaceI]];
+ forAll(f, fp)
+ {
+ // If the edge is part of this face...
+ if (f[fp] == endVertex)
+ {
+ // If deleting vertex would create a triangle, don't!
+ if (newFaceVertexCount[pFaces[pFaceI]] < 4)
+ {
+ flagDegenerateFace = true;
+ }
+ else
+ {
+ newFaceVertexCount[pFaces[pFaceI]]--;
+ }
+ }
+ // If the edge is not part of this face...
+ else
+ {
+ // Deleting vertex of a triangle...
+ if (newFaceVertexCount[pFaces[pFaceI]] < 3)
+ {
+ flagDegenerateFace = true;
+ }
+ }
+ }
+ }
+
+ // This if statement determines whether a point should be deleted.
+ if
+ (
+ pointsToRemove[meshPoints[startVertex]] == -1
+ && pointsToRemove[meshPoints[endVertex]] == -1
+ && !boundaryPointFlags[meshPoints[startVertex]]
+ && !flagDegenerateFace
+ )
+ {
+ pointsToRemove[meshPoints[startVertex]] = meshPoints[endVertex];
+ ++nPointsToRemove;
+ }
+ }
+ }
+
+ label totalNewPoints = points.size() - nPointsToRemove;
+
+ pointField newPoints(totalNewPoints, vector::zero);
+ labelList newPointNumbers(points.size(), -1);
+ label numberRemoved = 0;
+
+ forAll(points, pointI)
+ {
+ // If the point is NOT going to be removed.
+ if (pointsToRemove[pointI] == -1)
+ {
+ newPoints[pointI - numberRemoved] = points[pointI];
+ newPointNumbers[pointI] = pointI - numberRemoved;
+ }
+ else
+ {
+ numberRemoved++;
+ }
+ }
+
+ // Need a new faceList
+ faceList newFaces(faces.size());
+ label newFaceI = 0;
+
+ labelList newFace;
+ label newFaceSize = 0;
+
+ // Now need to iterate over the faces and remove points. Global index.
+ forAll(faces, faceI)
+ {
+ const face& f = faces[faceI];
+
+ newFace.clear();
+ newFace.setSize(f.size());
+ newFaceSize = 0;
+
+ forAll(f, fp)
+ {
+ label pointI = f[fp];
+ // If not removing the point, then add it to the new face.
+ if (pointsToRemove[pointI] == -1)
+ {
+ newFace[newFaceSize++] = newPointNumbers[pointI];
+ }
+ else
+ {
+ label newPointI = pointsToRemove[pointI];
+ // Replace deleted point with point that it is being
+ // collapsed to.
+ if
+ (
+ f.nextLabel(fp) != newPointI
+ && f.prevLabel(fp) != newPointI
+ )
+ {
+ label pChain = newPointI;
+ label totalChain = 0;
+ for (label nChain = 0; nChain <= totalChain; ++nChain)
+ {
+ if (newPointNumbers[pChain] != -1)
+ {
+ newFace[newFaceSize++] = newPointNumbers[pChain];
+ newPointNumbers[pointI] = newPointNumbers[pChain];
+ maxChain = max(totalChain, maxChain);
+ }
+ else
+ {
+ WarningIn("shortEdgeFilter")
+ << "Point " << pChain
+ << " marked for deletion as well as point "
+ << pointI << nl
+ << " Incrementing maxChain by 1 from "
+ << totalChain << " to " << totalChain + 1
+ << endl;
+ totalChain++;
+ }
+ pChain = pointsToRemove[pChain];
+ }
+ }
+ else
+ {
+ if (newPointNumbers[newPointI] != -1)
+ {
+ newPointNumbers[pointI] = newPointNumbers[newPointI];
+ }
+ }
+ }
+ }
+
+ newFace.setSize(newFaceSize);
+
+ if (newFace.size() > 2)
+ {
+ newFaces[newFaceI++] = face(newFace);
+ }
+ else
+ {
+ FatalErrorIn("shortEdgeFilter")
+ << "Only " << newFace.size() << " in face " << faceI
+ << exit(FatalError);
+ }
+ }
+
+ newFaces.setSize(newFaceI);
+
+ MeshedSurface fMesh
+ (
+ xferMove(newPoints),
+ xferMove(newFaces),
+ xferCopy(List())
+ );
+
+ const Map& fMeshPointMap = fMesh.meshPointMap();
+
+ // Reset patchSizes_
+ patchSizes_.clear();
+ patchSizes_.setSize(patchNames_.size(), 0);
+
+ label equalEdges = 0;
+ label notFound = 0;
+ label matches = 0;
+ label negativeLabels = 0;
+
+ forAll(newPointNumbers, pointI)
+ {
+ if (newPointNumbers[pointI] == -1)
+ {
+ WarningIn("shortEdgeFilter")
+ << pointI << " will be deleted and " << newPointNumbers[pointI]
+ << ", so it will not be replaced. "
+ << "This will cause edges to be deleted." << endl;
+ }
+ }
+
+ // Create new EdgeMap.
+ Info<< "Creating new EdgeMap." << endl;
+ EdgeMap newMapEdgesRegion(mapEdgesRegion_.size());
+
+ for
+ (
+ label bEdgeI = ms_.nInternalEdges();
+ bEdgeI < edges.size();
+ ++bEdgeI
+ )
+ {
+ label p1 = meshPoints[edges[bEdgeI][0]];
+ label p2 = meshPoints[edges[bEdgeI][1]];
+
+ edge e(p1, p2);
+
+ if (mapEdgesRegion_.found(e))
+ {
+ if
+ (
+ newPointNumbers[p1] != -1
+ && newPointNumbers[p2] != -1
+ )
+ {
+ if (newPointNumbers[p1] != newPointNumbers[p2])
+ {
+ label region = mapEdgesRegion_.find(e)();
+ newMapEdgesRegion.insert
+ (
+ edge
+ (
+ fMeshPointMap[newPointNumbers[p1]],
+ fMeshPointMap[newPointNumbers[p2]]
+ ),
+ region
+ );
+ patchSizes_[region]++;
+ matches++;
+ }
+ else
+ {
+ equalEdges++;
+ }
+ }
+ else
+ {
+ negativeLabels++;
+ }
+ }
+ else
+ {
+ notFound++;
+ }
+ }
+
+ if (debug)
+ {
+ Info<< "EdgeMapping :" << nl
+ << " Matches : " << matches << nl
+ << " Equal : " << equalEdges << nl
+ << " Negative : " << negativeLabels << nl
+ << " Not Found: " << notFound << endl;
+ }
+
+ mapEdgesRegion_.transfer(newMapEdgesRegion);
+
+ ms_.transfer(fMesh);
+
+ Info<< " Maximum number of chained collapses = " << maxChain << endl;
+
+ if (debug)
+ {
+ writeInfo(Info);
+ }
+}
+
+
+void Foam::shortEdgeFilter2D::writeInfo(Ostream& os)
+{
+ os << "Short Edge Filtering Information:" << nl
+ << " shortEdgeFilterFactor: " << shortEdgeFilterFactor_ << nl
+ << " edgeAttachedToBoundaryFactor: " << edgeAttachedToBoundaryFactor_
+ << endl;
+
+ forAll(patchNames_, patchI)
+ {
+ os << " Patch " << patchNames_[patchI]
+ << ", size " << patchSizes_[patchI] << endl;
+ }
+
+ os << " There are " << mapEdgesRegion_.size()
+ << " boundary edges." << endl;
+
+ os << " Mesh Info:" << nl
+ << " Points: " << ms_.nPoints() << nl
+ << " Faces: " << ms_.size() << nl
+ << " Edges: " << ms_.nEdges() << nl
+ << " Internal: " << ms_.nInternalEdges() << nl
+ << " External: " << ms_.nEdges() - ms_.nInternalEdges()
+ << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H
new file mode 100644
index 0000000000..53b248fc3d
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H
@@ -0,0 +1,133 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 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 ms_;
+
+ const scalar shortEdgeFilterFactor_;
+
+ const scalar edgeAttachedToBoundaryFactor_;
+
+ wordList patchNames_;
+
+ labelList patchSizes_;
+
+ EdgeMap mapEdgesRegion_;
+
+ EdgeMap indirectPatchEdge_;
+
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ shortEdgeFilter2D(const shortEdgeFilter2D&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const shortEdgeFilter2D&);
+
+
+public:
+
+ //- Runtime type information
+ ClassName("shortEdgeFilter2D");
+
+ // Constructors
+
+ shortEdgeFilter2D(const CV2D& cv2Dmesh, const dictionary& dict);
+
+
+ //- Destructor
+ ~shortEdgeFilter2D();
+
+
+ // Access Functions
+
+ const wordList& patchNames() const
+ {
+ return patchNames_;
+ }
+
+ const labelList& patchSizes() const
+ {
+ return patchSizes_;
+ }
+
+ const EdgeMap& mapEdgesRegion() const
+ {
+ return mapEdgesRegion_;
+ }
+
+ const MeshedSurface& fMesh() const
+ {
+ return ms_;
+ }
+
+
+ // Member Functions
+
+ void filter();
+
+ void writeInfo(Ostream& os);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/Allwclean b/applications/utilities/mesh/generation/cvMesh/Allwclean
new file mode 100755
index 0000000000..3eff84006a
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/Allwclean
@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+set -x
+
+wclean libso conformalVoronoiMesh
+wclean
+wclean cvMeshSurfaceSimplify
+wclean cvMeshBackgroundMesh
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cvMesh/Allwmake b/applications/utilities/mesh/generation/cvMesh/Allwmake
new file mode 100755
index 0000000000..d88e8cee4a
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/Allwmake
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+set -x
+
+wmake libso conformalVoronoiMesh
+wmake
+#wmake cvMeshBackgroundMesh
+(cd cvMeshSurfaceSimplify && ./Allwmake)
+wmake cellSizeAndAlignmentGrid
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cvMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/Make/files
new file mode 100644
index 0000000000..425e0f5b83
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/Make/files
@@ -0,0 +1,3 @@
+cvMesh.C
+
+EXE = $(FOAM_APPBIN)/cvMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
new file mode 100644
index 0000000000..9461fa3725
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/Make/options
@@ -0,0 +1,38 @@
+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} \
+ -IconformalVoronoiMesh/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+ -I$(LIB_SRC)/edgeMesh/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -IvectorTools
+
+EXE_LIBS = \
+ $(CGAL_LIBS) \
+ -lboost_thread \
+ -lmpfr \
+ -lconformalVoronoiMesh \
+ -lmeshTools \
+ -ldecompositionMethods \
+ -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
+ -ledgeMesh \
+ -lfileFormats \
+ -ltriSurface \
+ -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..454416299d
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C
@@ -0,0 +1,716 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012-2013 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()] = vit->alignment();
+ }
+
+ 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;
+}
+
+
+void refine
+(
+ cellShapeControlMesh& mesh,
+ const conformationSurfaces& geometryToConformTo,
+ const label maxRefinementIterations,
+ const scalar defaultCellSize
+)
+{
+ 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,
+ Vb::vtInternal
+ );
+ }
+ }
+}
+
+
+int main(int argc, char *argv[])
+{
+ #include "setRootCase.H"
+ #include "createTime.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ label maxRefinementIterations = 2;
+ label maxSmoothingIterations = 200;
+ scalar minResidual = 0;
+ scalar defaultCellSize = 0.001;
+ 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(),
+ Vb::vtInternalNearBoundary
+ );
+ }
+ }
+ else
+ {
+ CellSizeDelaunay::Vertex_handle vh = mesh.insert
+ (
+ points[pI],
+ defaultCellSize,
+ pointAlignment(),
+ Vb::vtInternalNearBoundary
+ );
+ }
+ }
+ }
+
+
+ // Refine the mesh
+ refine
+ (
+ mesh,
+ geometryToConformTo,
+ maxRefinementIterations,
+ defaultCellSize
+ );
+
+
+ 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->nearBoundary())
+ {
+ fixedAlignments[vit->index()] = vit->alignment();
+ }
+ }
+
+ 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];
+
+ // Enforce the boundary conditions
+ const triad& fixedAlignment = fixedAlignments[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();
+
+ label nFixed = 0;
+
+ forAll(fixedAlignment, dirI)
+ {
+ if (fixedAlignment.set(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];
+ }
+ }
+ }
+
+ for (direction dir = 0; dir < 3; ++dir)
+ {
+ if
+ (
+ newTriad.set(dir)
+ && oldTriad.set(dir)
+ //&& !fixedAlignment.set(dir)
+ )
+ {
+ scalar dotProd = (oldTriad[dir] & newTriad[dir]);
+ 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/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