Reintroducing CV meshers after removal by a commit in master, originating in the dsmc branch.

applications/utilities/mesh/generation/CV2DMesher/CGALTriangulation2Ddefs.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Typedefs
+    CGALTriangulation2Ddefs
+
+Description
+    CGAL data structures used for 2D Delaunay meshing.
+
+    Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel
+    otherwise the more robust but much less efficient
+    Exact_predicates_exact_constructions will be used.
+
+    Define CGAL_HIERARCHY to use hierarchical Delaunay triangulation which is
+    faster but uses more memory than the standard Delaunay triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CGALTriangulation2Ddefs_H
+#define CGALTriangulation2Ddefs_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CGAL/Delaunay_triangulation_2.h"
+
+#include "indexedVertex.H"
+#include "indexedFace.H"
+
+#ifdef CGAL_INEXACT
+    // Fast kernel using a double as the storage type but the triangulation
+    // may fail
+    #include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
+    typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+#else
+    // Very robust but expensive kernel
+    #include "CGAL/Exact_predicates_exact_constructions_kernel.h"
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+#endif
+
+typedef CGAL::indexedVertex<K>                           Vb;
+typedef CGAL::indexedFace<K>                             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<Vb>  Vbh;
+    typedef CGAL::Triangulation_data_structure_2<Vbh, Fb>    Tds;
+    typedef CGAL::Delaunay_triangulation_2<K, Tds>           Triangulation;
+    typedef CGAL::Triangulation_hierarchy_2<Triangulation>   HTriangulation;
+#else
+    // Data structures for standard Delaunay triangulation
+    typedef CGAL::Triangulation_data_structure_2<Vb, Fb>     Tds;
+    typedef CGAL::Delaunay_triangulation_2<K, Tds>           Triangulation;
+    typedef Triangulation                                    HTriangulation;
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/CV2D.C

@@ -0,0 +1,567 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+#include "Random.H"
+#include "transform.H"
+#include "IFstream.H"
+#include "uint.H"
+#include "ulong.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV2D::insertBoundingBox()
+{
+    Info<< "insertBoundingBox: creating bounding mesh" << endl;
+    scalar bigSpan = 10*tols_.span;
+    insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT);
+    insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT);
+    insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT);
+    insertPoint(point2D(bigSpan, bigSpan), Vb::FAR_POINT);
+}
+
+
+void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh)
+{
+    int i;
+    Face_handle f = vh->face(), next, start(f);
+
+    do
+    {
+        i=f->index(vh);
+        if (!is_infinite(f))
+        {
+            if (!internal_flip(f, cw(i))) external_flip(f, i);
+            if (f->neighbor(i) == start) start = f;
+        }
+        f = f->neighbor(cw(i));
+    } while (f != start);
+}
+
+void Foam::CV2D::external_flip(Face_handle& f, int i)
+{
+    Face_handle n = f->neighbor(i);
+
+    if
+    (
+        CGAL::ON_POSITIVE_SIDE
+     != side_of_oriented_circle(n, f->vertex(i)->point())
+    ) return;
+
+    flip(f, i);
+    i = n->index(f->vertex(i));
+    external_flip(n, i);
+}
+
+bool Foam::CV2D::internal_flip(Face_handle& f, int i)
+{
+    Face_handle n = f->neighbor(i);
+
+    if
+    (
+        CGAL::ON_POSITIVE_SIDE
+     != side_of_oriented_circle(n, f->vertex(i)->point())
+    )
+    {
+        return false;
+    }
+
+    flip(f, i);
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CV2D::CV2D
+(
+    const dictionary& controlDict,
+    const querySurface& qSurf
+)
+:
+    HTriangulation(),
+    qSurf_(qSurf),
+    controls_(controlDict),
+    tols_(controlDict, controls_.minCellSize, qSurf.bb()),
+    z_((1.0/3.0)*(qSurf_.bb().min().z() + qSurf_.bb().max().z())),
+    startOfInternalPoints_(0),
+    startOfSurfacePointPairs_(0),
+    startOfBoundaryConformPointPairs_(0)
+{
+    insertBoundingBox();
+    insertFeaturePoints();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::CV2D::~CV2D()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::CV2D::insertPoints
+(
+    const point2DField& points,
+    const scalar nearness
+)
+{
+    Info<< "insertInitialPoints(const point2DField& points): ";
+
+    startOfInternalPoints_ = number_of_vertices();
+    label nVert = startOfInternalPoints_;
+
+    // Add the points and index them
+    forAll(points, i)
+    {
+        const point2D& p = points[i];
+
+        if (qSurf_.wellInside(toPoint3D(p), nearness))
+        {
+            insert(toPoint(p))->index() = nVert++;
+        }
+        else
+        {
+            Warning
+                << "Rejecting point " << p << " outside surface" << endl;
+        }
+    }
+
+    Info<< nVert << " vertices inserted" << endl;
+
+    if (controls_.writeInitialTriangulation)
+    {
+        // Checking validity of triangulation
+        assert(is_valid());
+
+        writeTriangles("initial_triangles.obj", true);
+        writeFaces("initial_faces.obj", true);
+    }
+}
+
+
+void Foam::CV2D::insertPoints(const fileName& pointFileName)
+{
+    IFstream pointsFile(pointFileName);
+
+    if (pointsFile.good())
+    {
+        insertPoints(point2DField(pointsFile), 0.5*controls_.minCellSize2);
+    }
+    else
+    {
+        FatalErrorIn("insertInitialPoints")
+            << "Could not open pointsFile " << pointFileName
+            << exit(FatalError);
+    }
+}
+
+
+void Foam::CV2D::insertGrid()
+{
+    Info<< "insertInitialGrid: ";
+
+    startOfInternalPoints_ = number_of_vertices();
+    label nVert = startOfInternalPoints_;
+
+    scalar x0 = qSurf_.bb().min().x();
+    scalar xR = qSurf_.bb().max().x() - x0;
+    int ni = int(xR/controls_.minCellSize) + 1;
+    scalar deltax = xR/ni;
+
+    scalar y0 = qSurf_.bb().min().y();
+    scalar yR = qSurf_.bb().max().y() - y0;
+    int nj = int(yR/controls_.minCellSize) + 1;
+    scalar deltay = yR/nj;
+
+    Random rndGen(1321);
+    scalar pert = controls_.randomPurturbation*min(deltax, deltay);
+
+    for (int i=0; i<ni; i++)
+    {
+        for (int j=0; j<nj; j++)
+        {
+            point p(x0 + i*deltax, y0 + j*deltay, 0);
+
+            if (controls_.randomiseInitialGrid)
+            {
+                p.x() += pert*(rndGen.scalar01() - 0.5);
+                p.y() += pert*(rndGen.scalar01() - 0.5);
+            }
+
+            if (qSurf_.wellInside(p, 0.5*controls_.minCellSize2))
+            {
+                insert(Point(p.x(), p.y()))->index() = nVert++;
+            }
+        }
+    }
+
+    Info<< nVert << " vertices inserted" << endl;
+
+    if (controls_.writeInitialTriangulation)
+    {
+        // 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 (controls_.insertSurfaceNearestPointPairs)
+    {
+        insertSurfaceNearestPointPairs();
+    }
+
+    if (controls_.writeNearestTriangulation)
+    {
+        writeFaces("near_allFaces.obj", false);
+        writeFaces("near_faces.obj", true);
+        writeTriangles("near_triangles.obj", true);
+    }
+
+    // Insertion of point-pais for near-points may cause protrusions
+    // so insertBoundaryConformPointPairs must be executed last
+    if (controls_.insertSurfaceNearPointPairs)
+    {
+        insertSurfaceNearPointPairs();
+    }
+
+    startOfBoundaryConformPointPairs_ = number_of_vertices();
+}
+
+
+void Foam::CV2D::boundaryConform()
+{
+    if (!controls_.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<=controls_.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;
+}
+
+
+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)
+{
+    Info<< "newPointsFromVertices: ";
+
+    const vectorField& faceNormals = qSurf_.faceNormals();
+
+    // Initialise the total displacement and its distance for writing out
+    vector2D totalDisp = vector2D::zero;
+    scalar totalDist = 0;
+
+    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 (controls_.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 (controls_.nearWallAlignedDist > 0)
+            {
+                pointIndexHit pHit = qSurf_.tree().findNearest
+                (
+                    toPoint3D(defVert0),
+                    controls_.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 (controls_.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 (controls_.writeFinalTriangulation)
+    {
+        writeFaces("allFaces.obj", false);
+        writeFaces("faces.obj", true);
+        writeTriangles("allTriangles.obj", false);
+        writeTriangles("triangles.obj", true);
+        writePatch("patch.pch");
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/CV2D.H

@@ -0,0 +1,514 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    CV2D
+
+Description
+    Conformal-Voronoi 2D automatic mesher with grid or read initial points
+    and point position relaxation with optional "squarification".
+
+    There are a substantial number of options to this mesher read from 
+    CV2DMesherDict file e.g.:
+
+    // Min cell size used in tolerances when inserting points for
+    // boundary conforming.
+    // Also used to as the grid spacing usind in insertGrid.
+    minCellSize  0.05;
+
+    // Feature angle used to inser feature points
+    // 0 = all features, 180 = no features
+    featureAngle 45;
+
+    // Maximum quadrant angle allowed at a concave corner before
+    // additional "mitering" lines are added
+    maxQuadAngle 110;
+
+    // Should the mesh be square-dominated or of unbiased hexagons
+    squares      yes;
+
+    // Near-wall region where cells are aligned with the wall specified as a 
+    // number of cell layers
+    nearWallAlignedDist 3;
+
+    // Chose if the cell orientation should relax during the iterations
+    //  or remain fixed to the x-y directions
+    relaxOrientation    no;
+
+    // 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;
+
+    // Choose if to randomise the initial grid created by insertGrid.
+    randomiseInitialGrid yes;
+
+    // Perturbation fraction, 1 = cell-size.
+    randomPurturbation   0.1;
+
+    // Number of relaxation iterations.
+    nIterations      5;
+
+    // Relaxation factor at the start of the iteration sequence.
+    // 0.5 is a sensible maximum and < 0.2 converges better.
+    relaxationFactorStart 0.8;
+
+    // Relaxation factor at the end of the iteration sequence.
+    // Should be <= relaxationFactorStart
+    relaxationFactorEnd   0;
+
+    writeInitialTriangulation no;
+    writeFeatureTriangulation no;
+    writeNearestTriangulation no;
+    writeInsertedPointPairs   no;
+    writeFinalTriangulation   yes;
+
+    // Maximum number of iterations used in boundaryConform.
+    maxBoundaryConformingIter 5;
+
+    minEdgeLenCoeff           0.5;
+    maxNotchLenCoeff          0.3;
+    minNearPointDistCoeff     0.25;
+    ppDistCoeff               0.05;
+
+SourceFiles
+    CGALTriangulation2Ddefs.H
+    indexedVertex.H
+    indexedFace.H
+    CV2DI.H
+    CV2D.C
+    CV2DIO.C
+    tolerances.C
+    controls.C
+    insertFeaturePoints.C
+    insertSurfaceNearestPointPairs.C
+    insertSurfaceNearPointPairs.C
+    insertBoundaryConformPointPairs.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CV2D_H
+#define CV2D_H
+
+#define CGAL_INEXACT
+#define CGAL_HIERARCHY
+
+#include "CGALTriangulation2Ddefs.H"
+
+#include "querySurface.H"
+#include "point2DFieldFwd.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class CV2D Declaration
+\*---------------------------------------------------------------------------*/
+
+class CV2D
+:
+    public HTriangulation
+{
+public:
+
+    class controls
+    {
+    public:
+
+        //- Minimum cell size below which protusions through the surface are
+        //  not split
+        scalar minCellSize;
+
+        //- Square of minCellSize
+        scalar minCellSize2;
+
+        //- The feature angle used to select corners to be
+        //  explicitly represented in the mesh.
+        //  0 = all features, 180 = no features
+        scalar featAngle;
+
+        //- Maximum quadrant angle allowed at a concave corner before
+        //  additional "mitering" lines are added
+        scalar maxQuadAngle;
+
+        //- Should the mesh be square-dominated or of unbiased hexagons
+        Switch squares;
+
+        //- Near-wall region where cells are aligned with the wall
+        scalar nearWallAlignedDist;
+
+        //- Square of nearWallAlignedDist
+        scalar nearWallAlignedDist2;
+
+        //- Chose if the cell orientation should relax during the iterations
+        //  or remain fixed to the x-y directions
+        Switch relaxOrientation;
+
+        //- Insert near-boundary point mirror or point-pairs
+        Switch insertSurfaceNearestPointPairs;
+
+        //- Mirror near-boundary points rather than insert point-pairs
+        Switch mirrorPoints;
+
+        //- Insert point-pairs vor dual-cell vertices very near the surface
+        Switch insertSurfaceNearPointPairs;
+
+        Switch writeInitialTriangulation;
+        Switch writeFeatureTriangulation;
+        Switch writeNearestTriangulation;
+        Switch writeInsertedPointPairs;
+        Switch writeFinalTriangulation;
+
+        Switch randomiseInitialGrid;
+        scalar randomPurturbation;
+
+        label maxBoundaryConformingIter;
+
+        //- Relaxation factor at the start of the iteration
+        scalar relaxationFactorStart;
+
+        //- Relaxation factor at the end of the iteration
+        scalar relaxationFactorEnd;
+
+        controls(const dictionary& controlDict);
+    };
+
+
+    class tolerances
+    {
+    public:
+
+        //- Maximum cartesian span of the geometry
+        scalar span;
+
+        //- Square of span
+        scalar span2;
+
+        //- Minumum edge-length of the cell size below which protusions
+        //  through the surface are not split
+        scalar minEdgeLen;
+
+        //- Square of minEdgeLen
+        scalar minEdgeLen2;
+
+        //- Maximum notch size below which protusions into the surface are
+        //  not filled
+        scalar maxNotchLen;
+
+        //- Square of maxNotchLen
+        scalar maxNotchLen2;
+
+        //- The minimum distance alowed between a dual-cell vertex
+        // and the surface before a point-pair is introduced
+        scalar minNearPointDist;
+
+        //- Square of minNearPoint
+        scalar minNearPointDist2;
+
+        //- Distance between boundary conforming point-pairs
+        scalar ppDist;
+
+        //- Square of ppDist
+        scalar ppDist2;
+
+        tolerances
+        (
+            const dictionary& controlDict,
+            scalar minCellSize,
+            const boundBox&
+        );
+    };
+
+
+private:
+
+    // Private data
+
+        //- The surface to mesh
+        const querySurface& qSurf_;
+
+        //- Meshing controls
+        controls controls_;
+
+        //- Meshing tolerances
+        tolerances tols_;
+
+        //- 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_;
+
+        //- 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
+        );
+
+        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();
+
+        //- Insert point groups at the feature points.
+        void insertFeaturePoints();
+
+        //- 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<point2D>& nearSurfacePoints,
+            const DynamicList<point2D>& surfacePoints,
+            const DynamicList<label>& surfaceTris,
+            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 duak-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
+        );
+
+        //- 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);
+
+
+public:
+
+    // Constructors
+
+        //- Construct for given surface
+        CV2D(const dictionary& controlDict, const querySurface& qSurf);
+
+
+    // Destructor
+
+        ~CV2D();
+
+
+    // Member Functions
+
+        // Access
+
+            const controls& meshingControls() const
+            {
+                return controls_;
+            }
+
+
+        // Conversion functions between point2D, point and Point
+
+            inline const point2D& toPoint2D(const point&) 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 point2D& 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(const scalar relaxation);
+
+
+        // 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) 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
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/CV2DI.H

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline Foam::label Foam::CV2D::insertPoint
+(
+    const point2D& p,
+    const label type
+)
+{
+    uint nVert = number_of_vertices();
+
+    Vertex_handle vh = insert(toPoint(p));
+
+    if (nVert == number_of_vertices())
+    {
+        WarningIn("Foam::CV2D::insertPoint")
+            << "Failed to insert point " << p << endl;
+    }
+    else
+    {
+        vh->index() = nVert;
+        vh->type() = type;
+    }
+
+    return vh->index();
+}
+
+
+inline bool Foam::CV2D::insertMirrorPoint
+(
+    const point2D& nearSurfPt,
+    const point2D& surfPt
+)
+{
+    point2D mirrorPoint(2*surfPt - nearSurfPt);
+
+    if (qSurf_.outside(toPoint3D(mirrorPoint)))
+    {
+        insertPoint(mirrorPoint, Vb::MIRROR_POINT);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+inline void Foam::CV2D::insertPointPair
+(
+    const scalar ppDist,
+    const point2D& surfPt,
+    const vector2D& n
+)
+{
+    vector2D ppDistn = ppDist*n;
+
+    label master = insertPoint
+    (
+        surfPt - ppDistn,
+        number_of_vertices() + 1
+    );
+    
+    insertPoint(surfPt + ppDistn, master);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::point2D& Foam::CV2D::toPoint2D(const point& p) const
+{
+    return reinterpret_cast<const point2D&>(p);
+}
+
+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<point2DFromPoint>(P);
+}
+
+inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const
+{
+    return reinterpret_cast<PointFromPoint2D>(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 point2D& p)
+{
+    vh->set_point(toPoint(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)
+    );
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/CV2DIO.C

@@ -0,0 +1,285 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::CV2D::writePoints(const fileName& fName, bool internalOnly) const
+{
+    Info<< "Writing points to " << fName << nl << endl;
+    OFstream str(fName);
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!internalOnly || vit->internalOrBoundaryPoint())
+        {
+            meshTools::writeOBJ(str, toPoint3D(vit->point()));
+        }
+    }
+}
+
+
+void Foam::CV2D::writeTriangles(const fileName& fName, bool internalOnly) const
+{
+    Info<< "Writing triangles to " << fName << nl << endl;
+    OFstream str(fName);
+
+    labelList vertexMap(number_of_vertices());
+    label verti = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!internalOnly || !vit->farPoint())
+        {
+            vertexMap[vit->index()] = verti++;
+            meshTools::writeOBJ(str, toPoint3D(vit->point()));
+        }
+    }
+
+    for
+    (
+        Triangulation::Finite_faces_iterator fit = finite_faces_begin();
+        fit != finite_faces_end();
+        ++fit
+    )
+    {
+        if
+        (
+            !internalOnly 
+         || (
+                fit->vertex(0)->internalOrBoundaryPoint()
+             || fit->vertex(1)->internalOrBoundaryPoint()
+             || fit->vertex(2)->internalOrBoundaryPoint()
+            )
+        )
+        {
+            str << "f " << vertexMap[fit->vertex(0)->index()] + 1
+                << ' ' << vertexMap[fit->vertex(1)->index()] + 1
+                << ' ' << vertexMap[fit->vertex(2)->index()] + 1
+                << nl;
+        }
+    }
+}
+
+
+void Foam::CV2D::writeFaces(const fileName& fName, bool internalOnly) const
+{
+    Info<< "Writing dual faces to " << fName << nl << endl;
+    OFstream str(fName);
+
+    label dualVerti = 0;
+
+    for
+    (
+        Triangulation::Finite_faces_iterator fit = finite_faces_begin();
+        fit != finite_faces_end();
+        ++fit
+    )
+    {
+        if
+        (
+            !internalOnly
+         || (
+                fit->vertex(0)->internalOrBoundaryPoint()
+             || fit->vertex(1)->internalOrBoundaryPoint()
+             || fit->vertex(2)->internalOrBoundaryPoint()
+            )
+        )
+        {
+            fit->faceIndex() = dualVerti++;
+            meshTools::writeOBJ(str, toPoint3D(circumcenter(fit)));
+        }
+        else
+        {
+            fit->faceIndex() = -1;
+        }
+    }
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!internalOnly || vit->internalOrBoundaryPoint())
+        {
+            Face_circulator fcStart = incident_faces(vit);
+            Face_circulator fc = fcStart;
+
+            str<< 'f';
+
+            do
+            {
+                if (!is_infinite(fc))
+                {
+                    if (fc->faceIndex() < 0)
+                    {
+                        FatalErrorIn
+                        (
+                            "Foam::CV2D::writeFaces"
+                            "(const fileName& fName, bool internalOnly)"
+                        )<< "Dual face uses vertex defined by a triangle"
+                            " defined by an external point"
+                            << exit(FatalError);
+                    }
+
+                    str<< ' ' << fc->faceIndex() + 1;
+                }
+            } while (++fc != fcStart);
+
+            str<< nl;
+        }
+    }
+}
+
+
+void Foam::CV2D::calcDual(point2DField& dualPoints, faceList& dualFaces) 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);
+
+
+    // 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;
+
+    calcDual(dual2DPoints, dualFaces);
+    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;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/CV2DMesher.C

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 1991-2007 OpenCFD Ltd.
+    \\/      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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    CV2DMesher
+
+Description
+    Conformal-Voronoi 2D automatic mesher with grid or read initial points
+    and point position relaxation with optional "squarification".
+
+\*---------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+#include "argList.H"
+#include "IFstream.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.clear();
+    argList::validArgs.append("surface");
+    argList::validOptions.insert("pointsFile", "<filename>");
+
+    argList args(argc, argv);
+
+    // Read control dictionary
+    // ~~~~~~~~~~~~~~~~~~~~~~~
+    dictionary controlDict(IFstream(args.executable() + "Dict")());
+
+    label nIterations(readLabel(controlDict.lookup("nIterations")));
+
+
+    // Read the surface to conform to
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    querySurface surf(args.args()[1]);
+    surf.writeTreeOBJ();
+
+    Info<< nl
+        << "Read surface with " << surf.size() << " triangles from file "
+        << args.args()[1] << nl << endl;
+
+
+    // Read and triangulation
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    CV2D mesh(controlDict, surf);
+    if (args.options().found("pointsFile"))
+    {
+        fileName pointFileName(IStringStream(args.options()["pointsFile"])());
+        mesh.insertPoints(pointFileName);
+        mesh.insertSurfacePointPairs();
+        mesh.boundaryConform();
+    }
+    else
+    {
+        mesh.insertGrid();
+        mesh.insertSurfacePointPairs();
+        mesh.boundaryConform();
+    }
+
+    for (int iter=1; iter<=nIterations; iter++)
+    {
+        Info<< nl
+            << "Relaxation iteration " << iter << nl
+            << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
+
+        scalar relax =
+           mesh.meshingControls().relaxationFactorStart
+         +
+         (
+            mesh.meshingControls().relaxationFactorEnd
+          - mesh.meshingControls().relaxationFactorStart
+         )
+         *scalar(iter)/scalar(nIterations);
+
+        mesh.newPoints(relax);
+        mesh.removeSurfacePointPairs();
+        mesh.insertSurfacePointPairs();
+        mesh.boundaryConform();
+    }
+
+    mesh.write();
+
+    Info<< nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/Make/files

@@ -0,0 +1,14 @@
+#include CGAL_FILES
+
+querySurface.C
+CV2D.C
+controls.C
+tolerances.C
+insertFeaturePoints.C
+insertSurfaceNearestPointPairs.C
+insertSurfaceNearPointPairs.C
+insertBoundaryConformPointPairs.C
+CV2DIO.C
+CV2DMesher.C
+
+EXE = $(FOAM_APPBIN)/CV2DMesher

applications/utilities/mesh/generation/CV2DMesher/Make/options

@@ -0,0 +1,16 @@
+//EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_NDEBUG = -DNDEBUG
+
+include $(GENERAL_RULES)/CGAL
+FFLAGS = -DCGAL_FILES='"${CGAL_PATH}/CGAL/files"'
+
+EXE_INC = \
+    ${EXE_NDEBUG} \
+    ${CGAL_INC} \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -ltriSurface
+

applications/utilities/mesh/generation/CV2DMesher/controls.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CV2D::controls::controls(const dictionary& controlDict)
+:
+    minCellSize(readScalar(controlDict.lookup("minCellSize"))),
+    minCellSize2(Foam::sqr(minCellSize)),
+
+    featAngle(readScalar(controlDict.lookup("featureAngle"))),
+    maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))),
+    squares(controlDict.lookup("squares")),
+
+    nearWallAlignedDist
+    (
+        readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize
+    ),
+    nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)),
+
+    relaxOrientation(controlDict.lookup("relaxOrientation")),
+
+    insertSurfaceNearestPointPairs
+    (
+        controlDict.lookup("insertSurfaceNearestPointPairs")
+    ),
+    mirrorPoints(controlDict.lookup("mirrorPoints")),
+    insertSurfaceNearPointPairs
+    (
+        controlDict.lookup("insertSurfaceNearPointPairs")
+    ),
+    writeInitialTriangulation(controlDict.lookup("writeInitialTriangulation")),
+    writeFeatureTriangulation(controlDict.lookup("writeFeatureTriangulation")),
+    writeNearestTriangulation(controlDict.lookup("writeNearestTriangulation")),
+    writeInsertedPointPairs(controlDict.lookup("writeInsertedPointPairs")),
+    writeFinalTriangulation(controlDict.lookup("writeFinalTriangulation")),
+    randomiseInitialGrid(controlDict.lookup("randomiseInitialGrid")),
+    randomPurturbation(readScalar(controlDict.lookup("randomPurturbation"))),
+    maxBoundaryConformingIter
+    (
+        readLabel(controlDict.lookup("maxBoundaryConformingIter"))
+    ),
+    relaxationFactorStart
+    (
+        readScalar(controlDict.lookup("relaxationFactorStart"))
+    ),
+    relaxationFactorEnd(readScalar(controlDict.lookup("relaxationFactorEnd")))
+{}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/indexedFace.H

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    indexedFace
+
+Description
+    An indexed form of CGAL::Triangulation_face_base_2<K> used to keep
+    track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedFace_H
+#define indexedFace_H
+
+#include <CGAL/Triangulation_2.h>
+#include "indexedVertex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class indexedFace Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Gt, class Fb=CGAL::Triangulation_face_base_2<Gt> >
+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<TDS2>::Other Fb2;
+        typedef indexedFace<Gt, Fb2>  Other;
+    };
+
+
+    indexedFace()
+    :
+        Fb(),
+        index_(CHANGED)
+    {}
+
+    indexedFace(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2)
+    :
+        Fb(v0, v1, v2),
+        index_(CHANGED)
+    {}
+
+    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)
+    {}
+
+
+    void set_vertex(int i, Vertex_handle v)
+    {
+        index_ = CHANGED;
+        Fb::set_vertex(i, v);
+    }
+
+    void set_vertices()
+    {
+        index_ = CHANGED;
+        Fb::set_vertices();
+    }
+
+    void set_vertices(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2)
+    {
+        index_ = CHANGED;
+        Fb::set_vertices(v0, v1, v2);
+    }
+
+
+    int& faceIndex()
+    {
+        return index_;
+    }
+
+    int faceIndex() const
+    {
+        return index_;
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/indexedVertex.H

@@ -0,0 +1,261 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    indexedVertex
+
+Description
+    An indexed form of CGAL::Triangulation_vertex_base_2<K> used to keep
+    track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedVertex_H
+#define indexedVertex_H
+
+#include <CGAL/Triangulation_2.h>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class indexedVertex Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Gt, class Vb=CGAL::Triangulation_vertex_base_2<Gt> >
+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<typename TDS2>
+    struct Rebind_TDS
+    {
+        typedef typename Vb::template Rebind_TDS<TDS2>::Other    Vb2;
+        typedef indexedVertex<Gt,Vb2>                           Other;
+    };
+
+
+    indexedVertex()
+    :
+        Vb(),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT)
+    {}
+
+    indexedVertex(const Point& p)
+    :
+        Vb(p),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT)
+    {}
+
+    indexedVertex(const Point& p, Face_handle f)
+    :
+        Vb(f, p),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT)
+    {}
+
+    indexedVertex(Face_handle f)
+    :
+        Vb(f),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT)
+    {}
+
+
+    int& index()
+    {
+        return index_;
+    }
+
+    int index() const
+    {
+        return index_;
+    }
+
+
+    int& type()
+    {
+        return type_;
+    }
+
+    int type() const
+    {
+        return type_;
+    }
+
+
+    //- Is point a far-point
+    inline bool farPoint() const
+    {
+        return type_ == FAR_POINT;
+    }
+
+    //- Is point internal, i.e. not on boundary
+    inline bool internalPoint() const
+    {
+        return type_ <= INTERNAL_POINT;
+    }
+
+    //- Is point internal and near the boundary
+    inline bool nearBoundary() const
+    {
+        return type_ == NEAR_BOUNDARY_POINT;
+    }
+
+    //- Set the point to be near the boundary
+    inline void setNearBoundary()
+    {
+        type_ = NEAR_BOUNDARY_POINT;
+    }
+
+    //- Is point a mirror point
+    inline bool mirrorPoint() const
+    {
+        return type_ == MIRROR_POINT;
+    }
+
+    //- Either master or slave of pointPair.
+    inline bool pairPoint() const
+    {
+        return type_ >= 0;
+    }
+
+    //- Master of a pointPair is the lowest numbered one.
+    inline bool ppMaster() const
+    {
+        if (type_ > index_)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    //- Slave of a pointPair is the highest numbered one.
+    inline bool ppSlave() const
+    {
+        if (type_ >= 0 && type_ < index_)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    //- Either original internal point or master of pointPair.
+    inline bool internalOrBoundaryPoint() const
+    {
+        return internalPoint() || ppMaster();
+    }
+
+    //- Is point near the boundary or part of the boundary definition
+    inline bool nearOrOnBoundary() const
+    {
+        return pairPoint() || mirrorPoint() || nearBoundary();
+    }
+
+    //- Do the two given vertices consitute a boundary point-pair
+    inline friend bool pointPair
+    (
+        const indexedVertex& v0,
+        const indexedVertex& v1
+    )
+    {
+        return v0.index_ == v1.type_ || v1.index_ == v0.type_;
+    }
+
+    //- Do the three given vertices consitute a boundary triangle
+    inline friend bool 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));
+    }
+
+    //- Do the three given vertices consitute an outside triangle
+    inline friend bool outsideTriangle
+    (
+        const indexedVertex& v0,
+        const indexedVertex& v1,
+        const indexedVertex& v2
+    )
+    {
+        return (v0.farPoint() || v0.ppSlave())
+            || (v1.farPoint() || v1.ppSlave())
+            || (v2.farPoint() || v2.ppSlave());
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/insertBoundaryConformPointPairs.C

@@ -0,0 +1,290 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV2D::insertPointPair
+(
+    Triangulation::Finite_vertices_iterator& vit,
+    const point2D& p,
+    const label trii
+)
+{
+    if
+    (
+        !controls_.mirrorPoints
+     || !insertMirrorPoint(toPoint2D(vit->point()), p)
+    )
+    {
+        insertPointPair
+        (
+            tols_.ppDist,
+            p,
+            toPoint2D(qSurf_.faceNormals()[trii])
+        );
+    }
+
+    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;
+    scalar interDist2 = 0;
+
+    Face_circulator fcStart = incident_faces(vit);
+    Face_circulator fc = fcStart;
+    label vi = 0;
+
+    do
+    {
+        if (!is_infinite(fc))
+        {
+            pointIndexHit pHit = qSurf_.tree().findLine
+            (
+                toPoint3D(defVert),
+                toPoint3D(vertices[vi])
+            );
+
+            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 = tols_.maxNotchLen2;
+                    }
+
+                    if (dist2 > mps2)
+                    {
+                        found = true;
+                        interPoint = toPoint2D(pHit.hitPoint());
+                        interTri = pHit.index();
+                        interDist2 = dist2;
+                    }
+                }
+            }
+
+            vi++;
+        }
+    } while (++fc != fcStart);
+
+    if (found)
+    {
+        insertPointPair(vit, interPoint, interTri);
+        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 = tols_.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<nVerts; vi++)
+            {
+                vector2D edge(vertices[vi+1] - vertices[vi]);
+                perimeter += mag(edge);
+                vector2D otherEdge = defVert - vertices[vi];
+                areaT2 += mag(edge.x()*otherEdge.y() - edge.y()*otherEdge.x());
+            }
+
+            // If the dual-cell is very small reject refinement
+            if (areaT2 < tols_.minEdgeLen2) continue;
+
+            // Estimate the cell width
+            scalar cellWidth = areaT2/perimeter;
+
+
+            // Check dimensions of dual-cell
+            /*
+            // Quick rejection of dual-cell refinement based on it's perimeter
+            if (perimeter < 2*tols_.minCellSize) continue;
+
+            // Also check the area of the cell and reject refinement
+            // if it is less than that allowed
+            if (areaT2 < tols_.minCellSize2) continue;
+
+            // Estimate the cell width and reject refinement if it is less than
+            // that allowed
+            if (cellWidth < 0.5*tols_.minEdgeLen) continue;
+            */
+
+            if
+            (
+                perimeter > 2*controls_.minCellSize
+             && areaT2 > controls_.minCellSize2
+             && cellWidth > 0.5*tols_.minEdgeLen
+            )
+            {
+                maxProtSize2 = 0.25*tols_.maxNotchLen2;
+            }
+        }
+
+        if (insertPointPairAtIntersection(vit, defVert, vertices, maxProtSize2))
+        {
+            nIntersections++;
+        }
+    }
+
+    return nIntersections;
+}
+
+
+void Foam::CV2D::markNearBoundaryPoints()
+{
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        if (vit->internalPoint())
+        {
+            point vert(toPoint3D(vit->point()));
+
+            pointIndexHit pHit =
+                qSurf_.tree().findNearest(vert, 4*controls_.minCellSize2);
+
+            if (pHit.hit())
+            {
+                vit->setNearBoundary();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/insertFeaturePoints.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+#include "plane.H"
+#include "triSurfaceTools.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Create feature points/edges by creating a triplet in the corner.
+// (this triplet will have as its circumcentre the feature)
+void Foam::CV2D::insertFeaturePoints()
+{
+    labelList featEdges(qSurf_.extractFeatures2D(controls_.featAngle));
+
+    const pointField& localPts = qSurf_.localPoints();
+
+    forAll(featEdges, i)
+    {
+        label edgeI = featEdges[i];
+        const edge& featEdge = qSurf_.edges()[edgeI];
+
+        // Get the feature point as the mid-point of the edge and convert to 2D
+        point2D featPt = toPoint2D(featEdge.centre(qSurf_.localPoints()));
+
+        // Pick up the two faces adjacent to the feature edge
+        const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
+
+        label faceA = eFaces[0];
+        vector2D nA = toPoint2D(qSurf_.faceNormals()[faceA]);
+
+        label faceB = eFaces[1];
+        vector2D nB = toPoint2D(qSurf_.faceNormals()[faceB]);
+
+        // Intersect planes parallel to faceA and faceB offset by ppDist.
+        plane planeA(toPoint3D(featPt - tols_.ppDist*nA), toPoint3D(nA));
+        plane planeB(toPoint3D(featPt - tols_.ppDist*nB), toPoint3D(nB));
+        plane::ray interLine(planeA.planeIntersect(planeB));
+
+        // The reference point is where this line intersects the z_ plane
+        point2D refPt = toPoint2D
+        (
+            interLine.refPoint()
+          + ((z_ - interLine.refPoint().z())/interLine.dir().z())
+           *interLine.dir()
+        );
+
+        point2D faceAVert = toPoint2D
+        (
+            localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)]
+        );
+
+        // Determine convex or concave angle
+        if (((faceAVert - featPt) & nB) < 0)
+        {
+            // Convex. So refPt will be inside domain and hence a master point
+
+            // Insert the master point refering the the first slave
+            label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
+
+            // Insert the slave points by reflecting refPt in both faces.
+            // with each slave refering to the master
+
+            point2D reflectedA = refPt + 2*((featPt - refPt) & nA)*nA;
+            insertPoint(reflectedA, masterPtIndex);
+
+            point2D reflectedB = refPt + 2*((featPt - refPt) & nB)*nB;
+            insertPoint(reflectedB, masterPtIndex);
+        }
+        else
+        {
+            // Concave. master and reflected points inside the domain.
+            // Generate reflected master to be outside.
+            point2D reflMasterPt = refPt + 2*(featPt - refPt);
+
+            // Reflect refPt in both faces.
+            point2D reflectedA = 
+                reflMasterPt + 2*((featPt - reflMasterPt) & nA)*nA;
+
+            point2D reflectedB =
+                reflMasterPt + 2*((featPt - reflMasterPt) & nB)*nB;
+
+            // Total angle around the concave feature
+//             scalar totalAngle = 
+//                 180*(2.0*mathematicalConstant::pi - acos(mag(nA & nB)))
+//                /mathematicalConstant::pi;
+
+            scalar totalAngle = 
+                180*(mathematicalConstant::pi + acos(mag(nA & nB)))
+               /mathematicalConstant::pi;
+
+            // Number of quadrants the angle should be split into
+            int nQuads = int(totalAngle/controls_.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 = number_of_vertices() + 2 + nAddPoints;
+
+            // Master A is inside.
+            label reflectedAI = insertPoint(reflectedA, reflectedMaster);
+
+            // Master B is inside.
+            insertPoint(reflectedB, reflectedMaster);
+
+            if (nAddPoints == 1)
+            {
+                // One additinal point is the reflection of the slave point,
+                // i.e. the original reference point
+                insertPoint(refPt, reflectedMaster);
+            }
+            else if (nAddPoints == 2)
+            {
+                point2D reflectedAa = refPt - ((featPt - reflMasterPt) & nB)*nB;
+                insertPoint(reflectedAa, reflectedMaster);
+
+                point2D reflectedBb = refPt - ((featPt - reflMasterPt) & nA)*nA;
+                insertPoint(reflectedBb, reflectedMaster);
+            }
+
+            // Slave is outside.
+            insertPoint(reflMasterPt, reflectedAI);
+        }
+    }
+
+    if (controls_.writeFeatureTriangulation)
+    {
+        writePoints("feat_allPoints.obj", false);
+        writeFaces("feat_allFaces.obj", false);
+        writeFaces("feat_faces.obj", true);
+        writeTriangles("feat_triangles.obj", true);
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/insertSurfaceNearPointPairs.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV2D::insertSurfaceNearPointPairs()
+{
+    Info<< "insertSurfaceNearPointPairs: ";
+
+    label nNearPoints = 0;
+
+    for
+    (
+        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        eit++
+    )
+    {
+        Vertex_handle v0h = eit->first->vertex(cw(eit->second));
+        Vertex_handle v1h = eit->first->vertex(ccw(eit->second));
+
+        if (v0h->ppMaster() && v1h->ppMaster())
+        {
+            point2DFromPoint v0(toPoint2D(v0h->point()));
+            point2DFromPoint v1(toPoint2D(v1h->point()));
+
+            // Check that the two triangle vertices are further apart than the
+            // minimum cell size
+            if (magSqr(v1 - v0) > controls_.minCellSize2)
+            {
+                point2D e0(toPoint2D(circumcenter(eit->first)));
+
+                point2D e1
+                (
+                    toPoint2D(circumcenter(eit->first->neighbor(eit->second)))
+                );
+
+                // Calculate the length^2 of the edge normal to the surface
+                scalar edgeLen2 = magSqr(e0 - e1);
+
+                if (edgeLen2 < tols_.minNearPointDist2)
+                {
+                    pointIndexHit pHit = qSurf_.tree().findNearest
+                    (
+                        toPoint3D(e0),
+                        tols_.minEdgeLen2
+                    );
+
+                    if (pHit.hit())
+                    {
+                        insertPointPair
+                        (
+                            tols_.ppDist,
+                            toPoint2D(pHit.hitPoint()),
+                            toPoint2D(qSurf_.faceNormals()[pHit.index()])
+                        );
+
+                        nNearPoints++;
+
+                        // Correct the edge iterator for the change in the 
+                        // number od edges following the point-pair insertion
+                        eit = Finite_edges_iterator
+                        (
+                            finite_edges_end().base(),
+                            eit.predicate(),
+                            eit.base()
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    Info<< nNearPoints << " point-pairs inserted" << endl;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/insertSurfaceNearestPointPairs.C

@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::CV2D::dualCellSurfaceIntersection
+(
+    const Triangulation::Finite_vertices_iterator& vit
+) const
+{
+    Triangulation::Edge_circulator ecStart = incident_edges(vit);
+    Triangulation::Edge_circulator ec = ecStart;
+
+    do
+    {
+        if (!is_infinite(ec))
+        {
+            point e0 = toPoint3D(circumcenter(ec->first));
+
+            // If edge end is outside bounding box then edge cuts boundary
+            if (!qSurf_.bb().contains(e0))
+            {
+                return true;
+            }
+
+            point e1 = toPoint3D(circumcenter(ec->first->neighbor(ec->second)));
+
+            // If other edge end is ouside bounding box then edge cuts boundary
+            if (!qSurf_.bb().contains(e1))
+            {
+                return true;
+            }
+
+            if (magSqr(e1 - e0) > tols_.minEdgeLen2)
+            {
+                pointIndexHit pHit = qSurf_.tree().findLineAny(e0, e1);
+
+                if (pHit.hit())
+                {
+                    return true;
+                }
+            }
+        }
+
+    } while (++ec != ecStart);
+
+    return false;
+}
+
+
+void Foam::CV2D::insertPointPairs
+(
+    const DynamicList<point2D>& nearSurfacePoints,
+    const DynamicList<point2D>& surfacePoints,
+    const DynamicList<label>& surfaceTris,
+    const fileName fName
+)
+{
+    if (controls_.mirrorPoints)
+    {
+        forAll(surfacePoints, ppi)
+        {
+            insertMirrorPoint
+            (
+                nearSurfacePoints[ppi],
+                surfacePoints[ppi]
+            );
+        }
+    }
+    else
+    {
+        forAll(surfacePoints, ppi)
+        {
+            insertPointPair
+            (
+                tols_.ppDist,
+                surfacePoints[ppi],
+                toPoint2D(qSurf_.faceNormals()[surfaceTris[ppi]])
+            );
+        }
+    }
+
+    Info<< surfacePoints.size() << " point-pairs inserted" << endl;
+
+    if (controls_.writeInsertedPointPairs)
+    {
+        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<point2D> nearSurfacePoints(nSurfacePointsEst);
+    DynamicList<point2D> surfacePoints(nSurfacePointsEst);
+    DynamicList<label> surfaceTris(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 = qSurf_.tree().findNearest
+            (
+                toPoint3D(vert),
+                4*controls_.minCellSize2
+            );
+
+            if (pHit.hit())
+            {
+                vit->setNearBoundary();
+
+                // Reference to the nearest triangle
+                const labelledTri& f = qSurf_[pHit.index()];
+
+                // 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(), 1e-6, 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[pHit.index()][nearLabel];
+                    const labelList& eFaces = edgeFaces[edgeI];
+
+                    if
+                    (
+                        eFaces.size() == 2
+                     && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
+                       < -0.2
+                    )
+                    {
+                        internalFeatureEdge = true;
+                    }
+                }
+
+                if (!internalFeatureEdge && dualCellSurfaceIntersection(vit))
+                {
+                    nearSurfacePoints.append(vert);
+                    surfacePoints.append(toPoint2D(pHit.hitPoint()));
+                    surfaceTris.append(pHit.index());
+                }
+            }
+        }
+    }
+
+    insertPointPairs
+    (
+        nearSurfacePoints,
+        surfacePoints,
+        surfaceTris,
+        "surfaceNearestIntersections.obj"
+    );
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/querySurface.C

@@ -0,0 +1,229 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "querySurface.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::querySurface::querySurface(const fileName& surfaceFileName)
+:
+    triSurface(surfaceFileName),
+    rndGen_(12345),
+    bb_(localPoints()),
+    tree_
+    (
+        treeDataTriSurface(*this),
+        bb_.extend(rndGen_, 1e-3), // slightly randomize bb
+        8,                         // maxLevel
+        4, //10,                   // leafsize
+        10.0 //3.0                 // duplicity
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::querySurface::~querySurface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::querySurface::extractFeatures2D
+(
+    const scalar featAngle
+) const
+{
+    scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
+
+    const labelListList& edgeFaces = this->edgeFaces();
+    const pointField& localPoints = this->localPoints();
+    const edgeList& edges = this->edges();
+    const vectorField& faceNormals = this->faceNormals();
+
+    DynamicList<label> featEdges(edges.size());
+
+    forAll(edgeFaces, edgeI)
+    {
+        const edge& e = edges[edgeI];
+
+        if (magSqr(e.vec(localPoints) & vector(1,1,0)) < SMALL)
+        {
+            const labelList& eFaces = edgeFaces[edgeI];
+
+            if
+            (
+                eFaces.size() == 2
+             && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
+            )
+            {
+                featEdges.append(edgeI);
+            }
+        }
+    }
+
+    return featEdges.shrink();
+}
+
+Foam::indexedOctree<Foam::treeDataTriSurface>::volumeType
+Foam::querySurface::insideOutside
+(
+    const scalar searchSpan2,
+    const point& pt
+) const
+{
+    if (!bb_.contains(pt))
+    {
+        return indexedOctree<treeDataTriSurface>::OUTSIDE;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, searchSpan2);
+
+    if (!pHit.hit())
+    {
+        return tree_.getVolumeType(pt);
+    }
+    else
+    {
+        return indexedOctree<treeDataTriSurface>::MIXED;
+    }
+}
+
+
+// Check if point is inside surface
+bool Foam::querySurface::inside(const point& pt) const
+{
+    if (!bb_.contains(pt))
+    {
+        return false;
+    }
+
+    return
+    (
+        tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::INSIDE
+    );
+}
+
+
+// Check if point is outside surface
+bool Foam::querySurface::outside(const point& pt) const
+{
+    if (!bb_.contains(pt))
+    {
+        return true;
+    }
+
+    return
+    (
+        tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::OUTSIDE
+    );
+}
+
+
+// Check if point is inside surface by at least dist2
+bool Foam::querySurface::wellInside(const point& pt, const scalar dist2) const
+{
+    if (!bb_.contains(pt))
+    {
+        return false;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, dist2);
+
+    if (pHit.hit())
+    {
+        return false;
+    }
+    else
+    {
+        return
+            tree_.getVolumeType(pt)
+         == indexedOctree<treeDataTriSurface>::INSIDE;
+    }
+}
+
+
+// Check if point is outside surface by at least dist2
+bool Foam::querySurface::wellOutside(const point& pt, const scalar dist2) const
+{
+    if (!bb_.contains(pt))
+    {
+        return true;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, dist2);
+
+    if (pHit.hit())
+    {
+        return false;
+    }
+    else
+    {
+        return
+            tree_.getVolumeType(pt)
+         == indexedOctree<treeDataTriSurface>::OUTSIDE;
+    }
+}
+
+
+void Foam::querySurface::writeTreeOBJ() const
+{
+    OFstream str("tree.obj");
+    label vertI = 0;
+
+    const List<indexedOctree<treeDataTriSurface>::node>& nodes = tree_.nodes();
+
+    forAll(nodes, nodeI)
+    {
+        const indexedOctree<treeDataTriSurface>::node& nod = nodes[nodeI];
+
+        const treeBoundBox& bb = nod.bb_;
+
+        const pointField points(bb.points());
+
+        label startVertI = vertI;
+
+        forAll(points, i)
+        {
+            meshTools::writeOBJ(str, points[i]);
+            vertI++;
+        }
+
+        const edgeList edges(treeBoundBox::edges);
+
+        forAll(edges, i)
+        {
+            const edge& e = edges[i];
+
+            str << "l " << e[0]+startVertI+1 << ' ' << e[1]+startVertI+1
+                << nl;
+        }
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/querySurface.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    querySurface
+
+Description
+    Searchable triSurface using an octree to speed-up queries.
+
+SourceFiles
+    querySurface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef querySurface_H
+#define querySurface_H
+
+#include "triSurface.H"
+#include "treeDataTriSurface.H"
+#include "indexedOctree.H"
+#include "Random.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class querySurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class querySurface
+:
+    public triSurface
+{
+    // Private data
+
+        Random rndGen_;
+
+        // Bounding box of surface. Used for relative tolerances.
+        treeBoundBox bb_;
+
+        // Search engine on surface
+        indexedOctree<treeDataTriSurface> tree_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        querySurface(const querySurface&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const querySurface&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given file name of the surface
+        querySurface(const fileName& surfaceFileName);
+
+
+    // Destructor
+
+        ~querySurface();
+
+
+    // Member Functions
+
+        // Access
+
+            const treeBoundBox& bb() const
+            {
+                return bb_;
+            }
+
+            const indexedOctree<treeDataTriSurface>& tree() const
+            {
+                return tree_;
+            }
+
+
+        // Query
+
+            //- Extract feature edges/points(2D)
+            //  using the given feature angle in deg
+            labelList extractFeatures2D(const scalar featAngle) const;
+
+            //- Returns inside, outside or mixed (= close to surface)
+            indexedOctree<Foam::treeDataTriSurface>::volumeType insideOutside
+            (
+                const scalar searchSpan2,
+                const point& pt
+            ) const;
+
+            //- Check if point is inside surface
+            bool inside(const point& pt) const;
+
+            //- Check if point is outside surface
+            bool outside(const point& pt) const;
+
+            //- Check if point is inside surface by at least dist2
+            bool wellInside(const point& pt, const scalar dist2) const;
+
+            //- Check if point is outside surface by at least dist2
+            bool wellOutside(const point& pt, const scalar dist2) const;
+
+
+        // Write
+
+            void writeTreeOBJ() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//#include "querySurfaceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV2DMesher/tolerances.C

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV2D.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CV2D::tolerances::tolerances
+(
+    const dictionary& controlDict,
+    const scalar minCellSize,
+    const boundBox& bb
+)
+:
+    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(controlDict.lookup("minEdgeLenCoeff"))*minCellSize),
+    minEdgeLen2(Foam::sqr(minEdgeLen)),
+
+    maxNotchLen(readScalar(controlDict.lookup("maxNotchLenCoeff"))*minCellSize),
+    maxNotchLen2(Foam::sqr(maxNotchLen)),
+
+    minNearPointDist
+    (
+        readScalar(controlDict.lookup("minNearPointDistCoeff"))*minCellSize
+    ),
+    minNearPointDist2(Foam::sqr(minNearPointDist)),
+
+    ppDist(readScalar(controlDict.lookup("ppDistCoeff"))*minCellSize),
+    ppDist2(Foam::sqr(ppDist))
+{}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/CGALTriangulation3Ddefs.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Typedefs
+    CGALTriangulation3Ddefs
+
+Description
+    CGAL data structures used for 3D Delaunay meshing.
+
+    Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel
+    otherwise the more robust but much less efficient
+    Exact_predicates_exact_constructions will be used.
+
+    Define CGAL_HIERARCHY to use hierarchical Delaunay triangulation which is
+    faster but uses more memory than the standard Delaunay triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CGALTriangulation3Ddefs_H
+#define CGALTriangulation3Ddefs_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CGAL/Delaunay_triangulation_3.h"
+
+#include "indexedVertex.H"
+#include "indexedCell.H"
+
+#ifdef CGAL_INEXACT
+    // Fast kernel using a double as the storage type but the triangulation
+    // may fail
+    #include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
+    typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+#else
+    // Very robust but expensive kernel
+    #include "CGAL/Exact_predicates_exact_constructions_kernel.h"
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+#endif
+
+typedef CGAL::indexedVertex<K>                           Vb;
+typedef CGAL::indexedCell<K>                             Cb;
+
+#ifdef CGAL_HIERARCHY
+    // Data structures for hierarchical Delaunay triangulation which is more
+    // efficient but also uses more storage
+    #include "CGAL/Triangulation_hierarchy_3.h"
+    typedef CGAL::Triangulation_hierarchy_vertex_base_3<Vb>  Vbh;
+    typedef CGAL::Triangulation_data_structure_3<Vbh, Cb>    Tds;
+    typedef CGAL::Delaunay_triangulation_3<K, Tds>           Triangulation;
+    typedef CGAL::Triangulation_hierarchy_3<Triangulation>   HTriangulation;
+#else
+    // Data structures for standard Delaunay triangulation
+    typedef CGAL::Triangulation_data_structure_3<Vb, Cb>     Tds;
+    typedef CGAL::Delaunay_triangulation_3<K, Tds>           Triangulation;
+    typedef Triangulation                                    HTriangulation;
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/CV3D.H

@@ -0,0 +1,460 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    CV3D
+
+Description
+
+SourceFiles
+    CV3DI.H
+    CV3D.C
+    CV3DIO.C
+    controls.C
+    tolerances.C
+    insertFeaturePoints.C
+    insertSurfaceNearestPointPairs.C
+    insertSurfaceNearPointPairs.C
+    insertBoundaryConformPointPairs.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CV3D_H
+#define CV3D_H
+
+//#define CGAL_INEXACT
+#define CGAL_HIERARCHY
+
+#include "CGALTriangulation3Ddefs.H"
+
+#include "querySurface.H"
+#include "IOdictionary.H"
+#include "pointIOField.H"
+#include "argList.H"
+#include "DynamicList.H"
+#include "Switch.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "SortableList.H"
+#include "plane.H"
+#include "triSurfaceTools.H"
+#include "mathematicalConstants.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+/*---------------------------------------------------------------------------*\
+                           Class CV3D Declaration
+\*---------------------------------------------------------------------------*/
+
+class CV3D
+:
+    public HTriangulation
+{
+public:
+
+    class controls
+    {
+    public:
+
+        //- Relaxation factor at the start of the iteration
+        scalar relaxationFactorStart;
+
+        //- Relaxation factor at the end of the iteration
+        scalar relaxationFactorEnd;
+
+        //- Minimum cell size below which protusions through the surface are
+        //  not split
+        scalar minCellSize;
+
+        //- Square of minCellSize
+        scalar minCellSize2;
+
+        //- The feature angle used to select corners to be
+        //  explicitly represented in the mesh.
+        //  0 = no features, 180 = all features
+        scalar includedAngle;
+
+        //- Maximum quadrant angle allowed at a concave corner before
+        //  additional "mitering" lines are added
+        scalar maxQuadAngle;
+
+        //- Should the mesh be square-dominated or of unbiased hexagons
+        Switch squares;
+
+        //- Near-wall region where cells are aligned with the wall
+        scalar nearWallAlignedDist;
+
+        //- Square of nearWallAlignedDist
+        scalar nearWallAlignedDist2;
+
+        //- Chose if the cell orientation should relax during the iterations
+        //  or remain fixed to the x-y directions
+        Switch relaxOrientation;
+
+        //- Insert near-boundary point mirror or point-pairs
+        Switch insertSurfaceNearestPointPairs;
+
+        //- Mirror near-boundary points rather than insert point-pairs
+        Switch mirrorPoints;
+
+        //- Insert point-pairs for dual-cell vertices very near the surface
+        Switch insertSurfaceNearPointPairs;
+
+        Switch writeInitialTriangulation;
+        Switch writeFeatureTriangulation;
+        Switch writeNearestTriangulation;
+        Switch writeInsertedPointPairs;
+        Switch writeFinalTriangulation;
+
+        Switch randomiseInitialGrid;
+        scalar randomPerturbation;
+
+        controls(const IOdictionary& controlDict);
+    };
+
+    class tolerances
+    {
+    public:
+
+        //- Maximum cartesian span of the geometry
+        scalar span;
+
+        //- Square of span
+        scalar span2;
+
+        //- Minumum edge-length of the cell size below which protusions
+        //  through the surface are not split
+        scalar minEdgeLen;
+
+        //- Square of minEdgeLen
+        scalar minEdgeLen2;
+
+        //- Distance between boundary conforming point-pairs
+        scalar ppDist;
+
+        //- Square of ppDist
+        scalar ppDist2;
+
+        //- Maximum distance between adjacent control points on a feature edge
+        //- before they are considered a group
+        scalar edgeGroupSpacing;
+
+        //- Guard distance from a feature point within which an edge control
+        //- point may not be placed
+        scalar featurePointGuard;
+
+        //- Guard distance from a feature edge within which a surface control
+        //- point may not be placed
+        scalar featureEdgeGuard;
+
+        //- minimum and maximum distances along a feature edge allowed between
+        //- pairs of points.  Eventually these should be able to adapt to local
+        //- grading requirements.
+        scalar minEdgeSpacing;
+        scalar maxEdgeSpacing;
+
+        tolerances
+        (
+            const IOdictionary& controlDict,
+            scalar minCellSize,
+            const boundBox&
+        );
+    };
+
+private:
+
+// Private data
+
+        //- The surface to mesh
+        const querySurface& qSurf_;
+
+        //- The time registry of the application
+        const Time& runTime_;
+
+        //- Meshing controls
+        controls controls_;
+
+        //- Meshing tolerances
+        tolerances tols_;
+
+        //- 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 constraining points to be reinserted after a
+        //- triangulation clear
+        List<Vb> featureConstrainingVertices_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        CV3D(const CV3D&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const CV3D&);
+
+        //- Insert point and return it's index
+        inline label insertPoint
+        (
+            const point& pt,
+            const label type
+        );
+
+        inline bool insertMirrorPoint
+        (
+            const point& nearSurfPt,
+            const point& 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 point& surfPt,
+            const vector& n
+        );
+
+        inline void insertVb(const Vb& v);
+
+        inline void movePoint(const Vertex_handle& vh, const point& p);
+
+        //- Create the initial mesh from the bounding-box
+        void insertBoundingBox();
+
+        //- Insert point groups at the feature points.
+        void insertFeaturePoints();
+
+        //- Reinsert stored feature point defining points.
+        void reinsertFeaturePoints();
+
+        //- Reinsert points that have been saved to clear the mesh.
+        void reinsertPoints (const pointField& points);
+
+        //- 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<point>& nearSurfacePoints,
+            const DynamicList<point>& surfacePoints,
+            const DynamicList<label>& surfaceTris,
+            const fileName fName
+        );
+
+        //- Insert edge control groups at the given edge points using the
+        //  normals of the faces attached to the edge.
+        void insertEdgePointGroups
+        (
+            const DynamicList<point>& edgePoints,
+            const DynamicList<label>& edgeLabels,
+            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;
+
+        //- Smooths the points positioned along a feature edge by collecting
+        //- bunches together and deleting excess points.
+        void smoothEdgePositions
+        (
+            DynamicList<point>& edgePoints,
+            DynamicList<label>& edgeLabels
+        ) const;
+
+        void smoothEdge
+        (
+            List<point>& edgePoints,
+            List<scalar>& edgeDistances,
+            const label edgeI
+        ) 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();
+
+        void markNearBoundaryPoints();
+
+        //- Dual calculation
+        void calcDualMesh
+        (
+            pointField& points,
+            faceList& faces,
+            labelList& owner,
+            labelList& neighbour,
+            wordList& patchNames,
+            labelList& patchSizes,
+            labelList& patchStarts
+        );
+
+        void setVertexAlignmentDirections();
+
+        scalar alignmentDistanceWeight(scalar dist) const;
+
+        scalar faceAreaWeight(scalar faceArea) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct for given surface
+        CV3D
+        (
+            const Time& runTime,
+            const querySurface& qSurf
+        );
+
+    // Destructor
+
+        ~CV3D();
+
+
+    // Member Functions
+
+        // Access
+
+            const controls& meshingControls() const
+            {
+                return controls_;
+            }
+
+            // Conversion functions between point and Point
+
+#           ifdef CGAL_INEXACT
+                typedef const point& pointFromPoint;
+                typedef const Point& PointFrompoint;
+#           else
+                typedef point pointFromPoint;
+                typedef Point PointFrompoint;
+#           endif
+
+            inline pointFromPoint topoint(const Point&) const;
+            inline PointFrompoint toPoint(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 pointField& 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 BCC lattice of points.
+            //  Points outside the geometry are ignored.
+            void insertGrid();
+
+
+        // Point motion
+
+            //- Calculate the displacements to relax the points
+            void relaxPoints(const scalar relaxation);
+
+            //- 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();
+
+        // Check
+
+        // Edit
+
+        // Write
+
+            //- Write Delaunay points to .obj file
+            void writePoints(const fileName& fName, bool internalOnly) const;
+
+            //- Write Delaunay triangles as .obj file
+            void writeTriangles(const fileName& fName, bool internalOnly) const;
+
+            //- Write dual points and faces as .obj file
+            void writeDual
+            (
+                const pointField& points,
+                const faceList& faces,
+                const fileName& fName
+            ) const;
+
+            //- Write polyMesh
+            void writeMesh(bool writeToTimestep = false);
+
+            void write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CV3DI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/CV3DI.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline Foam::label Foam::CV3D::insertPoint
+(
+    const point& p,
+    const label type
+)
+{
+    uint nVert = number_of_vertices();
+
+    Vertex_handle vh = insert(toPoint(p));
+
+    if (nVert == number_of_vertices())
+    {
+        WarningIn("Foam::CV3D::insertPoint")
+            << "Failed to insert point " << p << endl;
+    }
+    else
+    {
+        vh->index() = nVert;
+        vh->type() = type;
+    }
+
+    return vh->index();
+}
+
+
+inline bool Foam::CV3D::insertMirrorPoint
+(
+    const point& nearSurfPt,
+    const point& surfPt
+)
+{
+    point mirrorPoint(2*surfPt - nearSurfPt);
+
+    if (qSurf_.outside(mirrorPoint))
+    {
+        insertPoint(mirrorPoint, Vb::MIRROR_POINT);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+inline void Foam::CV3D::insertPointPair
+(
+    const scalar ppDist,
+    const point& surfPt,
+    const vector& n
+)
+{
+    vector ppDistn = ppDist*n;
+
+    label master = insertPoint
+    (
+        surfPt - ppDistn,
+        number_of_vertices() + 1
+    );
+
+    insertPoint(surfPt + ppDistn, master);
+}
+
+
+inline void Foam::CV3D::insertVb(const Vb& v)
+{
+    const Point& Pt(v.point());
+
+    uint nVert = number_of_vertices();
+
+    Vertex_handle vh = insert(Pt);
+
+    if (nVert == number_of_vertices())
+    {
+        FatalErrorIn("Foam::CV3D::insert(const Vb& v")
+        << "Failed to reinsert Vb at " << topoint(Pt)
+            << endl;
+    }
+    else
+    {
+        vh->index() = v.index();
+        vh->type() = v.type();
+    }
+}
+
+
+void Foam::CV3D::movePoint(const Vertex_handle& vh, const point& p)
+{
+    label vi = vh->index();
+
+    // remove(vh);
+    // insert(toPoint(p))->index() = vi;
+
+    move_point(vh, toPoint(p));
+    vh->index() = vi;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+#ifdef CGAL_INEXACT
+
+inline Foam::CV3D::pointFromPoint Foam::CV3D::topoint
+(
+    const Point& P
+) const
+{
+    return reinterpret_cast<pointFromPoint>(P);
+}
+
+inline Foam::CV3D::PointFrompoint Foam::CV3D::toPoint
+(
+    const point& p
+) const
+{
+    return reinterpret_cast<PointFrompoint>(p);
+}
+
+#else
+
+inline Foam::CV3D::pointFromPoint Foam::CV3D::topoint
+(
+    const Point& P
+) const
+{
+    return point
+    (
+        CGAL::to_double(P.x()),
+        CGAL::to_double(P.y()),
+        CGAL::to_double(P.z())
+    );
+}
+
+inline Foam::CV3D::PointFrompoint Foam::CV3D::toPoint
+(
+    const point& p
+) const
+{
+    return Point(p.x(), p.y(), p.z());
+}
+
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/CV3DIO.C

@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::CV3D::writePoints(const fileName& fName, bool internalOnly) const
+{
+    Info<< nl << "Writing points to " << fName << endl;
+    OFstream str(fName);
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!internalOnly || vit->internalOrBoundaryPoint())
+        {
+            meshTools::writeOBJ(str, topoint(vit->point()));
+        }
+    }
+}
+
+
+void Foam::CV3D::writeDual
+(
+    const pointField& points,
+    const faceList& faces,
+    const fileName& fName
+) const
+{
+    Info<< nl << "Writing dual points and faces to " << fName << endl;
+
+    OFstream str(fName);
+
+    forAll(points, p)
+    {
+        meshTools::writeOBJ(str, points[p]);
+    }
+
+    forAll (faces, f)
+    {
+        str<< 'f';
+
+        const face& fP = faces[f];
+
+        forAll(fP, p)
+        {
+            str<< ' ' << fP[p] + 1;
+        }
+
+        str<< nl;
+    }
+}
+
+
+void Foam::CV3D::writeTriangles(const fileName& fName, bool internalOnly) const
+{
+    Info<< nl << "Writing triangles to " << fName << endl;
+    OFstream str(fName);
+
+    labelList vertexMap(number_of_vertices());
+    label verti = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!internalOnly || !vit->farPoint())
+        {
+            vertexMap[vit->index()] = verti++;
+            meshTools::writeOBJ(str, topoint(vit->point()));
+        }
+    }
+
+    for
+    (
+        Triangulation::Finite_facets_iterator fit = finite_facets_begin();
+        fit != finite_facets_end();
+        ++fit
+    )
+    {
+        const Cell_handle& c(fit->first);
+
+        const int& oppositeVertex(fit->second);
+
+        List<label> facetIndices(3,-1);
+
+        bool writeFacet = true;
+
+        for(label i = 0, k = 0;i < 4; i++)
+        {
+            if(i != oppositeVertex)
+            {
+                if(!internalOnly || !c->vertex(i)->farPoint())
+                {
+                    facetIndices[k] = i;
+                    k++;
+                }
+                else
+                {
+                    writeFacet = false;
+                }
+            }
+        }
+
+        if(writeFacet)
+        {
+            str << "f "
+                << vertexMap[c->vertex(facetIndices[0])->index()] + 1
+                << ' ' << vertexMap[c->vertex(facetIndices[1])->index()] + 1
+                << ' ' << vertexMap[c->vertex(facetIndices[2])->index()] + 1
+                << nl;
+        }
+    }
+}
+
+void Foam::CV3D::writeMesh(bool writeToTimestep)
+{
+    pointField points(0);
+    faceList faces(0);
+    labelList owner(0);
+    labelList neighbour(0);
+    wordList patchNames(0);
+    labelList patchSizes(0);
+    labelList patchStarts(0);
+
+    calcDualMesh
+    (
+        points,
+        faces,
+        owner,
+        neighbour,
+        patchNames,
+        patchSizes,
+        patchStarts
+    );
+
+    writeDual(points, faces, "dualMesh.obj");
+
+    IOobject io
+    (
+        Foam::polyMesh::defaultRegion,
+        runTime_.constant(),
+        runTime_,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    if (writeToTimestep)
+    {
+        Info<< nl << "Writing polyMesh to time directory "
+            << runTime_.timeName() << endl;
+
+        io = IOobject
+        (
+            Foam::polyMesh::defaultRegion,
+            runTime_.path()/runTime_.timeName(),
+            runTime_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        );
+    }
+    else
+    {
+        Info<< nl << "Writing polyMesh to constant." << endl;
+    }
+
+    polyMesh pMesh
+    (
+        io,
+        xferMove(points),
+        xferMove(faces),
+        xferMove(owner),
+        xferMove(neighbour)
+    );
+    // polyMesh pMesh
+    // (
+    //     io,
+    //     points,
+    //     faces,
+    //     owner,
+    //     neighbour
+    // );
+
+    List<polyPatch*> patches(patchStarts.size());
+
+    forAll (patches, p)
+    {
+        patches[p] = new polyPatch
+        (
+            patchNames[p],
+            patchSizes[p],
+            patchStarts[p],
+            p,
+            pMesh.boundaryMesh()
+        );
+    }
+
+    pMesh.addPatches(patches);
+
+    if (!pMesh.write())
+    {
+        FatalErrorIn("CV3D::writeMesh()")
+            << "Failed writing polyMesh."
+            << exit(FatalError);
+    }
+
+
+}
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
+  \\    /    O peration      |
+   \\  /     A nd            | Copyright (C) 1991-2007 OpenCFD Ltd.
+    \\/      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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    CV3DMesher
+
+Description
+    Conformal-Voronoi 3D automatic mesher
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CV3D.H"
+#include "IFstream.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.clear();
+    argList::validArgs.append("surface");
+    argList::validOptions.insert("pointsFile", "<filename>");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    IOdictionary CV3DMesherDict
+    (
+        IOobject
+        (
+            "CV3DMesherDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    // Read the surface to conform to
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    querySurface surf
+    (
+        args.args()[1],
+        readScalar(CV3DMesherDict.lookup("includedAngle"))
+    );
+
+    Info<< nl
+        << "Read surface with " << surf.size() << " triangles from file "
+        << args.args()[1] << nl << endl;
+
+    // Read and triangulation
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    CV3D mesh(runTime, surf);
+
+    if (args.options().found("pointsFile"))
+    {
+        fileName pointFileName(args.options()["pointsFile"]);
+        mesh.insertPoints(pointFileName);
+        mesh.insertSurfacePointPairs();
+        mesh.boundaryConform();
+    }
+    else
+    {
+        mesh.insertGrid();
+        mesh.insertSurfacePointPairs();
+        mesh.boundaryConform();
+    }
+
+    scalar relaxation = mesh.meshingControls().relaxationFactorStart;
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        Info<< nl << "Time = " << runTime.timeName()
+            << nl << "relaxation = " << relaxation
+            << endl;
+
+        mesh.relaxPoints(relaxation);
+
+        mesh.insertSurfacePointPairs();
+
+        relaxation -= (relaxation - mesh.meshingControls().relaxationFactorEnd)
+        /max
+        (
+            (runTime.endTime().value() - runTime.timeOutputValue())
+            /runTime.deltaT().value(),
+            1
+        );
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    mesh.write();
+
+    Info << nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/Make/files

@@ -0,0 +1,15 @@
+#include CGAL_FILES
+
+querySurface.C
+CV3D.C
+controls.C
+tolerances.C
+CV3DIO.C
+CV3DMesher.C
+calcDualMesh.C
+insertFeaturePoints.C
+insertSurfaceNearestPointPairs.C
+insertSurfaceNearPointPairs.C
+insertBoundaryConformPointPairs.C
+
+EXE = $(FOAM_APPBIN)/CV3DMesher

applications/utilities/mesh/generation/CV3DMesher/Make/options

@@ -0,0 +1,21 @@
+//EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+
+include $(GENERAL_RULES)/CGAL
+FFLAGS = -DCGAL_FILES='"${CGAL_PATH}/CGAL/files"'
+
+EXE_INC = \
+    ${EXE_FROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${CGAL_INC} \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -ltriSurface \
+    -ldynamicMesh \
+    -lboost_thread-gcc43-mt-1_37

applications/utilities/mesh/generation/CV3DMesher/backup/hardCodedSimpleCubeForPolyTopoChange.H

@@ -0,0 +1,96 @@
+    IOobject io
+    (
+        Foam::polyMesh::defaultRegion,
+        runTime.constant(),
+        runTime,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    polyMesh pMesh(io);
+
+    IOobject io2
+    (
+        Foam::polyMesh::defaultRegion,
+        runTime.timeName(),
+        runTime,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    polyMesh pMesh2(io2, pointField(0), faceList(0), cellList(0));
+
+    polyTopoChange meshCreator(pMesh.boundaryMesh().size());
+
+    meshCreator.addPoint(point(-1,-1,-1), -1, -1, true);
+    meshCreator.addPoint(point(1,-1,-1), -1, -1, true);
+    meshCreator.addPoint(point(1,1,-1), -1, -1, true);
+    meshCreator.addPoint(point(-1,1,-1), -1, -1, true);
+    meshCreator.addPoint(point(-1,-1,1), -1, -1, true);
+    meshCreator.addPoint(point(1,-1,1), -1, -1, true);
+    meshCreator.addPoint(point(1,1,1), -1, -1, true);
+    meshCreator.addPoint(point(-1,1,1), -1, -1, true);
+    
+    meshCreator.addCell(-1, -1, -1, -1, -1);
+
+    labelList faceConstruct(4);
+
+    faceConstruct[0] = 1;
+    faceConstruct[1] = 2;
+    faceConstruct[2] = 6;
+    faceConstruct[3] = 5;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 0;
+    faceConstruct[1] = 4;
+    faceConstruct[2] = 7;
+    faceConstruct[3] = 3;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 2;
+    faceConstruct[1] = 3;
+    faceConstruct[2] = 7;
+    faceConstruct[3] = 6;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 0;
+    faceConstruct[1] = 1;
+    faceConstruct[2] = 5;
+    faceConstruct[3] = 4;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 0;
+    faceConstruct[1] = 3;
+    faceConstruct[2] = 2;
+    faceConstruct[3] = 1;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 4;
+    faceConstruct[1] = 5;
+    faceConstruct[2] = 6;
+    faceConstruct[3] = 7;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    Info<< meshCreator.points() << endl;
+    Info<< meshCreator.faces() << endl;
+    Info<< meshCreator.faceOwner() << endl;
+    Info<< meshCreator.faceNeighbour() << endl;
+
+//     calcDualMesh(meshCreator);
+
+    autoPtr<mapPolyMesh> map = meshCreator.changeMesh(pMesh, false);
+
+    pMesh.updateMesh(map);
+
+    if (!pMesh.write())
+    {
+        FatalErrorIn("CV3D::writeMesh(const Time& runTime)")
+            << "Failed writing polyMesh."
+            << exit(FatalError);
+    }

applications/utilities/mesh/generation/CV3DMesher/backup/indexedVertex_with_displacementSum.H

@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    indexedVertex
+
+Description
+    An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
+    track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedVertex_H
+#define indexedVertex_H
+
+#include <CGAL/Triangulation_3.h>
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class indexedVertex Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Gt, class Vb=CGAL::Triangulation_vertex_base_3<Gt> >
+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_;
+
+        Foam::vector displacementSum_;
+
+
+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::Cell_handle        Cell_handle;
+    typedef typename Vb::Point              Point;
+
+    template<typename TDS2>
+    struct Rebind_TDS
+    {
+        typedef typename Vb::template Rebind_TDS<TDS2>::Other   Vb2;
+        typedef indexedVertex<Gt,Vb2>                           Other;
+    };
+
+    indexedVertex()
+    :
+        Vb(),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        displacementSum_(Foam::vector::zero)
+    {}
+
+    indexedVertex(const Point& p)
+    :
+        Vb(p),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        displacementSum_(Foam::vector::zero)
+    {}
+
+    indexedVertex(const Point& p, Cell_handle f)
+    :
+        Vb(f, p),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        displacementSum_(Foam::vector::zero)
+    {}
+
+    indexedVertex(Cell_handle f)
+    :
+        Vb(f),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        displacementSum_(Foam::vector::zero)
+    {}
+
+
+    int& index()
+    {
+        return index_;
+    }
+
+    int index() const
+    {
+        return index_;
+    }
+
+
+    int& type()
+    {
+        return type_;
+    }
+
+    int type() const
+    {
+        return type_;
+    }
+
+
+    Foam::vector& displacementSum()
+    {
+        return displacementSum_;
+    }
+
+    const Foam::vector& displacementSum() const
+    {
+        return displacementSum_;
+    }
+
+    //- Is point a far-point
+    inline bool farPoint() const
+    {
+        return type_ == FAR_POINT;
+    }
+
+    //- Is point internal, i.e. not on boundary
+    inline bool internalPoint() const
+    {
+        return type_ <= INTERNAL_POINT;
+    }
+
+    //- Is point internal and near the boundary
+    inline bool nearBoundary() const
+    {
+        return type_ == NEAR_BOUNDARY_POINT;
+    }
+
+    //- Set the point to be near the boundary
+    inline void setNearBoundary()
+    {
+        type_ = NEAR_BOUNDARY_POINT;
+    }
+
+    //- Is point a mirror point
+    inline bool mirrorPoint() const
+    {
+        return type_ == MIRROR_POINT;
+    }
+
+    //- Either master or slave of pointPair.
+    inline bool pairPoint() const
+    {
+        return type_ >= 0;
+    }
+
+    //- Master of a pointPair is the lowest numbered one.
+    inline bool ppMaster() const
+    {
+        if (type_ > index_)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    //- Slave of a pointPair is the highest numbered one.
+    inline bool ppSlave() const
+    {
+        if (type_ >= 0 && type_ < index_)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    //- Either original internal point or master of pointPair.
+    inline bool internalOrBoundaryPoint() const
+    {
+        return internalPoint() || ppMaster();
+    }
+
+    //- Is point near the boundary or part of the boundary definition
+    inline bool nearOrOnBoundary() const
+    {
+        return pairPoint() || mirrorPoint() || nearBoundary();
+    }
+
+    //- Do the two given vertices consitute a boundary point-pair
+    inline friend bool pointPair
+    (
+        const indexedVertex& v0,
+        const indexedVertex& v1
+    )
+    {
+        return v0.index_ == v1.type_ || v1.index_ == v0.type_;
+    }
+
+    //- Do the three given vertices consitute a boundary triangle
+    inline friend bool 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));
+    }
+
+    //- Do the three given vertices consitute an outside triangle
+    inline friend bool outsideTriangle
+    (
+        const indexedVertex& v0,
+        const indexedVertex& v1,
+        const indexedVertex& v2
+    )
+    {
+        return (v0.farPoint() || v0.ppSlave())
+            || (v1.farPoint() || v1.ppSlave())
+            || (v2.farPoint() || v2.ppSlave());
+    }
+
+
+    // inline void operator=(const Triangulation::Finite_vertices_iterator vit)
+    // {
+    //     Vb::operator=indexedVertex(vit->point());
+
+    //     this->index_ = vit->index();
+
+    //     this->type_ = vit->type();
+
+    //     this->displacementSum_ = vit->displacementSum();
+    // }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C

@@ -0,0 +1,642 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV3D::calcDualMesh
+(
+    pointField& points,
+    faceList& faces,
+    labelList& owner,
+    labelList& neighbour,
+    wordList& patchNames,
+    labelList& patchSizes,
+    labelList& patchStarts
+)
+{
+    // ~~~~~~~~~~~ removing short edges by indexing dual vertices ~~~~~~~~~~~~~~
+
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        cit->cellIndex() = -1;
+    }
+
+    points.setSize(number_of_cells());
+
+    label dualVerti = 0;
+
+    // Scanning by number of short (dual) edges (nSE) attached to the
+    // circumcentre of each Delaunay tet.  A Delaunay tet may only have four
+    // dual edges emanating from its circumcentre, assigning positions and
+    // indices to those with 4 short edges attached first, then >= 3, then >= 2
+    // etc.
+    for (label nSE = 4; nSE >= 0; nSE--)
+    {
+        Info<< nl << "Scanning for dual vertices with >= "
+            << nSE
+            << " short edges attached." << endl;
+
+        for
+        (
+            Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+            cit != finite_cells_end();
+            ++cit
+        )
+        {
+            // If the Delaunay tet has an index already then it has either
+            // evaluated itself and taken action or has had its index dictated
+            // by a neighbouring tet with more short edges attached.
+
+            if (cit->cellIndex() == -1)
+            {
+                point dualVertex = topoint(dual(cit));
+
+                label shortEdges = 0;
+
+                List<bool> edgeIsShort(4, false);
+
+                List<bool> neighbourAlreadyIndexed(4, false);
+
+                // Loop over the four facets of the Delaunay tet
+                for (label f = 0; f < 4; f++)
+                {
+                    // Check that at least one of the vertices of the facet is
+                    // an internal or boundary point
+                    if
+                    (
+                        cit->vertex(vertex_triple_index(f, 0))->
+                        internalOrBoundaryPoint()
+                     || cit->vertex(vertex_triple_index(f, 1))->
+                        internalOrBoundaryPoint()
+                     || cit->vertex(vertex_triple_index(f, 2))->
+                        internalOrBoundaryPoint()
+                    )
+                    {
+                        point neighDualVertex;
+
+                        label cNI = cit->neighbor(f)->cellIndex();
+
+                        if (cNI == -1)
+                        {
+                            neighDualVertex = topoint(dual(cit->neighbor(f)));
+                        }
+                        else
+                        {
+                            neighDualVertex = points[cNI];
+                        }
+
+                        if
+                        (
+                            magSqr(dualVertex - neighDualVertex)
+                            < tols_.minEdgeLen2
+                        )
+                        {
+                            edgeIsShort[f] = true;
+
+                            if (cNI > -1)
+                            {
+                                neighbourAlreadyIndexed[f] = true;
+                            }
+
+                            shortEdges++;
+                        }
+                    }
+                }
+
+                if (nSE == 0 && shortEdges == 0)
+                {
+                    // Final iteration and no short edges are found, index
+                    // remaining dual vertices.
+
+                    if
+                    (
+                        cit->vertex(0)->internalOrBoundaryPoint()
+                     || cit->vertex(1)->internalOrBoundaryPoint()
+                     || cit->vertex(2)->internalOrBoundaryPoint()
+                     || cit->vertex(3)->internalOrBoundaryPoint()
+                    )
+                    {
+                        cit->cellIndex() = dualVerti;
+                        points[dualVerti] = dualVertex;
+                        dualVerti++;
+                    }
+                }
+                else if
+                (
+                    shortEdges >= nSE
+                )
+                {
+                    // Info<< neighbourAlreadyIndexed << ' '
+                    //     << edgeIsShort << endl;
+
+                    label numUnindexedNeighbours = 1;
+
+                    for (label f = 0; f < 4; f++)
+                    {
+                        if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+                        {
+                            dualVertex += topoint(dual(cit->neighbor(f)));
+
+                            numUnindexedNeighbours++;
+                        }
+                    }
+
+                    dualVertex /= numUnindexedNeighbours;
+
+                    label nearestExistingIndex = -1;
+
+                    point nearestIndexedNeighbourPos = vector::zero;
+
+                    scalar minDistSqrToNearestIndexedNeighbour = VGREAT;
+
+                    for (label f = 0; f < 4; f++)
+                    {
+                        if (edgeIsShort[f] && neighbourAlreadyIndexed[f])
+                        {
+                            label cNI = cit->neighbor(f)->cellIndex();
+
+                            point indexedNeighbourPos = points[cNI];
+
+                            if
+                            (
+                                magSqr(indexedNeighbourPos - dualVertex)
+                                < minDistSqrToNearestIndexedNeighbour
+                            )
+                            {
+                                nearestExistingIndex = cNI;
+
+                                nearestIndexedNeighbourPos =
+                                indexedNeighbourPos;
+
+                                minDistSqrToNearestIndexedNeighbour =
+                                magSqr(indexedNeighbourPos - dualVertex);
+                            }
+                        }
+                    }
+
+                    if
+                    (
+                        nearestExistingIndex > -1
+                     && minDistSqrToNearestIndexedNeighbour < tols_.minEdgeLen2
+                    )
+                    {
+                        points[nearestExistingIndex] =
+                        0.5*(dualVertex + nearestIndexedNeighbourPos);
+
+                        for (label f = 0; f < 4; f++)
+                        {
+                            if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+                            {
+                                cit->neighbor(f)->cellIndex() =
+                                nearestExistingIndex;
+                            }
+                        }
+
+                        cit->cellIndex() = nearestExistingIndex;
+                    }
+                    else
+                    {
+                        for (label f = 0; f < 4; f++)
+                        {
+                            if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+                            {
+                                cit->neighbor(f)->cellIndex() = dualVerti;
+                            }
+                        }
+
+                        cit->cellIndex() = dualVerti;
+
+                        points[dualVerti] = dualVertex;
+
+                        dualVerti++;
+                    }
+                }
+            }
+        }
+    }
+
+    points.setSize(dualVerti);
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~ dual cell indexing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // assigns an index to the Delaunay vertices which will be the dual cell
+    // index used for owner neighbour assignment.
+
+    // The indices of the points are reset which destroys the point-pair
+    // matching, so the type of each vertex are reset to avoid any ambiguity.
+
+    label dualCelli = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalOrBoundaryPoint())
+        {
+            vit->type() = Vb::INTERNAL_POINT;
+            vit->index() = dualCelli;
+            dualCelli++;
+        }
+        else
+        {
+            vit->type() = Vb::FAR_POINT;
+            vit->index() = -1;
+        }
+    }
+
+    // ~~~~~~~~~~~~ dual face and owner neighbour construction ~~~~~~~~~~~~~~~~~
+
+    label nPatches = qSurf_.patches().size() + 1;
+
+    label defaultPatchIndex = qSurf_.patches().size();
+
+    patchNames.setSize(nPatches);
+
+    const geometricSurfacePatchList& surfacePatches = qSurf_.patches();
+
+    forAll(surfacePatches, sP)
+    {
+        patchNames[sP] = surfacePatches[sP].name();
+    }
+
+    patchNames[defaultPatchIndex] = "CV3D_default_patch";
+
+    patchSizes.setSize(nPatches);
+
+    patchStarts.setSize(nPatches);
+
+    List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
+
+    List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
+
+    faces.setSize(number_of_edges());
+
+    owner.setSize(number_of_edges());
+
+    neighbour.setSize(number_of_edges());
+
+    label dualFacei = 0;
+
+    for
+    (
+        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            Cell_circulator ccStart = incident_cells(*eit);
+            Cell_circulator cc1 = ccStart;
+            Cell_circulator cc2 = cc1;
+
+            // Advance the second circulator so that it always stays on the next
+            // cell around the edge;
+            cc2++;
+
+            DynamicList<label> verticesOnFace;
+
+            do
+            {
+                label cc1I = cc1->cellIndex();
+
+                label cc2I = cc2->cellIndex();
+
+
+                if (cc1I < 0 || cc2I < 0)
+                {
+                    FatalErrorIn("Foam::CV3D::calcDualMesh")
+                    << "Dual face uses circumcenter defined by a "
+                        << "Delaunay tetrahedron with no internal "
+                        << "or boundary points."
+                        << exit(FatalError);
+                }
+
+                if (cc1I != cc2I)
+                {
+                    verticesOnFace.append(cc1I);
+                }
+
+                cc1++;
+
+                cc2++;
+            } while (cc1 != ccStart);
+
+            verticesOnFace.shrink();
+
+            if (verticesOnFace.size() >= 3)
+            {
+                face newDualFace(verticesOnFace);
+
+                label dcA = vA->index();
+
+                if (!vA->internalOrBoundaryPoint())
+                {
+                    dcA = -1;
+                }
+
+                label dcB = vB->index();
+
+                if (!vB->internalOrBoundaryPoint())
+                {
+                    dcB = -1;
+                }
+
+                label dcOwn = -1;
+                label dcNei = -1;
+
+                if (dcA == -1 && dcB == -1)
+                {
+                    FatalErrorIn("calcDualMesh")
+                    << "Attempting to create a face joining "
+                        << "two external dual cells "
+                        << exit(FatalError);
+                }
+                else if (dcA == -1 || dcB == -1)
+                {
+                    // boundary face, find which is the owner
+
+                    if (dcA == -1)
+                    {
+                        dcOwn = dcB;
+
+                        // reverse face order to correctly orientate normal
+                        reverse(newDualFace);
+                    }
+                    else
+                    {
+                        dcOwn = dcA;
+                    }
+
+                    // Find which patch this face is on by finding the
+                    // intersection with the surface of the Delaunay edge
+                    // generating the face and identify the region of the
+                    // intersection.
+
+                    point ptA = topoint(vA->point());
+
+                    point ptB = topoint(vB->point());
+
+                    pointIndexHit pHit = qSurf_.tree().findLineAny(ptA, ptB);
+
+                    label patchIndex = qSurf_[pHit.index()].region();
+
+                    if (patchIndex == -1)
+                    {
+                        patchIndex = defaultPatchIndex;
+
+                        WarningIn("Foam::CV3D::calcDualMesh.C")
+                            << "Dual face found that is not on a surface "
+                            << "patch. Adding to CV3D_default_patch."
+                            << endl;
+                    }
+
+                    patchFaces[patchIndex].append(newDualFace);
+                    patchOwners[patchIndex].append(dcOwn);
+                }
+                else
+                {
+                    // internal face, find the lower cell to be the owner
+
+                    if (dcB > dcA)
+                    {
+                        dcOwn = dcA;
+                        dcNei = dcB;
+                    }
+                    else
+                    {
+                        dcOwn = dcB;
+                        dcNei = dcA;
+
+                        // reverse face order to correctly orientate normal
+                        reverse(newDualFace);
+                    }
+
+                    faces[dualFacei] = newDualFace;
+
+                    owner[dualFacei] = dcOwn;
+
+                    neighbour[dualFacei] = dcNei;
+
+                    dualFacei++;
+                }
+            }
+            // else
+            // {
+            //     Info<< verticesOnFace.size()
+            //         << " size face not created." << endl;
+            // }
+        }
+    }
+
+    label nInternalFaces = dualFacei;
+
+    faces.setSize(nInternalFaces);
+
+    owner.setSize(nInternalFaces);
+
+    neighbour.setSize(nInternalFaces);
+
+    // ~~~~~~~~ sort owner, reordinging neighbour and faces to match ~~~~~~~~~~~
+    // two stage sort for upper triangular order:  sort by owner first, then for
+    // each block of owners sort by neighbour
+
+    labelList sortingIndices;
+
+    // Stage 1
+
+    {
+        SortableList<label> sortedOwner(owner);
+
+        sortingIndices = sortedOwner.indices();
+    }
+
+    {
+        labelList copyOwner(owner.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyOwner[sI] = owner[sortingIndices[sI]];
+        }
+
+        owner = copyOwner;
+    }
+
+    {
+        labelList copyNeighbour(neighbour.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyNeighbour[sI] = neighbour[sortingIndices[sI]];
+        }
+
+        neighbour = copyNeighbour;
+    }
+
+    {
+        faceList copyFaces(faces.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyFaces[sI] = faces[sortingIndices[sI]];
+        }
+
+        faces = copyFaces;
+    }
+
+    // Stage 2
+
+    sortingIndices = -1;
+
+    DynamicList<label> ownerCellJumps;
+
+    // Force first owner entry to be a jump
+    ownerCellJumps.append(0);
+
+    for (label o = 1; o < owner.size(); o++)
+    {
+        if (owner[o] > owner[o-1])
+        {
+            ownerCellJumps.append(o);
+        }
+    }
+
+    ownerCellJumps.shrink();
+
+    forAll(ownerCellJumps, oCJ)
+    {
+        label start = ownerCellJumps[oCJ];
+
+        label length;
+
+        if (oCJ == ownerCellJumps.size() - 1)
+        {
+            length = owner.size() - start;
+        }
+        else
+        {
+            length = ownerCellJumps[oCJ + 1] - start;
+        }
+
+        SubList<label> neighbourBlock(neighbour, length, start);
+
+        SortableList<label> sortedNeighbourBlock(neighbourBlock);
+
+        forAll(sortedNeighbourBlock, sNB)
+        {
+            sortingIndices[start + sNB] =
+                sortedNeighbourBlock.indices()[sNB] + start;
+        }
+    }
+
+    // Perform sort
+
+    {
+        labelList copyOwner(owner.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyOwner[sI] = owner[sortingIndices[sI]];
+        }
+
+        owner = copyOwner;
+    }
+
+    {
+        labelList copyNeighbour(neighbour.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyNeighbour[sI] = neighbour[sortingIndices[sI]];
+        }
+
+        neighbour = copyNeighbour;
+    }
+
+    {
+        faceList copyFaces(faces.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyFaces[sI] = faces[sortingIndices[sI]];
+        }
+
+        faces = copyFaces;
+    }
+
+    // ~~~~~~~~ add patch information ~~~~~~~~~~~
+
+    label nBoundaryFaces = 0;
+
+    forAll(patchFaces, p)
+    {
+        patchFaces[p].shrink();
+
+        patchOwners[p].shrink();
+
+        patchSizes[p] = patchFaces[p].size();
+
+        patchStarts[p] = nInternalFaces + nBoundaryFaces;
+
+        nBoundaryFaces += patchSizes[p];
+    }
+
+    faces.setSize(nInternalFaces + nBoundaryFaces);
+
+    owner.setSize(nInternalFaces + nBoundaryFaces);
+
+    forAll(patchFaces, p)
+    {
+        forAll(patchFaces[p], f)
+        {
+            faces[dualFacei] = patchFaces[p][f];
+
+            owner[dualFacei] = patchOwners[p][f];
+
+            dualFacei++;
+        }
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/controls.C

@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CV3D::controls::controls(const IOdictionary& controlDict)
+:
+    relaxationFactorStart
+    (
+        readScalar(controlDict.lookup("relaxationFactorStart"))
+    ),
+    relaxationFactorEnd(readScalar(controlDict.lookup("relaxationFactorEnd"))),
+    minCellSize(readScalar(controlDict.lookup("minCellSize"))),
+    minCellSize2(Foam::sqr(minCellSize)),
+    includedAngle(readScalar(controlDict.lookup("includedAngle"))),
+    maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))),
+    squares(controlDict.lookup("squares")),
+    nearWallAlignedDist
+    (
+        readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize
+    ),
+    nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)),
+    relaxOrientation(controlDict.lookup("relaxOrientation")),
+    insertSurfaceNearestPointPairs
+    (
+        controlDict.lookup("insertSurfaceNearestPointPairs")
+    ),
+    mirrorPoints(controlDict.lookup("mirrorPoints")),
+    insertSurfaceNearPointPairs
+    (
+        controlDict.lookup("insertSurfaceNearPointPairs")
+    ),
+    writeInitialTriangulation(controlDict.lookup("writeInitialTriangulation")),
+    writeFeatureTriangulation(controlDict.lookup("writeFeatureTriangulation")),
+    writeNearestTriangulation(controlDict.lookup("writeNearestTriangulation")),
+    writeInsertedPointPairs(controlDict.lookup("writeInsertedPointPairs")),
+    writeFinalTriangulation(controlDict.lookup("writeFinalTriangulation")),
+    randomiseInitialGrid(controlDict.lookup("randomiseInitialGrid")),
+    randomPerturbation(readScalar(controlDict.lookup("randomPerturbation")))
+{}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/indexedCell.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    indexedCell
+
+Description
+    An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep
+    track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedCell_H
+#define indexedCell_H
+
+#include <CGAL/Triangulation_3.h>
+#include "indexedVertex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class indexedCell Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Gt, class Cb=CGAL::Triangulation_cell_base_3<Gt> >
+class indexedCell
+:
+    public Cb
+{
+    // Private data
+
+        //- The index for this triangle Cell
+        //  -1: triangle and changed and associated data needs updating
+        // >=0: index of cells, set by external update algorithm
+        int index_;
+
+
+public:
+
+    enum cellTypes
+    {
+        UNCHANGED           = 0,
+        CHANGED             = -1,
+        SAVE_CHANGED        = -2
+    };
+
+    typedef typename Cb::Vertex_handle                   Vertex_handle;
+    typedef typename Cb::Cell_handle                     Cell_handle;
+
+    template < typename TDS2 >
+    struct Rebind_TDS
+    {
+        typedef typename Cb::template Rebind_TDS<TDS2>::Other Cb2;
+        typedef indexedCell<Gt, Cb2>  Other;
+    };
+
+
+    indexedCell()
+    :
+        Cb(),
+        index_(CHANGED)
+    {}
+
+    indexedCell
+    (
+        Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
+    )
+    :
+        Cb(v0, v1, v2, v3),
+        index_(CHANGED)
+    {}
+
+    indexedCell
+    (
+        Vertex_handle v0,
+        Vertex_handle v1,
+        Vertex_handle v2,
+        Vertex_handle v3,
+        Cell_handle n0,
+        Cell_handle n1,
+        Cell_handle n2,
+        Cell_handle n3
+    )
+    :
+        Cb(v0, v1, v2, v3, n0, n1, n2, n3),
+        index_(CHANGED)
+    {}
+
+    void set_vertex(int i, Vertex_handle v)
+    {
+        index_ = CHANGED;
+        Cb::set_vertex(i, v);
+    }
+
+    void set_vertices()
+    {
+        index_ = CHANGED;
+        Cb::set_vertices();
+    }
+
+    void set_vertices
+    (
+        Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
+    )
+    {
+        index_ = CHANGED;
+        Cb::set_vertices(v0, v1, v2, v3);
+    }
+
+    int& cellIndex()
+    {
+        return index_;
+    }
+
+    int cellIndex() const
+    {
+        return index_;
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/indexedVertex.H

@@ -0,0 +1,315 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    indexedVertex
+
+Description
+    An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
+    track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedVertex_H
+#define indexedVertex_H
+
+#include <CGAL/Triangulation_3.h>
+#include "List.H"
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class indexedVertex Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Gt, class Vb=CGAL::Triangulation_vertex_base_3<Gt> >
+class indexedVertex
+:
+    public Vb
+{
+    // Private data
+
+        //- The index for this triangle vertex
+        int index_;
+
+        //- type 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_;
+
+        Foam::List<Foam::vector> alignmentDirections_;
+
+        Foam::scalar distanceToClosestSurface_;
+
+        Foam::label indexOfClosestPatch_;
+
+
+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::Cell_handle        Cell_handle;
+    typedef typename Vb::Point              Point;
+
+    template<typename TDS2>
+    struct Rebind_TDS
+    {
+        typedef typename Vb::template Rebind_TDS<TDS2>::Other   Vb2;
+        typedef indexedVertex<Gt,Vb2>                           Other;
+    };
+
+    indexedVertex()
+    :
+        Vb(),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        alignmentDirections_(0),
+        distanceToClosestSurface_(-1),
+        indexOfClosestPatch_(-1)
+    {}
+
+    indexedVertex(const Point& p)
+    :
+        Vb(p),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        alignmentDirections_(0),
+        distanceToClosestSurface_(-1),
+        indexOfClosestPatch_(-1)
+    {}
+
+    indexedVertex(const Point& p, Cell_handle f)
+    :
+        Vb(f, p),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        alignmentDirections_(0),
+        distanceToClosestSurface_(-1),
+        indexOfClosestPatch_(-1)
+    {}
+
+    indexedVertex(Cell_handle f)
+    :
+        Vb(f),
+        index_(INTERNAL_POINT),
+        type_(INTERNAL_POINT),
+        alignmentDirections_(0),
+        distanceToClosestSurface_(-1),
+        indexOfClosestPatch_(-1)
+    {}
+
+
+    int& index()
+    {
+        return index_;
+    }
+
+    int index() const
+    {
+        return index_;
+    }
+
+    int& type()
+    {
+        return type_;
+    }
+
+    int type() const
+    {
+        return type_;
+    }
+
+    inline Foam::List<Foam::vector>& alignmentDirections()
+    {
+        return alignmentDirections_;
+    }
+
+    inline const Foam::List<Foam::vector>& alignmentDirections() const
+    {
+        return alignmentDirections_;
+    }
+
+    inline Foam::scalar& distanceToClosestSurface()
+    {
+        return distanceToClosestSurface_;
+    }
+
+    inline Foam::label& indexOfClosestPatch()
+    {
+        return indexOfClosestPatch_;
+    }
+
+    inline bool uninitialised() const
+    {
+        return type_ == INTERNAL_POINT && index_ == INTERNAL_POINT;
+    }
+
+    //- Is point a far-point
+    inline bool farPoint() const
+    {
+        return type_ == FAR_POINT;
+    }
+
+    //- Is point internal, i.e. not on boundary
+    inline bool internalPoint() const
+    {
+        return type_ <= INTERNAL_POINT;
+    }
+
+    //- Is point internal and near the boundary
+    inline bool nearBoundary() const
+    {
+        return type_ == NEAR_BOUNDARY_POINT;
+    }
+
+    //- Set the point to be near the boundary
+    inline void setNearBoundary()
+    {
+        type_ = NEAR_BOUNDARY_POINT;
+    }
+
+    //- Is point a mirror point
+    inline bool mirrorPoint() const
+    {
+        return type_ == MIRROR_POINT;
+    }
+
+    //- Either master or slave of pointPair.
+    inline bool pairPoint() const
+    {
+        return type_ >= 0;
+    }
+
+    //- Master of a pointPair is the lowest numbered one.
+    inline bool ppMaster() const
+    {
+        if (type_ > index_)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    //- Slave of a pointPair is the highest numbered one.
+    inline bool ppSlave() const
+    {
+        if (type_ >= 0 && type_ < index_)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    //- Either original internal point or master of pointPair.
+    inline bool internalOrBoundaryPoint() const
+    {
+        return internalPoint() || ppMaster();
+    }
+
+    //- Is point near the boundary or part of the boundary definition
+    inline bool nearOrOnBoundary() const
+    {
+        return pairPoint() || mirrorPoint() || nearBoundary();
+    }
+
+    //- Do the two given vertices consitute a boundary point-pair
+    inline friend bool pointPair
+    (
+        const indexedVertex& v0,
+        const indexedVertex& v1
+    )
+    {
+        return v0.index_ == v1.type_ || v1.index_ == v0.type_;
+    }
+
+    //- Do the three given vertices consitute a boundary triangle
+    inline friend bool 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));
+    }
+
+    //- Do the three given vertices consitute an outside triangle
+    inline friend bool outsideTriangle
+    (
+        const indexedVertex& v0,
+        const indexedVertex& v1,
+        const indexedVertex& v2
+    )
+    {
+        return (v0.farPoint() || v0.ppSlave())
+            || (v1.farPoint() || v1.ppSlave())
+            || (v2.farPoint() || v2.ppSlave());
+    }
+
+
+    // inline void operator=(const Triangulation::Finite_vertices_iterator vit)
+    // {
+    //     Vb::operator=indexedVertex(vit->point());
+
+    //     this->index_ = vit->index();
+
+    //     this->type_ = vit->type();
+
+    //     this->displacementSum_ = vit->displacementSum();
+    // }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/insertBoundaryConformPointPairs.C

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV3D::markNearBoundaryPoints()
+{
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        if (vit->internalPoint())
+        {
+            point vert(topoint(vit->point()));
+
+            pointIndexHit pHit =
+                qSurf_.tree().findNearest(vert, 4*controls_.minCellSize2);
+
+            if (pHit.hit())
+            {
+                vit->setNearBoundary();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C

@@ -0,0 +1,615 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV3D::insertFeaturePoints()
+{
+    Info<< nl << "Inserting feature points" << endl;
+
+    const edgeList& edges = qSurf_.edges();
+    const pointField& localPts = qSurf_.localPoints();
+
+    const labelList& featPoints = qSurf_.features().featurePoints();
+    labelListList featPointFeatEdges = qSurf_.featurePointFeatureEdges();
+
+    scalar planeErrorAngle = 0.1*(180.0 - controls_.includedAngle);
+
+    scalar planeErrorAngleCos =
+    cos(mathematicalConstant::pi*planeErrorAngle/180.0);
+
+    forAll(featPoints, i)
+    {
+        label ptI = featPoints[i];
+        const point& featPt = localPts[ptI];
+
+        Info<< nl <<"Feature at " << featPt << endl;
+
+        const labelList& featEdges = featPointFeatEdges[i];
+
+        DynamicList<label> convexEdges(0);
+        DynamicList<label> concaveEdges(0);
+
+        List<vector> planeNormals(2*featEdges.size());
+
+        // Classify the edges around the feature point.
+        forAll(featEdges, fE)
+        {
+            label edgeI = featEdges[fE];
+
+            // Pick up the two faces adjacent to the feature edge
+            const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
+
+            label faceA = eFaces[0];
+            vector nA = qSurf_.faceNormals()[faceA];
+
+            label faceB = eFaces[1];
+            vector nB = qSurf_.faceNormals()[faceB];
+
+            point faceAVert =
+                localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
+
+            // Determine convex or concave angle
+            if (((faceAVert - featPt) & nB) < 0)
+            {
+                // Convex feature edge
+                convexEdges.append(edgeI);
+            }
+            else
+            {
+                // Concave feature edge
+                concaveEdges.append(edgeI);
+            }
+
+            planeNormals[2*fE] = nA;
+            planeNormals[2*fE + 1] = nB;
+        }
+
+        convexEdges.shrink();
+        concaveEdges.shrink();
+
+        // Identify the normals of the faces attached to feature edges to
+        // identify the unique planes to be reconstructed.  The triangulated
+        // surface will usually not mean that two feature edges that should
+        // bound a plane are attached to the same face.
+
+        List<vector> uniquePlaneNormals(featEdges.size());
+
+        List<bool> planeMatchedStatus(2*featEdges.size(), bool(false));
+
+        label uniquePlaneNormalI = 0;
+
+        // Examine the plane normals to identify unique planes.
+        forAll(planeNormals, nA)
+        {
+            const vector& normalA = planeNormals[nA];
+
+            scalar maxNormalDotProduct = -SMALL;
+
+            label matchingNormal = -1;
+
+            if (!planeMatchedStatus[nA])
+            {
+                forAll(planeNormals, nB)
+                {
+                    if (nA == nB)
+                    {
+                        continue;
+                    }
+
+                    if (!planeMatchedStatus[nB])
+                    {
+                        const vector& normalB = planeNormals[nB];
+
+                        scalar normalDotProduct = normalA & normalB;
+
+                        if (normalDotProduct > maxNormalDotProduct)
+                        {
+                            maxNormalDotProduct = normalDotProduct;
+
+                            matchingNormal = nB;
+                        }
+                    }
+                }
+            }
+
+            if (matchingNormal >= 0)
+            {
+                if (maxNormalDotProduct < planeErrorAngleCos)
+                {
+                    FatalErrorIn("insertFeaturePoints")
+                        << "Matching planes are not similar enough "
+                        << "at point located at "
+                        << featPt << nl
+                        << exit(FatalError);
+                }
+
+                const vector& normalB = planeNormals[matchingNormal];
+
+                uniquePlaneNormals[uniquePlaneNormalI] =
+                (normalA + normalB)/(mag(normalA + normalB) + VSMALL);
+
+                uniquePlaneNormalI++;
+
+                planeMatchedStatus[nA] = true;
+                planeMatchedStatus[matchingNormal] = true;
+            }
+        }
+
+        if (concaveEdges.size() == 0)
+        {
+            Info<< tab << "Convex feature, "
+                << convexEdges.size() << " edges." << endl;
+
+            vector cornerNormal = sum(uniquePlaneNormals);
+            cornerNormal /= mag(cornerNormal);
+
+            point internalPt =  featPt - tols_.ppDist*cornerNormal;
+            label internalPtIndex =
+            insertPoint(internalPt, number_of_vertices() + 1);
+
+            forAll (uniquePlaneNormals, uPN)
+            {
+                const vector& n = uniquePlaneNormals[uPN];
+
+                plane planeN = plane(featPt, n);
+
+                point externalPt =
+                internalPt + 2.0 * planeN.distance(internalPt) * n;
+
+                insertPoint(externalPt, internalPtIndex);
+            }
+
+        }
+        else if (convexEdges.size() == 0)
+        {
+            Info<< tab << "Concave feature, "
+                << concaveEdges.size() << " edges." << endl;
+
+            vector cornerNormal = sum(uniquePlaneNormals);
+            cornerNormal /= mag(cornerNormal);
+
+            point externalPt = featPt +  tols_.ppDist*cornerNormal;
+
+            label externalPtIndex = number_of_vertices() + concaveEdges.size();
+
+            label internalPtIndex = -1;
+
+            forAll (uniquePlaneNormals, uPN)
+            {
+                const vector& n = uniquePlaneNormals[uPN];
+
+                plane planeN = plane(featPt, n);
+
+                point internalPt = externalPt
+                - 2.0 * planeN.distance(externalPt) * n;
+
+                internalPtIndex = insertPoint(internalPt, externalPtIndex);
+            }
+
+            insertPoint(externalPt,internalPtIndex);
+        }
+        else
+        {
+            Info<< tab << "Mixed feature: convex, concave = "
+                << convexEdges.size() << ", " << concaveEdges.size() << endl;
+
+            if (convexEdges.size() + concaveEdges.size() > 3)
+            {
+                Info<< concaveEdges.size() + convexEdges.size()
+                    << " mixed edge feature."
+                    << " NOT IMPLEMENTED." << endl;
+            }
+            else if (convexEdges.size() > concaveEdges.size())
+            {
+                // Find which planes are joined to the concave edge
+
+                List<label> concaveEdgePlanes(2,label(-1));
+
+                label concaveEdgeI = concaveEdges[0];
+
+                // Pick up the two faces adjacent to the concave feature edge
+                const labelList& eFaces = qSurf_.edgeFaces()[concaveEdgeI];
+
+                label faceA = eFaces[0];
+
+                vector nA = qSurf_.faceNormals()[faceA];
+
+                scalar maxNormalDotProduct = -SMALL;
+
+                forAll(uniquePlaneNormals, uPN)
+                {
+                    scalar normalDotProduct = nA & uniquePlaneNormals[uPN];
+
+                    if (normalDotProduct > maxNormalDotProduct)
+                    {
+                        maxNormalDotProduct = normalDotProduct;
+
+                        concaveEdgePlanes[0] = uPN;
+                    }
+                }
+
+                label faceB = eFaces[1];
+                vector nB = qSurf_.faceNormals()[faceB];
+
+                maxNormalDotProduct = -SMALL;
+
+                forAll(uniquePlaneNormals, uPN)
+                {
+                    scalar normalDotProduct = nB & uniquePlaneNormals[uPN];
+
+                    if (normalDotProduct > maxNormalDotProduct)
+                    {
+                        maxNormalDotProduct = normalDotProduct;
+
+                        concaveEdgePlanes[1] = uPN;
+                    }
+                }
+
+                const vector& concaveEdgePlaneANormal =
+                uniquePlaneNormals[concaveEdgePlanes[0]];
+
+                const vector& concaveEdgePlaneBNormal =
+                uniquePlaneNormals[concaveEdgePlanes[1]];
+
+                label convexEdgesPlaneI;
+
+                if (findIndex(concaveEdgePlanes, 0) == -1)
+                {
+                    convexEdgesPlaneI = 0;
+                }
+                else if (findIndex(concaveEdgePlanes, 1) == -1)
+                {
+                    convexEdgesPlaneI = 1;
+                }
+                else
+                {
+                    convexEdgesPlaneI = 2;
+                }
+
+                const vector& convexEdgesPlaneNormal =
+                    uniquePlaneNormals[convexEdgesPlaneI];
+
+                const edge& concaveEdge = edges[concaveEdgeI];
+
+                // Check direction of edge, if the feature point is at the end()
+                // the reverse direction.
+
+                scalar edgeDirection = 1.0;
+
+                if (ptI == concaveEdge.end())
+                {
+                    edgeDirection *= -1.0;
+                }
+
+                // Intersect planes parallel to the concave edge planes offset
+                // by ppDist and the plane defined by featPt and the edge
+                // vector.
+                plane planeA
+                (
+                    featPt + tols_.ppDist*concaveEdgePlaneANormal,
+                    concaveEdgePlaneANormal
+                );
+
+                plane planeB
+                (
+                    featPt + tols_.ppDist*concaveEdgePlaneBNormal,
+                    concaveEdgePlaneBNormal
+                );
+
+                point concaveEdgeLocalFeatPt = featPt
+                + tols_.ppDist*edgeDirection
+                * concaveEdge.vec(localPts)/concaveEdge.mag(localPts);
+
+                // Finding the nearest point on the intersecting line to the
+                // edge point.  Floating point errors often encountered using
+                // planePlaneIntersect
+
+                plane planeF(concaveEdgeLocalFeatPt, concaveEdge.vec(localPts));
+
+                point concaveEdgeExternalPt =
+                planeF.planePlaneIntersect(planeA,planeB);
+
+                label concaveEdgeExternalPtI = number_of_vertices() + 4;
+
+                // Redefine planes to be on the feature surfaces to project
+                // through
+
+                planeA = plane(featPt, concaveEdgePlaneANormal);
+
+                planeB = plane(featPt, concaveEdgePlaneBNormal);
+
+                point internalPtA = concaveEdgeExternalPt
+                - 2*planeA.distance(concaveEdgeExternalPt)
+                *concaveEdgePlaneANormal;
+
+                label internalPtAI =
+                insertPoint(internalPtA, concaveEdgeExternalPtI);
+
+                point internalPtB = concaveEdgeExternalPt
+                - 2*planeB.distance(concaveEdgeExternalPt)
+                *concaveEdgePlaneBNormal;
+
+                label internalPtBI =
+                insertPoint(internalPtB, concaveEdgeExternalPtI);
+
+                plane planeC(featPt, convexEdgesPlaneNormal);
+
+                point externalPtD = internalPtA +
+                2*planeC.distance(internalPtA) * convexEdgesPlaneNormal;
+
+                insertPoint(externalPtD, internalPtAI);
+
+                point externalPtE = internalPtB +
+                2*planeC.distance(internalPtB) * convexEdgesPlaneNormal;
+
+                insertPoint(externalPtE, internalPtBI);
+
+                insertPoint(concaveEdgeExternalPt, internalPtAI);
+
+                scalar totalAngle = 180/mathematicalConstant::pi *
+                (
+                    mathematicalConstant::pi +
+                    acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
+                );
+
+                if (totalAngle > controls_.maxQuadAngle)
+                {
+                    // Add additional mitering points
+
+                    vector concaveEdgeNormal =
+                    edgeDirection*concaveEdge.vec(localPts)
+                    /concaveEdge.mag(localPts);
+
+                    scalar angleSign = 1.0;
+
+                    if
+                    (
+                        qSurf_.outside
+                        (
+                            featPt - convexEdgesPlaneNormal*tols_.ppDist
+                        )
+                    )
+                    {
+                        angleSign = -1.0;
+                    }
+
+                    scalar phi = angleSign
+                    *acos(concaveEdgeNormal & -convexEdgesPlaneNormal);
+
+                    scalar guard =
+                    (
+                        1 + sin(phi)*tols_.ppDist/mag
+                        (
+                            concaveEdgeLocalFeatPt - concaveEdgeExternalPt
+                        )
+                    )/cos(phi) - 1;
+
+                    point internalPtF = concaveEdgeExternalPt + (2 + guard)
+                    *(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
+
+                    label internalPtFI =
+                    insertPoint(internalPtF, number_of_vertices() + 1);
+
+                    point externalPtG = internalPtF +
+                    2*planeC.distance(internalPtF) * convexEdgesPlaneNormal;
+                    insertPoint(externalPtG, internalPtFI);
+                }
+            }
+            else
+            {
+                // Find which planes are joined to the convex edge
+
+                List<label> convexEdgePlanes(2,label(-1));
+
+                label convexEdgeI = convexEdges[0];
+
+                // Pick up the two faces adjacent to the convex feature edge
+                const labelList& eFaces = qSurf_.edgeFaces()[convexEdgeI];
+
+                label faceA = eFaces[0];
+                vector nA = qSurf_.faceNormals()[faceA];
+
+                scalar maxNormalDotProduct = -SMALL;
+
+                forAll(uniquePlaneNormals, uPN)
+                {
+                    scalar normalDotProduct = nA & uniquePlaneNormals[uPN];
+
+                    if (normalDotProduct > maxNormalDotProduct)
+                    {
+                        maxNormalDotProduct = normalDotProduct;
+
+                        convexEdgePlanes[0] = uPN;
+                    }
+                }
+
+                label faceB = eFaces[1];
+                vector nB = qSurf_.faceNormals()[faceB];
+
+                maxNormalDotProduct = -SMALL;
+
+                forAll(uniquePlaneNormals, uPN)
+                {
+                    scalar normalDotProduct = nB & uniquePlaneNormals[uPN];
+
+                    if (normalDotProduct > maxNormalDotProduct)
+                    {
+                        maxNormalDotProduct = normalDotProduct;
+
+                        convexEdgePlanes[1] = uPN;
+                    }
+                }
+
+                const vector& convexEdgePlaneANormal =
+                    uniquePlaneNormals[convexEdgePlanes[0]];
+                const vector& convexEdgePlaneBNormal =
+                    uniquePlaneNormals[convexEdgePlanes[1]];
+
+                label concaveEdgesPlaneI;
+
+                if (findIndex(convexEdgePlanes, 0) == -1)
+                {
+                    concaveEdgesPlaneI = 0;
+                }
+                else if (findIndex(convexEdgePlanes, 1) == -1)
+                {
+                    concaveEdgesPlaneI = 1;
+                }
+                else
+                {
+                    concaveEdgesPlaneI = 2;
+                }
+
+                const vector& concaveEdgesPlaneNormal =
+                    uniquePlaneNormals[concaveEdgesPlaneI];
+
+                const edge& convexEdge = edges[convexEdgeI];
+
+                // Check direction of edge, if the feature point is at the end()
+                // the reverse direction.
+
+                scalar edgeDirection = 1.0;
+
+                if (ptI == convexEdge.end())
+                {
+                    edgeDirection *= -1.0;
+                }
+
+                // Intersect planes parallel to the convex edge planes offset by
+                // ppDist and the plane defined by featPt and the edge vector.
+                plane planeA
+                (
+                    featPt - tols_.ppDist*convexEdgePlaneANormal,
+                    convexEdgePlaneANormal
+                );
+
+                plane planeB
+                (
+                    featPt - tols_.ppDist*convexEdgePlaneBNormal,
+                    convexEdgePlaneBNormal
+                );
+
+                point convexEdgeLocalFeatPt = featPt
+                  + tols_.ppDist*edgeDirection
+                  * convexEdge.vec(localPts)/convexEdge.mag(localPts);
+
+
+                // Finding the nearest point on the intersecting line to the
+                // edge point.  Floating point errors often encountered using
+                // planePlaneIntersect
+
+                plane planeF(convexEdgeLocalFeatPt, convexEdge.vec(localPts));
+
+                point convexEdgeInternalPt =
+                planeF.planePlaneIntersect(planeA,planeB);
+
+                planeA = plane(featPt, convexEdgePlaneANormal);
+
+                planeB = plane(featPt, convexEdgePlaneBNormal);
+
+                point externalPtA = convexEdgeInternalPt
+                + 2*planeA.distance(convexEdgeInternalPt)
+                * convexEdgePlaneANormal;
+                label externalPtAI = number_of_vertices() + 3;
+
+                point externalPtB = convexEdgeInternalPt
+                + 2*planeB.distance(convexEdgeInternalPt)
+                * convexEdgePlaneBNormal;
+
+                label externalPtBI = number_of_vertices() + 4;
+
+                label convexEdgeInternalPtI = insertPoint
+                (
+                    convexEdgeInternalPt,
+                    externalPtAI
+                );
+
+                plane planeC(featPt, concaveEdgesPlaneNormal);
+
+                point internalPtD = externalPtA -
+                    2*planeC.distance(externalPtA) * concaveEdgesPlaneNormal;
+                insertPoint(internalPtD, externalPtAI);
+
+                point internalPtE = externalPtB -
+                    2*planeC.distance(externalPtB) * concaveEdgesPlaneNormal;
+                insertPoint(internalPtE, externalPtBI);
+
+                insertPoint(externalPtA, convexEdgeInternalPtI);
+
+                insertPoint(externalPtB, convexEdgeInternalPtI);
+            }
+        }
+    }
+
+    featureConstrainingVertices_.setSize(number_of_vertices());
+
+    label featPtI = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        featureConstrainingVertices_[featPtI] = Vb(vit->point());
+
+        featureConstrainingVertices_[featPtI].index() = vit->index();
+
+        featureConstrainingVertices_[featPtI].type() = vit->type();
+
+        featPtI++;
+    }
+
+    if (controls_.writeFeatureTriangulation)
+    {
+        writePoints("feat_allPoints.obj", false);
+        writePoints("feat_points.obj", true);
+        writeTriangles("feat_triangles.obj", true);
+    }
+}
+
+
+void Foam::CV3D::reinsertFeaturePoints()
+{
+    if (featureConstrainingVertices_.size())
+    {
+        forAll(featureConstrainingVertices_, f)
+        {
+            insertVb(featureConstrainingVertices_[f]);
+        }
+    }
+    else
+    {
+        WarningIn("Foam::CV3D::reinsertFeaturePoints")
+            << "No stored feature points to reinsert." << endl;
+    }
+}
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearPointPairs.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV3D::insertSurfaceNearPointPairs()
+{
+    Info<< "insertSurfaceNearPointPairs: " << nl << endl;
+
+    label nNearPoints = 0;
+
+//     This needs a rethink from scratch in 3D.  
+//     Probably iterating around the facets  rather than edges
+//     Remember about the incident_edges, cells etc
+
+    Info<< nNearPoints << " point-pairs inserted" << endl;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C

@@ -0,0 +1,988 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::CV3D::dualCellSurfaceIntersection
+(
+    const Triangulation::Finite_vertices_iterator& vit
+) const
+{
+    std::list<Facet>  facets;
+    incident_facets(vit, std::back_inserter(facets));
+
+    for
+    (
+        std::list<Facet>::iterator fit=facets.begin();
+        fit != facets.end();
+        ++fit
+    )
+    {
+        if
+        (
+            is_infinite(fit->first)
+         || is_infinite(fit->first->neighbor(fit->second))
+        )
+        {
+            return true;
+        }
+
+        point dE0 = topoint(dual(fit->first));
+
+        // If edge end is outside bounding box then edge cuts boundary
+        if (!qSurf_.bb().contains(dE0))
+        {
+            return true;
+        }
+
+        point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
+
+        // If other edge end is outside bounding box then edge cuts boundary
+        if (!qSurf_.bb().contains(dE1))
+        {
+            return true;
+        }
+
+        if (magSqr(dE1 - dE0) > tols_.minEdgeLen2)
+        {
+            pointIndexHit pHit = qSurf_.tree().findLineAny(dE0, dE1);
+
+            if (pHit.hit())
+            {
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+
+void Foam::CV3D::smoothEdgePositions
+(
+    DynamicList<point>& edgePoints,
+    DynamicList<label>& edgeLabels
+) const
+{
+    const pointField& localPts = qSurf_.localPoints();
+    const edgeList& edges = qSurf_.edges();
+
+    // Sort by edge label then by distance from the start of the edge
+
+    SortableList<label> sortedEdgeLabels(edgeLabels);
+
+    labelList sortingIndices = sortedEdgeLabels.indices();
+
+    {
+        labelList copyEdgeLabels(edgeLabels.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgeLabels[sI] = edgeLabels[sortingIndices[sI]];
+        }
+
+        edgeLabels.transfer(copyEdgeLabels);
+    }
+
+    {
+        List<point> copyEdgePoints(edgePoints.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgePoints[sI] = edgePoints[sortingIndices[sI]];
+        }
+
+        edgePoints.transfer(copyEdgePoints);
+    }
+
+    List<scalar> edgeDistances(edgePoints.size());
+
+    forAll(edgeDistances, eD)
+    {
+        const point& edgeStart = localPts[edges[edgeLabels[eD]].start()];
+
+        edgeDistances[eD] = mag(edgeStart - edgePoints[eD]);
+    }
+
+    // Sort by edgeDistances in blocks of edgeLabel
+
+    DynamicList<label> edgeLabelJumps;
+
+    // Force first edgeLabel to be a jump
+    edgeLabelJumps.append(0);
+
+    for (label eL = 1; eL < edgeLabels.size(); eL++)
+    {
+        if (edgeLabels[eL] > edgeLabels[eL-1])
+        {
+            edgeLabelJumps.append(eL);
+        }
+    }
+
+    edgeLabelJumps.shrink();
+
+    forAll(edgeLabelJumps, eLJ)
+    {
+        label start = edgeLabelJumps[eLJ];
+
+        label length;
+
+        if (eLJ == edgeLabelJumps.size() - 1)
+        {
+            length = edgeLabels.size() - start;
+        }
+        else
+        {
+            length = edgeLabelJumps[eLJ + 1] - start;
+        }
+
+        SubList<scalar> edgeDistanceBlock(edgeDistances, length, start);
+
+        SortableList<scalar> sortedEdgeDistanceBlock(edgeDistanceBlock);
+
+        forAll(sortedEdgeDistanceBlock, sEDB)
+        {
+            sortingIndices[start + sEDB] =
+                sortedEdgeDistanceBlock.indices()[sEDB] + start;
+        }
+    }
+
+    {
+        List<point> copyEdgePoints(edgePoints.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgePoints[sI] = edgePoints[sortingIndices[sI]];
+        }
+
+        edgePoints.transfer(copyEdgePoints);
+    }
+
+    {
+        List<scalar> copyEdgeDistances(edgeDistances.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyEdgeDistances[sI] = edgeDistances[sortingIndices[sI]];
+        }
+
+        edgeDistances.transfer(copyEdgeDistances);
+    }
+
+    // Create a List to hold each edge, process them individually, then
+    // recombine the edges into a single list.
+
+    List<List<point> > edgePointIndividualLists(edgeLabelJumps.size());
+
+    List<List<scalar> > edgeDistanceIndividualLists(edgeLabelJumps.size());
+
+    List<label> edgeLabelIndividualList(edgeLabelJumps.size());
+
+    forAll(edgeLabelJumps, eLJ)
+    {
+        label start = edgeLabelJumps[eLJ];
+
+        label edgeI = edgeLabels[start];
+
+        label length;
+
+        if (eLJ == edgeLabelJumps.size() - 1)
+        {
+            length = edgeLabels.size() - start;
+        }
+        else
+        {
+            length = edgeLabelJumps[eLJ + 1] - start;
+        }
+
+        edgePointIndividualLists[eLJ] = SubList<point>
+        (
+            edgePoints,
+            length,
+            start
+        );
+
+        edgeDistanceIndividualLists[eLJ] = SubList<scalar>
+        (
+            edgeDistances,
+            length,
+            start
+        );
+
+        edgeLabelIndividualList[eLJ] = edgeI;
+    }
+
+    edgePoints.clear();
+    edgeDistances.clear();
+    edgeLabels.clear();
+
+    forAll(edgeLabelIndividualList, e)
+    {
+        label edgeI = edgeLabelIndividualList[e];
+
+        smoothEdge
+        (
+            edgePointIndividualLists[e],
+            edgeDistanceIndividualLists[e],
+            edgeI
+        );
+
+        const List<point>& tempEdgePoints = edgePointIndividualLists[e];
+
+        forAll(tempEdgePoints, tEP)
+        {
+            edgePoints.append(tempEdgePoints[tEP]);
+
+            edgeLabels.append(edgeI);
+        }
+    }
+
+    edgePoints.shrink();
+
+    edgeLabels.shrink();
+}
+
+
+
+void Foam::CV3D::smoothEdge
+(
+    List<point>& edgePoints,
+    List<scalar>& edgeDistances,
+    const label edgeI
+) const
+{
+    const pointField& localPts = qSurf_.localPoints();
+    const edgeList& edges = qSurf_.edges();
+
+    // Process the points along each edge (in blocks of edgeLabel) performing 3
+    // functions:
+    // 1: move points away from feature points
+    // 2: aggregate tight groups of points into one point
+    // 3: adjust the spacing of remaining points on a pair by pair basis to
+    //    remove excess points and add points to long uncontrolled spans.
+
+
+    const edge& e(edges[edgeI]);
+
+    const point& eStart(localPts[e.start()]);
+
+    const point& eEnd(localPts[e.end()]);
+
+    scalar edgeLength = mag(eStart - eEnd);
+
+    if (edgeLength < 2*tols_.featurePointGuard)
+    {
+        Info<< "edge " << edgeI
+            << " is too short with respect to the featurePointGuard "
+            << "distance to allow edge control points to be placed."
+            << nl << "Edge length = " << edgeLength
+            << nl <<endl;
+
+        return;
+    }
+
+    // part 1
+    {
+        DynamicList<point> tempEdgePoints;
+
+        bool startGuardPlaced = false;
+
+        bool endGuardPlaced = false;
+
+        forAll (edgePoints, eP)
+        {
+            const point& edgePoint = edgePoints[eP];
+
+            const scalar& edgeDist = edgeDistances[eP];
+
+            if
+            (
+                edgeDist < tols_.featurePointGuard
+             && !startGuardPlaced
+            )
+            {
+                tempEdgePoints.append
+                (
+                    eStart + (edgePoint - eStart)
+                    * tols_.featurePointGuard/edgeDist
+                );
+
+                startGuardPlaced = true;
+            }
+            else if
+            (
+                edgeDist > (edgeLength - tols_.featurePointGuard)
+             && !endGuardPlaced
+            )
+            {
+                tempEdgePoints.append
+                (
+                    eEnd + (edgePoint - eEnd)
+                    * tols_.featurePointGuard/(edgeLength - edgeDist)
+                );
+
+                endGuardPlaced = true;
+            }
+            else if
+            (
+                edgeDist > tols_.featurePointGuard
+             && edgeDist < (edgeLength - tols_.featurePointGuard)
+            )
+            {
+                tempEdgePoints.append(edgePoint);
+            }
+        }
+
+        edgePoints.transfer(tempEdgePoints.shrink());
+    }
+
+    // Recalculate edge distances.
+
+    edgeDistances.setSize(edgePoints.size());
+
+    forAll(edgeDistances, eD)
+    {
+        edgeDistances[eD] = mag(eStart - edgePoints[eD]);
+    }
+
+    // part 2
+    {
+        DynamicList<point> tempEdgePoints;
+
+        label groupSize = 0;
+
+        point newEdgePoint(vector::zero);
+
+        // if (edgePoints.size() == 1)
+        // {
+        //     tempEdgePoints.append(edgePoints[0]);
+        // }
+        // else if
+        // (
+        //     (edgeDistances[1] - edgeDistances[0]) > tols_.edgeGroupSpacing
+        // )
+        // {
+        //     tempEdgePoints.append(edgePoints[0]);
+        // }
+
+        if (edgePoints.size() > 1)
+        {
+            if ((edgeDistances[1] - edgeDistances[0]) < tols_.edgeGroupSpacing)
+            {
+                // ...the first two points on the edge start a group
+
+                newEdgePoint += edgePoints[0];
+
+                groupSize++;
+            }
+            else
+            {
+                tempEdgePoints.append(edgePoints[0]);
+            }
+        }
+        else
+        {
+            // ...add the first point by default
+
+            tempEdgePoints.append(edgePoints[0]);
+        }
+
+        for (label eP = 1; eP < edgePoints.size(); eP++)
+        {
+            const scalar& edgeDist = edgeDistances[eP];
+            const scalar& previousEdgeDist = edgeDistances[eP - 1];
+
+            if ((edgeDist - previousEdgeDist) < tols_.edgeGroupSpacing)
+            {
+                newEdgePoint += edgePoints[eP];
+
+                groupSize++;
+            }
+            else if (groupSize > 0)
+            {
+                // A point group has been formed and has finished
+
+                newEdgePoint /= groupSize;
+
+                tempEdgePoints.append(newEdgePoint);
+
+                newEdgePoint = vector::zero;
+
+                groupSize = 0;
+            }
+            else
+            {
+                tempEdgePoints.append(edgePoints[eP]);
+            }
+        }
+
+        if (groupSize > 0)
+        {
+            // A point group has been formed at the end of the edge and needs to
+            // be finished.
+
+            newEdgePoint /= groupSize;
+
+            tempEdgePoints.append(newEdgePoint);
+        }
+
+        edgePoints.transfer(tempEdgePoints.shrink());
+    }
+
+    // Recalculate edge distances.
+
+    edgeDistances.setSize(edgePoints.size());
+
+    forAll(edgeDistances, eD)
+    {
+        edgeDistances[eD] = mag(eStart - edgePoints[eD]);
+    }
+
+    // part 3
+    {
+        // Special treatment for gaps between closest point to start
+
+        DynamicList<point> tempEdgePoints;
+
+        if (edgeDistances[0] - tols_.featurePointGuard > tols_.maxEdgeSpacing)
+        {
+            scalar gap = edgeDistances[0] - tols_.featurePointGuard;
+
+            label nInsertions = label(gap/tols_.maxEdgeSpacing);
+
+            // Info<< "Gap at start of edge of "  << gap
+            //     << ". Inserting " << nInsertions << " points" << endl;
+
+            scalar spacing = gap / (nInsertions + 1);
+
+            for (label nI = 1; nI <= nInsertions; nI++)
+            {
+                tempEdgePoints.append
+                (
+                    eStart + (eEnd - eStart)
+                    * (nI * spacing + tols_.featurePointGuard) /edgeLength
+                );
+            }
+        }
+
+        // Identify gaps in middle of edges.
+        // Insert first point by default.
+
+        tempEdgePoints.append(edgePoints[0]);
+
+        for (label eP = 1; eP < edgePoints.size(); eP++)
+        {
+            const scalar& edgeDist = edgeDistances[eP];
+            const scalar& previousEdgeDist = edgeDistances[eP - 1];
+
+            if ((edgeDist - previousEdgeDist) > tols_.maxEdgeSpacing)
+            {
+                scalar gap = edgeDist - previousEdgeDist;
+
+                label nInsertions = label(gap/tols_.maxEdgeSpacing);
+
+                // Info<< "Gap in edge of "  << gap
+                //     << ". Inserting " << nInsertions << " points" << endl;
+
+                scalar spacing = gap / (nInsertions + 1);
+
+                for (label nI = 1; nI<= nInsertions; nI++)
+                {
+                    tempEdgePoints.append
+                    (
+                        eStart + (eEnd - eStart)
+                        * (nI * spacing + previousEdgeDist) /edgeLength
+                    );
+                }
+
+
+            }
+
+            tempEdgePoints.append(edgePoints[eP]);
+        }
+
+        // Special treatment for gaps between closest point to end
+
+        if
+        (
+            (edgeLength - edgeDistances[edgeDistances.size() - 1]
+                - tols_.featurePointGuard)
+                > tols_.maxEdgeSpacing
+        )
+        {
+            scalar lastPointDist =  edgeDistances[edgeDistances.size() - 1];
+
+            const point& lastPoint = edgePoints[edgePoints.size() - 1];
+
+            scalar gap = edgeLength - lastPointDist - tols_.featurePointGuard;
+
+            label nInsertions = label(gap/tols_.maxEdgeSpacing);
+
+            // Info<< "Gap at end of edge of "  << gap
+            //     << ". Inserting " << nInsertions << " points" << endl;
+
+            scalar spacing = gap / (nInsertions + 1);
+
+            for (label nI = 1; nI <= nInsertions; nI++)
+            {
+                tempEdgePoints.append
+                (
+                    lastPoint + (eEnd - lastPoint)
+                        * nI * spacing / gap
+                );
+            }
+        }
+
+        edgePoints.transfer(tempEdgePoints.shrink());
+
+        // Remove pairs of points that are too close together.
+
+        label nPointsRemoved = 1;
+
+        while (nPointsRemoved > 0)
+        {
+            nPointsRemoved = 0;
+
+            // Recalculate edge distances.
+
+            edgeDistances.setSize(edgePoints.size());
+
+            forAll(edgeDistances, eD)
+            {
+                edgeDistances[eD] = mag(eStart - edgePoints[eD]);
+            }
+
+            // Insert first point
+            tempEdgePoints.append(edgePoints[0]);
+
+            bool previousPointMustBeKept = false;
+
+            for (label eP = 1; eP < edgePoints.size(); eP++)
+            {
+                const scalar& edgeDist = edgeDistances[eP];
+                const scalar& previousEdgeDist = edgeDistances[eP - 1];
+
+                if ((edgeDist - previousEdgeDist) < tols_.minEdgeSpacing)
+                {
+                    if (!previousPointMustBeKept)
+                    {
+                        tempEdgePoints.remove();
+                    }
+
+                    const point& currentPoint = edgePoints[eP];
+
+                    const point& previousPoint = edgePoints[eP - 1];
+
+                    tempEdgePoints.append(0.5*(previousPoint + currentPoint));
+
+                    nPointsRemoved++;
+
+                    previousPointMustBeKept = true;
+                }
+                else
+                {
+                    tempEdgePoints.append(edgePoints[eP]);
+
+                    previousPointMustBeKept = false;
+                }
+            }
+
+            // Info<< edgeI << tab
+            //<< nPointsRemoved << " points removed." << endl;
+
+            edgePoints.transfer(tempEdgePoints.shrink());
+        }
+    }
+}
+
+
+void Foam::CV3D::insertPointPairs
+(
+    const DynamicList<point>& nearSurfacePoints,
+    const DynamicList<point>& surfacePoints,
+    const DynamicList<label>& surfaceTris,
+    const fileName fName
+)
+{
+    if (controls_.mirrorPoints)
+    {
+        forAll(surfacePoints, ppi)
+        {
+            insertMirrorPoint
+            (
+                nearSurfacePoints[ppi],
+                surfacePoints[ppi]
+            );
+        }
+    }
+    else
+    {
+        forAll(surfacePoints, ppi)
+        {
+            insertPointPair
+            (
+                tols_.ppDist,
+                surfacePoints[ppi],
+                qSurf_.faceNormals()[surfaceTris[ppi]]
+            );
+        }
+    }
+
+    Info<< surfacePoints.size() << " point-pairs inserted" << endl;
+
+    if (controls_.writeInsertedPointPairs)
+    {
+        OFstream str(fName);
+        label vertI = 0;
+
+        forAll(surfacePoints, ppi)
+        {
+            meshTools::writeOBJ(str, surfacePoints[ppi]);
+            vertI++;
+        }
+
+        Info<< "insertPointPairs: Written " << surfacePoints.size()
+            << " inserted point-pair locations to file "
+            << str.name() << endl;
+    }
+}
+
+
+void Foam::CV3D::insertEdgePointGroups
+(
+    const DynamicList<point>& edgePoints,
+    const DynamicList<label>& edgeLabels,
+    const fileName fName
+)
+{
+    const pointField& localPts = qSurf_.localPoints();
+
+    forAll(edgePoints, eP)
+    {
+        const point& edgePt = edgePoints[eP];
+
+        const label edgeI = edgeLabels[eP];
+
+        // Pick up the two faces adjacent to the feature edge
+        const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
+
+        label faceA = eFaces[0];
+        vector nA = qSurf_.faceNormals()[faceA];
+
+        label faceB = eFaces[1];
+        vector nB = qSurf_.faceNormals()[faceB];
+
+        // Intersect planes parallel to faceA and faceB offset by ppDist
+        // and the plane defined by edgePt and the edge vector.
+        plane planeA(edgePt - tols_.ppDist*nA, nA);
+        plane planeB(edgePt - tols_.ppDist*nB, nB);
+
+        // Finding the nearest point on the intersecting line to the edge point.
+        // Floating point errors often encountered using planePlaneIntersect
+        // plane planeF(edgePt, (nA^nB));
+        // point refPt = planeF.planePlaneIntersect(planeA,planeB);
+
+        plane::ray planeIntersect(planeA.planeIntersect(planeB));
+
+        pointHit refPtHit = linePointRef
+        (
+            planeIntersect.refPoint() + 2*tols_.span*planeIntersect.dir(),
+            planeIntersect.refPoint() - 2*tols_.span*planeIntersect.dir()
+        ).nearestDist(edgePt);
+
+        point refPt = refPtHit.hitPoint();
+
+        point faceAVert =
+        localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
+
+        // Determine convex or concave angle
+        if (((faceAVert - edgePt) & nB) < 0)
+        {
+            // Convex. So refPt will be inside domain and hence a master point
+
+            // Insert the master point refering the the first slave
+            label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
+
+            // Insert the slave points by reflecting refPt in both faces.
+            // with each slave refering to the master
+
+            point reflectedA = refPt + 2*((edgePt - refPt) & nA)*nA;
+            insertPoint(reflectedA, masterPtIndex);
+
+            point reflectedB = refPt + 2*((edgePt - refPt) & nB)*nB;
+            insertPoint(reflectedB, masterPtIndex);
+        }
+        else
+        {
+            // Concave. master and reflected points inside the domain.
+            // Generate reflected master to be outside.
+            point reflMasterPt = refPt + 2*(edgePt - refPt);
+
+            // Reflect refPt in both faces.
+            point reflectedA =
+            reflMasterPt + 2*((edgePt - reflMasterPt) & nA)*nA;
+
+            point reflectedB =
+            reflMasterPt + 2*((edgePt - reflMasterPt) & nB)*nB;
+
+            scalar totalAngle =
+            180*(mathematicalConstant::pi + acos(mag(nA & nB)))
+            /mathematicalConstant::pi;
+
+            // Number of quadrants the angle should be split into
+            int nQuads = int(totalAngle/controls_.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 = number_of_vertices() + 2 + nAddPoints;
+
+            // Master A is inside.
+            label reflectedAI = insertPoint(reflectedA, reflectedMaster);
+
+            // Master B is inside.
+            insertPoint(reflectedB, reflectedMaster);
+
+            if (nAddPoints == 1)
+            {
+                // One additinal point is the reflection of the slave point,
+                // i.e. the original reference point
+                insertPoint(refPt, reflectedMaster);
+            }
+            else if (nAddPoints == 2)
+            {
+                point reflectedAa = refPt
+                - ((edgePt - reflMasterPt) & nB)*nB;
+                insertPoint(reflectedAa, reflectedMaster);
+
+                point reflectedBb = refPt
+                - ((edgePt - reflMasterPt) & nA)*nA;
+                insertPoint(reflectedBb, reflectedMaster);
+            }
+
+            // Slave is outside.
+            insertPoint(reflMasterPt, reflectedAI);
+        }
+    }
+
+    Info<< edgePoints.size() << " edge-control locations inserted" << endl;
+
+    if (controls_.writeInsertedPointPairs)
+    {
+        OFstream str(fName);
+        label vertI = 0;
+
+        forAll(edgePoints, eP)
+        {
+            meshTools::writeOBJ(str, edgePoints[eP]);
+            vertI++;
+        }
+
+        Info<< "insertEdgePointGroups: Written " << edgePoints.size()
+            << " inserted edge-control locations to file "
+            << str.name() << endl;
+    }
+}
+
+
+void Foam::CV3D::insertSurfaceNearestPointPairs()
+{
+    Info<< nl << "insertSurfaceNearestPointPairs: " << endl;
+
+    label nSurfacePointsEst = number_of_vertices();
+
+    scalar distanceFactor = 8.0;
+    scalar distanceFactor2 = Foam::sqr(distanceFactor);
+
+    if (qSurf_.features().featureEdges().size())
+    {
+        DynamicList<point> nearSurfacePoints(nSurfacePointsEst);
+
+        for
+        (
+            Triangulation::Finite_vertices_iterator vit =
+            finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            if (vit->internalPoint())
+            {
+                point vert(topoint(vit->point()));
+
+                pointIndexHit pHit = qSurf_.tree().findNearest
+                (
+                    vert,
+                    distanceFactor2*controls_.minCellSize2
+                );
+
+                if (pHit.hit())
+                {
+                    vit->setNearBoundary();
+
+                    if (dualCellSurfaceIntersection(vit))
+                    {
+                        nearSurfacePoints.append(vert);
+                    }
+                }
+            }
+        }
+
+        pointField nearSurfacePointsForEdges(nearSurfacePoints.shrink());
+
+        labelList allEdgeLabels;
+        labelList allEdgeEndPoints;
+        pointField allEdgePoints;
+
+        qSurf_.features().nearestSurfEdge
+        (
+            qSurf_.features().featureEdges(),
+            nearSurfacePointsForEdges,
+            vector::one * distanceFactor * controls_.minCellSize,
+            allEdgeLabels,
+            allEdgeEndPoints,
+            allEdgePoints
+        );
+
+        DynamicList<label> edgeLabels(allEdgeLabels.size());
+        DynamicList<vector> edgePoints(allEdgePoints.size());
+
+        forAll(allEdgePoints, eP)
+        {
+            if (allEdgeLabels[eP] >= 0 && allEdgeEndPoints[eP] < 0)
+            {
+                edgeLabels.append(allEdgeLabels[eP]);
+                edgePoints.append(allEdgePoints[eP]);
+            }
+        }
+
+        edgePoints.shrink();
+        edgeLabels.shrink();
+
+        if (edgePoints.size())
+        {
+            Warning<< "Edge point insertion disabled." << endl;
+
+            // smoothEdgePositions(edgePoints, edgeLabels);
+
+            // insertEdgePointGroups
+            // (
+            //     edgePoints,
+            //     edgeLabels,
+            //     "surfaceNearestEdgePoints.obj"
+            // );
+        }
+    }
+
+    DynamicList<point> allNearSurfacePoints(nSurfacePointsEst);
+    DynamicList<point> allSurfacePoints(nSurfacePointsEst);
+    DynamicList<label> allSurfaceTris(nSurfacePointsEst);
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        if (vit->internalPoint())
+        {
+            point vert(topoint(vit->point()));
+
+            pointIndexHit pHit = qSurf_.tree().findNearest
+            (
+                vert,
+                distanceFactor2*controls_.minCellSize2
+            );
+
+            if (pHit.hit())
+            {
+                vit->setNearBoundary();
+
+                if (dualCellSurfaceIntersection(vit))
+                {
+                    allNearSurfacePoints.append(vert);
+                    allSurfacePoints.append(pHit.hitPoint());
+                    allSurfaceTris.append(pHit.index());
+                }
+            }
+        }
+    }
+
+    pointField surfacePointsForEdges(allSurfacePoints.shrink());
+
+    labelList allEdgeLabels;
+    labelList allEdgeEndPoints;
+    pointField allEdgePoints;
+
+    qSurf_.features().nearestSurfEdge
+    (
+        qSurf_.features().featureEdges(),
+        surfacePointsForEdges,
+        vector::one * 2.0 * tols_.featureEdgeGuard,
+        allEdgeLabels,
+        allEdgeEndPoints,
+        allEdgePoints
+    );
+
+    DynamicList<point> nearSurfacePoints(nSurfacePointsEst);
+    DynamicList<point> surfacePoints(nSurfacePointsEst);
+    DynamicList<label> surfaceTris(nSurfacePointsEst);
+
+    forAll(allEdgePoints, eP)
+    {
+        if (allEdgeLabels[eP] == -1)
+        {
+            nearSurfacePoints.append(allNearSurfacePoints[eP]);
+            surfacePoints.append(allSurfacePoints[eP]);
+            surfaceTris.append(allSurfaceTris[eP]);
+        }
+    }
+
+    Info<< nl << "Number of surface conformation points not placed because "
+        << nl << "they were too close to a feature edge = "
+        << allSurfacePoints.size() - surfacePoints.size()
+        << endl;
+
+    nearSurfacePoints.shrink();
+    surfacePoints.shrink();
+    surfaceTris.shrink();
+
+    insertPointPairs
+    (
+        nearSurfacePoints,
+        surfacePoints,
+        surfaceTris,
+        "surfaceNearestIntersections.obj"
+    );
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/querySurface.C

@@ -0,0 +1,248 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "querySurface.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::querySurface::querySurface
+(
+    const fileName& surfaceFileName,
+    const scalar& includedAngle
+)
+:
+    triSurface(surfaceFileName),
+    rndGen_(12345),
+    bb_(localPoints()),
+    tree_
+    (
+        treeDataTriSurface(*this),
+        bb_.extend(rndGen_, 1e-3), // slightly randomize bb
+        8,                         // maxLevel
+        4, //10,                   // leafsize
+        10.0 //3.0                 // duplicity
+    ),
+    sFeat_(*this, includedAngle)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::querySurface::~querySurface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelListList Foam::querySurface::featurePointFeatureEdges() const
+{
+    const labelList& featPoints = features().featurePoints();
+
+    const labelList& featEdges = features().featureEdges();
+
+    const edgeList& edges = this->edges();
+
+    List<DynamicList<label> > tempFeatPointFeatEdges(featPoints.size());
+
+    forAll(featPoints, pI)
+    {
+        label fP = featPoints[pI];
+
+        forAll(featEdges, eI)
+        {
+            label fE = featEdges[eI];
+
+            const edge& e(edges[fE]);
+
+            if (e.start() == fP || e.end() == fP)
+            {
+                tempFeatPointFeatEdges[pI].append(fE);
+            }
+        }
+    }
+
+    labelListList featPointFeatEdges(tempFeatPointFeatEdges.size());
+
+    forAll(featPointFeatEdges, fPFE)
+    {
+        featPointFeatEdges[fPFE].transfer(tempFeatPointFeatEdges[fPFE].shrink());
+    }
+
+    return featPointFeatEdges;
+}
+
+Foam::indexedOctree<Foam::treeDataTriSurface>::volumeType
+Foam::querySurface::insideOutside
+(
+    const scalar searchSpan2,
+    const point& pt
+) const
+{
+    if (!bb_.contains(pt))
+    {
+        return indexedOctree<treeDataTriSurface>::OUTSIDE;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, searchSpan2);
+
+    if (!pHit.hit())
+    {
+        return tree_.getVolumeType(pt);
+    }
+    else
+    {
+        return indexedOctree<treeDataTriSurface>::MIXED;
+    }
+}
+
+
+// Check if point is inside surface
+bool Foam::querySurface::inside(const point& pt) const
+{
+    if (!bb_.contains(pt))
+    {
+        return false;
+    }
+
+    return
+    (
+        tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::INSIDE
+    );
+}
+
+
+// Check if point is outside surface
+bool Foam::querySurface::outside(const point& pt) const
+{
+    if (!bb_.contains(pt))
+    {
+        return true;
+    }
+
+    return
+    (
+        tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::OUTSIDE
+    );
+}
+
+
+// Check if point is inside surface by at least dist2
+bool Foam::querySurface::wellInside(const point& pt, const scalar dist2) const
+{
+    if (!bb_.contains(pt))
+    {
+        return false;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, dist2);
+
+    if (pHit.hit())
+    {
+        return false;
+    }
+    else
+    {
+        return
+            tree_.getVolumeType(pt)
+         == indexedOctree<treeDataTriSurface>::INSIDE;
+    }
+}
+
+
+// Check if point is outside surface by at least dist2
+bool Foam::querySurface::wellOutside(const point& pt, const scalar dist2) const
+{
+    if (!bb_.contains(pt))
+    {
+        return true;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, dist2);
+
+    if (pHit.hit())
+    {
+        return false;
+    }
+    else
+    {
+        return
+            tree_.getVolumeType(pt)
+         == indexedOctree<treeDataTriSurface>::OUTSIDE;
+    }
+}
+
+
+bool Foam::querySurface::featurePoint(const label ptI) const
+{
+    if (findIndex(features().featurePoints(), ptI) >= 0)
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::querySurface::writeTreeOBJ() const
+{
+    OFstream str("tree.obj");
+    label vertI = 0;
+
+    const List<indexedOctree<treeDataTriSurface>::node>& nodes = tree_.nodes();
+
+    forAll(nodes, nodeI)
+    {
+        const indexedOctree<treeDataTriSurface>::node& nod = nodes[nodeI];
+
+        const treeBoundBox& bb = nod.bb_;
+
+        const pointField points(bb.points());
+
+        label startVertI = vertI;
+
+        forAll(points, i)
+        {
+            meshTools::writeOBJ(str, points[i]);
+            vertI++;
+        }
+
+        const edgeList edges(treeBoundBox::edges);
+
+        forAll(edges, i)
+        {
+            const edge& e = edges[i];
+
+            str << "l " << e[0]+startVertI+1 << ' ' << e[1]+startVertI+1
+                << nl;
+        }
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/querySurface.H

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    querySurface
+
+Description
+    Searchable triSurface using an octree to speed-up queries.
+
+SourceFiles
+    querySurface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef querySurface_H
+#define querySurface_H
+
+#include "triSurface.H"
+#include "treeDataTriSurface.H"
+#include "indexedOctree.H"
+#include "surfaceFeatures.H"
+#include "Random.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class querySurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class querySurface
+:
+    public triSurface
+{
+    // Private data
+
+        Random rndGen_;
+
+        // Bounding box of surface. Used for relative tolerances.
+        treeBoundBox bb_;
+
+        // Search engine on surface
+        indexedOctree<treeDataTriSurface> tree_;
+
+        surfaceFeatures sFeat_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        querySurface(const querySurface&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const querySurface&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given file name of the surface
+        querySurface
+        (
+            const fileName& surfaceFileName,
+            const scalar& includedAngle
+        );
+
+
+    // Destructor
+
+        ~querySurface();
+
+
+    // Member Functions
+
+        // Access
+
+            const treeBoundBox& bb() const
+            {
+                return bb_;
+            }
+
+            const indexedOctree<treeDataTriSurface>& tree() const
+            {
+                return tree_;
+            }
+
+            const surfaceFeatures& features() const
+            {
+                return sFeat_;
+            }
+
+
+        // Query
+
+            labelListList featurePointFeatureEdges() const;
+
+            //- Returns inside, outside or mixed (= close to surface)
+            indexedOctree<Foam::treeDataTriSurface>::volumeType insideOutside
+            (
+                const scalar searchSpan2,
+                const point& pt
+            ) const;
+
+            //- Check if point is inside surface
+            bool inside(const point& pt) const;
+
+            //- Check if point is outside surface
+            bool outside(const point& pt) const;
+
+            //- Check if point is inside surface by at least dist2
+            bool wellInside(const point& pt, const scalar dist2) const;
+
+            //- Check if point is outside surface by at least dist2
+            bool wellOutside(const point& pt, const scalar dist2) const;
+
+            //- Check if a point index is a feature point
+            bool featurePoint(const label ptI) const;
+
+
+        // Write
+
+            void writeTreeOBJ() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//#include "querySurfaceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/CV3DMesher/tolerances.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "CV3D.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::CV3D::tolerances::tolerances
+(
+    const IOdictionary& controlDict,
+    const scalar minCellSize,
+    const boundBox& bb
+)
+:
+    span
+    (
+        max(mag(bb.max().x()), mag(bb.min().x()))
+      + max(mag(bb.max().y()), mag(bb.min().y()))
+      + max(mag(bb.max().z()), mag(bb.min().z()))
+    ),
+    span2(Foam::sqr(span)),
+
+    minEdgeLen(readScalar(controlDict.lookup("minEdgeLenCoeff"))*minCellSize),
+    minEdgeLen2(Foam::sqr(minEdgeLen)),
+
+    ppDist(readScalar(controlDict.lookup("ppDistCoeff"))*minCellSize),
+    ppDist2(Foam::sqr(ppDist)),
+
+    edgeGroupSpacing
+    (
+        readScalar(controlDict.lookup("edgeGroupSpacingCoeff"))*minCellSize
+    ),
+
+    featurePointGuard
+    (
+        readScalar(controlDict.lookup("featurePointGuardCoeff"))*minCellSize
+    ),
+
+    featureEdgeGuard
+    (
+        readScalar(controlDict.lookup("featureEdgeGuardCoeff"))*minCellSize
+    ),
+
+    minEdgeSpacing
+    (
+        readScalar(controlDict.lookup("minEdgeSpacingCoeff"))*minCellSize
+    ),
+    maxEdgeSpacing
+    (
+        readScalar(controlDict.lookup("maxEdgeSpacingCoeff"))*minCellSize
+    )
+{}
+
+// ************************************************************************* //