Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C

@@ -31,12 +31,6 @@ License
 #include "labelIOField.H"
 #include "pointConversion.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Triangulation>
@@ -61,6 +55,8 @@ Foam::DelaunayMesh<Triangulation>::DelaunayMesh
     cellCount_(0),
     runTime_(runTime)
 {
+    Info<< "Reading " << meshName << " from " << runTime.timeName() << endl;
+
     pointIOField pts
     (
         IOobject
@@ -74,6 +70,8 @@ Foam::DelaunayMesh<Triangulation>::DelaunayMesh
         )
     );
 
+    if (pts.headerOk())
+    {
         labelIOField types
         (
             IOobject
@@ -82,23 +80,24 @@ Foam::DelaunayMesh<Triangulation>::DelaunayMesh
                 runTime.timeName(),
                 meshName,
                 runTime,
-            IOobject::READ_IF_PRESENT,
+                IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
         );
 
-    labelIOField indices
+// Do not read in indices
-    (
+//        labelIOField indices
-        IOobject
+//        (
-        (
+//            IOobject
-            "indices",
+//            (
-            runTime.timeName(),
+//                "indices",
-            meshName,
+//                runTime.timeName(),
-            runTime,
+//                meshName,
-            IOobject::READ_IF_PRESENT,
+//                runTime,
-            IOobject::NO_WRITE
+//                IOobject::MUST_READ,
-        )
+//                IOobject::NO_WRITE
-    );
+//            )
+//        );
 
         labelIOField processorIndices
         (
@@ -108,49 +107,34 @@ Foam::DelaunayMesh<Triangulation>::DelaunayMesh
                 runTime.timeName(),
                 meshName,
                 runTime,
-            IOobject::READ_IF_PRESENT,
+                IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
         );
 
-    if (pts.headerOk())
+        List<Vb> pointsToInsert(pts.size());
-    {
-        forAll(pts, ptI)
-        {
-            Vertex_handle vh = this->insert(toPoint<Point>(pts[ptI]));
 
-            if (indices.headerOk())
+        forAll(pointsToInsert, pI)
         {
-                vh->index() = indices[ptI];
+            pointsToInsert[pI] =
-                vertexCount_++;
+                Vb
-            }
-            else
-            {
-                vh->index() = getNewVertexIndex();
-            }
-
-            if (processorIndices.headerOk())
-            {
-                vh->procIndex() = processorIndices[ptI];
-            }
-            else
-            {
-                vh->procIndex() = Pstream::myProcNo();
-            }
-
-            if (types.headerOk())
-            {
-                vh->type() =
-                    static_cast<Foam::indexedVertexEnum::vertexType>
                 (
-                        types[ptI]
+                    toPoint(pts[pI]),
+                    pI,
+                    static_cast<indexedVertexEnum::vertexType>(types[pI]),
+                    processorIndices[pI]
                 );
         }
-            else
+
-            {
+        rangeInsertWithInfo
-                vh->type() = Vb::vtUnassigned;
+        (
-            }
+            pointsToInsert.begin(),
-        }
+            pointsToInsert.end(),
+            false,
+            false
+        );
+
+        vertexCount_ = Triangulation::number_of_vertices();
     }
 }
 
@@ -219,7 +203,7 @@ void Foam::DelaunayMesh<Triangulation>::insertPoints(const List<Vb>& vertices)
     (
         vertices.begin(),
         vertices.end(),
-        true
+        false
     );
 }
 
@@ -288,7 +272,8 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 (
     PointIterator begin,
     PointIterator end,
-    bool printErrors
+    bool printErrors,
+    bool reIndex
 )
 {
     typedef DynamicList
@@ -347,8 +332,15 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
             }
         }
         else
+        {
+            if (reIndex)
             {
                 hint->index() = getNewVertexIndex();
+            }
+            else
+            {
+                hint->index() = vert.index();
+            }
             hint->type() = vert.type();
             hint->procIndex() = vert.procIndex();
             hint->targetCellSize() = vert.targetCellSize();
@@ -358,15 +350,6 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 }
 
 
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "DelaunayMeshIO.C"

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H

@@ -90,6 +90,7 @@ public:
             FixedList<label, 2>::Hash<>
         > labelTolabelPairHashTable;
 
+
 private:
 
     // Private data
@@ -184,30 +185,49 @@ public:
 
     // Member Functions
 
+        // Access
+
+            //- Return a reference to the Time object
             inline const Time& time() const;
 
+
+        // Check
+
+            //- Write the cpuTime to screen
             inline void timeCheck
             (
                 const string& description,
                 const bool check = true
             ) const;
 
-        inline label getNewVertexIndex() const;
 
+        // Indexing functions
+
+            //- Create a new unique cell index and return
             inline label getNewCellIndex() const;
 
+            //- Create a new unique vertex index and return
+            inline label getNewVertexIndex() const;
+
+            //- Return the cell count (the next unique cell index)
             inline label cellCount() const;
 
-        inline void resetCellCount();
+            //- Return the vertex count (the next unique vertex index)
-
             inline label vertexCount() const;
 
+            //- Set the cell count to zero
+            inline void resetCellCount();
+
+            //- Set the vertex count to zero
             inline void resetVertexCount();
 
 
-        //- Remove the entire triangulation
+        // Triangulation manipulation functions
+
+            //- Clear the entire triangulation
             void reset();
 
+            //- Insert the list of vertices (calls rangeInsertWithInfo)
             void insertPoints(const List<Vb>& vertices);
 
             //- Function inserting points into a triangulation and setting the
@@ -221,14 +241,17 @@ public:
             (
                 PointIterator begin,
                 PointIterator end,
-            bool printErrors = true
+                bool printErrors = false,
+                bool reIndex = true
             );
 
 
-    // Queries
+        // Write
 
+            //- Write mesh statistics to stream
             void printInfo(Ostream& os) const;
 
+            //- Write vertex statistics in the form of a table to stream
             void printVertexInfo(Ostream& os) const;
 
             //- Create an fvMesh from the triangulation.

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshI.H

@@ -25,17 +25,6 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 template<class Triangulation>
 inline const Foam::Time& Foam::DelaunayMesh<Triangulation>::time() const
 {
@@ -70,23 +59,6 @@ void Foam::DelaunayMesh<Triangulation>::timeCheck
 }
 
 
-template<class Triangulation>
-inline Foam::label Foam::DelaunayMesh<Triangulation>::getNewVertexIndex() const
-{
-    label id = vertexCount_++;
-
-    if (id == labelMax)
-    {
-        WarningIn
-        (
-            "Foam::DelaunayMesh<Triangulation>::getNewVertexIndex() const"
-        )   << "Vertex counter has overflowed." << endl;
-    }
-
-    return id;
-}
-
-
 template<class Triangulation>
 inline Foam::label Foam::DelaunayMesh<Triangulation>::getNewCellIndex() const
 {
@@ -105,16 +77,26 @@ inline Foam::label Foam::DelaunayMesh<Triangulation>::getNewCellIndex() const
 
 
 template<class Triangulation>
-Foam::label Foam::DelaunayMesh<Triangulation>::cellCount() const
+inline Foam::label Foam::DelaunayMesh<Triangulation>::getNewVertexIndex() const
 {
-    return cellCount_;
+    label id = vertexCount_++;
+
+    if (id == labelMax)
+    {
+        WarningIn
+        (
+            "Foam::DelaunayMesh<Triangulation>::getNewVertexIndex() const"
+        )   << "Vertex counter has overflowed." << endl;
+    }
+
+    return id;
 }
 
 
 template<class Triangulation>
-void Foam::DelaunayMesh<Triangulation>::resetCellCount()
+Foam::label Foam::DelaunayMesh<Triangulation>::cellCount() const
 {
-    cellCount_ = 0;
+    return cellCount_;
 }
 
 
@@ -125,6 +107,13 @@ Foam::label Foam::DelaunayMesh<Triangulation>::vertexCount() const
 }
 
 
+template<class Triangulation>
+void Foam::DelaunayMesh<Triangulation>::resetCellCount()
+{
+    cellCount_ = 0;
+}
+
+
 template<class Triangulation>
 void Foam::DelaunayMesh<Triangulation>::resetVertexCount()
 {
@@ -132,22 +121,4 @@ void Foam::DelaunayMesh<Triangulation>::resetVertexCount()
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMeshIO.C

@@ -146,7 +146,8 @@ void Foam::DelaunayMesh<Triangulation>::printInfo(Ostream& os) const
         ++vit
     )
     {
-        if (!vit->farPoint())
+        // Only internal or boundary vertices have a size
+        if (vit->internalOrBoundaryPoint())
         {
             minSize = min(vit->targetCellSize(), minSize);
             maxSize = max(vit->targetCellSize(), maxSize);
@@ -361,19 +362,19 @@ Foam::DelaunayMesh<Triangulation>::createMesh
     // Calculate pts and a map of point index to location in pts.
     label vertI = 0;
 
-    labelIOField indices
+//    labelIOField indices
-    (
+//    (
-        IOobject
+//        IOobject
-        (
+//        (
-            "indices",
+//            "indices",
-            time().timeName(),
+//            time().timeName(),
-            name,
+//            name,
-            time(),
+//            time(),
-            IOobject::NO_READ,
+//            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
+//            IOobject::AUTO_WRITE
-        ),
+//        ),
-        Triangulation::number_of_vertices()
+//        Triangulation::number_of_vertices()
-    );
+//    );
 
     labelIOField types
     (
@@ -414,7 +415,7 @@ Foam::DelaunayMesh<Triangulation>::createMesh
         {
             vertexMap(labelPair(vit->index(), vit->procIndex())) = vertI;
             points[vertI] = topoint(vit->point());
-            indices[vertI] = vit->index();
+//            indices[vertI] = vit->index();
             types[vertI] = static_cast<label>(vit->type());
             processorIndices[vertI] = vit->procIndex();
             vertI++;
@@ -422,7 +423,7 @@ Foam::DelaunayMesh<Triangulation>::createMesh
     }
 
     points.setSize(vertI);
-    indices.setSize(vertI);
+//    indices.setSize(vertI);
     types.setSize(vertI);
     processorIndices.setSize(vertI);
 
@@ -477,8 +478,8 @@ Foam::DelaunayMesh<Triangulation>::createMesh
         bool c1Real = false;
         if
         (
-            !c1->hasFarPoint()
+            !Triangulation::is_infinite(c1)
-         && !Triangulation::is_infinite(c1)
+         && !c1->hasFarPoint()
          && c1->real()
         )
         {
@@ -490,8 +491,8 @@ Foam::DelaunayMesh<Triangulation>::createMesh
         bool c2Real = false;
         if
         (
-            !c2->hasFarPoint()
+            !Triangulation::is_infinite(c2)
-         && !Triangulation::is_infinite(c2)
+         && !c2->hasFarPoint()
          && c2->real()
         )
         {
@@ -637,7 +638,7 @@ Foam::DelaunayMesh<Triangulation>::createMesh
 
     if (writeDelaunayData)
     {
-        indices.write();
+//        indices.write();
         types.write();
         processorIndices.write();
     }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C

@@ -747,6 +747,44 @@ void Foam::DistributedDelaunayMesh<Triangulation>::sync
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class Triangulation>
+Foam::scalar
+Foam::DistributedDelaunayMesh<Triangulation>::calculateLoadUnbalance() const
+{
+    label nRealVertices = 0;
+
+    for
+    (
+        Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
+        vit != Triangulation::finite_vertices_end();
+        ++vit
+    )
+    {
+        // Only store real vertices that are not feature vertices
+        if (vit->real() && !vit->featurePoint())
+        {
+            nRealVertices++;
+        }
+    }
+
+    scalar globalNRealVertices = returnReduce
+    (
+        nRealVertices,
+        sumOp<label>()
+    );
+
+    scalar unbalance = returnReduce
+    (
+        mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
+        maxOp<scalar>()
+    );
+
+    Info<< "    Processor unbalance " << unbalance << endl;
+
+    return unbalance;
+}
+
+
 template<class Triangulation>
 bool Foam::DistributedDelaunayMesh<Triangulation>::distribute
 (
@@ -776,7 +814,8 @@ template<class Triangulation>
 Foam::autoPtr<Foam::mapDistribute>
 Foam::DistributedDelaunayMesh<Triangulation>::distribute
 (
-    const backgroundMeshDecomposition& decomposition
+    const backgroundMeshDecomposition& decomposition,
+    List<Foam::point>& points
 )
 {
     if (!Pstream::parRun())
@@ -786,21 +825,6 @@ Foam::DistributedDelaunayMesh<Triangulation>::distribute
 
     distributeBoundBoxes(decomposition.procBounds());
 
-    DynamicList<point> points(Triangulation::number_of_vertices());
-
-    for
-    (
-        Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
-        vit != Triangulation::finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->real())
-        {
-            points.append(topoint(vit->point()));
-        }
-    }
-
     autoPtr<mapDistribute> mapDist = decomposition.distributePoints(points);
 
     return mapDist;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.H

@@ -153,6 +153,8 @@ public:
         //- Use DelaunayMesh timeCheck function
         using DelaunayMesh<Triangulation>::timeCheck;
 
+        scalar calculateLoadUnbalance() const;
+
 
     // Member Functions
 
@@ -164,7 +166,8 @@ public:
 
         autoPtr<mapDistribute> distribute
         (
-            const backgroundMeshDecomposition& decomposition
+            const backgroundMeshDecomposition& decomposition,
+            List<Foam::point>& points
         );
 
         //- Refer vertices so that the processor interfaces are consistent

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshTools.C

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "DelaunayMeshTools.H"
+#include "meshTools.H"
+#include "OFstream.H"
+#include "pointConversion.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::DelaunayMeshTools::writeOBJ
+(
+    const fileName& fName,
+    const List<Foam::point>& points
+)
+{
+    if (points.size())
+    {
+        OFstream str(fName);
+
+        Pout<< nl
+            << "Writing " << points.size() << " points from pointList to "
+            << str.name() << endl;
+
+        forAll(points, p)
+        {
+            meshTools::writeOBJ(str, points[p]);
+        }
+    }
+}
+
+
+void Foam::DelaunayMeshTools::writeOBJ
+(
+    const fileName& fName,
+    const List<Vb>& points
+)
+{
+    if (points.size())
+    {
+        OFstream str(fName);
+
+        Pout<< nl
+            << "Writing " << points.size() << " points from pointList to "
+            << str.name() << endl;
+
+        forAll(points, p)
+        {
+            meshTools::writeOBJ(str, topoint(points[p].point()));
+        }
+    }
+}
+
+
+void Foam::DelaunayMeshTools::writeObjMesh
+(
+    const fileName& fName,
+    const pointField& points,
+    const faceList& faces
+)
+{
+    OFstream str(fName);
+
+    Pout<< nl
+        << "Writing points and faces to " << str.name() << endl;
+
+    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;
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshTools.H

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::DelaunayMeshTools
+
+Description
+    Collection of functions for operating on a Delaunay mesh. Includes:
+
+      - Functions for writing to an OBJ file
+      - Functions for extracting fields from the Delaunay triangulation
+
+SourceFiles
+    DelaunayMeshToolsI.H
+    DelaunayMeshTools.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DelaunayMeshTools_H
+#define DelaunayMeshTools_H
+
+#include "fileName.H"
+#include "List.H"
+#include "point.H"
+#include "CGALTriangulation3Ddefs.H"
+#include "indexedVertexEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Namespace DelaunayMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace DelaunayMeshTools
+{
+
+// OBJ writing
+
+    //- Write list of points to file
+    void writeOBJ(const fileName& fName, const List<Foam::point>& points);
+
+    //- Write list of points to file
+    void writeOBJ(const fileName& fName, const List<Vb>& points);
+
+    //- Write an OBJ mesh consisting of points and faces
+    void writeObjMesh
+    (
+        const fileName& fName,
+        const pointField& points,
+        const faceList& faces
+    );
+
+    //- Write Delaunay points in the range between (and including)
+    //  type startPointType and endPointType to an OBJ file
+    template<typename Triangulation>
+    void writeOBJ
+    (
+        const fileName& fName,
+        const Triangulation& t,
+        const indexedVertexEnum::vertexType startPointType,
+        const indexedVertexEnum::vertexType endPointType
+    );
+
+    //- Write Delaunay points of type pointType to .obj file
+    template<typename Triangulation>
+    void writeOBJ
+    (
+        const fileName& fName,
+        const Triangulation& t,
+        const indexedVertexEnum::vertexType pointType
+    );
+
+    //- Write the fixed Delaunay points to an OBJ file
+    template<typename Triangulation>
+    void writeFixedPoints(const fileName& fName, const Triangulation& t);
+
+    //- Write the boundary Delaunay points to an OBJ file
+    template<typename Triangulation>
+    void writeBoundaryPoints(const fileName& fName, const Triangulation& t);
+
+    //- Write the processor interface to an OBJ file
+    template<typename Triangulation>
+    void writeProcessorInterface
+    (
+        const fileName& fName,
+        const Triangulation& t,
+        const faceList& faces
+    );
+
+    //- Write the internal Delaunay vertices of the tessellation as a
+    //  pointField that may be used to restart the meshing process
+    template<typename Triangulation>
+    void writeInternalDelaunayVertices
+    (
+        const fileName& instance,
+        const Triangulation& t
+    );
+
+    //- Draws a tet cell to an output stream. The offset is supplied as the tet
+    //  number to be drawn.
+    template<typename CellHandle>
+    void drawDelaunayCell(Ostream& os, const CellHandle& c, label offset = 0);
+
+
+// Field extraction
+
+    //- Extract all points in vertex-index order
+    template<typename Triangulation>
+    tmp<pointField> allPoints(const Triangulation& t);
+
+
+} // End namespace DelaunayMeshTools
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DelaunayMeshToolsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C

@@ -0,0 +1,305 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "DelaunayMeshTools.H"
+#include "meshTools.H"
+#include "OFstream.H"
+#include "pointConversion.H"
+#include "pointIOField.H"
+#include "indexedVertexOps.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<typename Triangulation>
+void Foam::DelaunayMeshTools::writeOBJ
+(
+    const fileName& fName,
+    const Triangulation& t,
+    const indexedVertexEnum::vertexType startPointType,
+    const indexedVertexEnum::vertexType endPointType
+)
+{
+    OFstream str(fName);
+
+    Pout<< nl
+        << "Writing points of types:" << nl;
+
+    forAllConstIter
+    (
+        HashTable<int>,
+        indexedVertexEnum::vertexTypeNames_,
+        iter
+    )
+    {
+        if (iter() >= startPointType && iter() <= endPointType)
+        {
+            Pout<< "    " << iter.key() << nl;
+        }
+    }
+
+    Pout<< "to " << str.name() << endl;
+
+    for
+    (
+        typename Triangulation::Finite_vertices_iterator vit =
+            t.finite_vertices_begin();
+        vit != t.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->type() >= startPointType && vit->type() <= endPointType)
+        {
+            meshTools::writeOBJ(str, topoint(vit->point()));
+        }
+    }
+}
+
+
+template<typename Triangulation>
+void Foam::DelaunayMeshTools::writeOBJ
+(
+    const fileName& fName,
+    const Triangulation& t,
+    const indexedVertexEnum::vertexType pointType
+)
+{
+    writeOBJ(fName, t, pointType, pointType);
+}
+
+
+template<typename Triangulation>
+void Foam::DelaunayMeshTools::writeFixedPoints
+(
+    const fileName& fName,
+    const Triangulation& t
+)
+{
+    OFstream str(fName);
+
+    Pout<< nl
+        << "Writing fixed points to " << str.name() << endl;
+
+    for
+    (
+        typename Triangulation::Finite_vertices_iterator vit =
+            t.finite_vertices_begin();
+        vit != t.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->fixed())
+        {
+            meshTools::writeOBJ(str, topoint(vit->point()));
+        }
+    }
+}
+
+
+template<typename Triangulation>
+void Foam::DelaunayMeshTools::writeBoundaryPoints
+(
+    const fileName& fName,
+    const Triangulation& t
+)
+{
+    OFstream str(fName);
+
+    Pout<< nl
+        << "Writing boundary points to " << str.name() << endl;
+
+    for
+    (
+        typename Triangulation::Finite_vertices_iterator vit =
+            t.finite_vertices_begin();
+        vit != t.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!vit->internalPoint())
+        {
+            meshTools::writeOBJ(str, topoint(vit->point()));
+        }
+    }
+}
+
+
+template<typename Triangulation>
+void Foam::DelaunayMeshTools::writeProcessorInterface
+(
+    const fileName& fName,
+    const Triangulation& t,
+    const faceList& faces
+)
+{
+    OFstream str(fName);
+
+    pointField points(t.number_of_finite_cells(), point::max);
+
+    for
+    (
+        typename Triangulation::Finite_cells_iterator cit =
+            t.finite_cells_begin();
+        cit != t.finite_cells_end();
+        ++cit
+    )
+    {
+        if (!cit->hasFarPoint() && !t.is_infinite(cit))
+        {
+            points[cit->cellIndex()] = cit->dual();
+        }
+    }
+
+    meshTools::writeOBJ(str, faces, points);
+}
+
+
+template<typename Triangulation>
+void Foam::DelaunayMeshTools::writeInternalDelaunayVertices
+(
+    const fileName& instance,
+    const Triangulation& t
+)
+{
+    pointField internalDelaunayVertices(t.number_of_vertices());
+
+    label vertI = 0;
+
+    for
+    (
+        typename Triangulation::Finite_vertices_iterator vit =
+            t.finite_vertices_begin();
+        vit != t.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalPoint())
+        {
+            internalDelaunayVertices[vertI++] = topoint(vit->point());
+        }
+    }
+
+    internalDelaunayVertices.setSize(vertI);
+
+    pointIOField internalDVs
+    (
+        IOobject
+        (
+            "internalDelaunayVertices",
+            instance,
+            t.time(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        internalDelaunayVertices
+    );
+
+    Info<< nl
+        << "Writing " << internalDVs.name()
+        << " to " << internalDVs.instance()
+        << endl;
+
+    internalDVs.write();
+}
+
+
+template<typename CellHandle>
+void Foam::DelaunayMeshTools::drawDelaunayCell
+(
+    Ostream& os,
+    const CellHandle& c,
+    label offset
+)
+{
+    // Supply offset as tet number
+    offset *= 4;
+
+    os  << "# cell index: " << label(c->cellIndex())
+        << " INT_MIN = " << INT_MIN
+        << endl;
+
+    os  << "# circumradius "
+        << mag(c->dual() - topoint(c->vertex(0)->point()))
+        << endl;
+
+    for (int i = 0; i < 4; i++)
+    {
+        os  << "# index / type / procIndex: "
+            << label(c->vertex(i)->index()) << " "
+            << label(c->vertex(i)->type()) << " "
+            << label(c->vertex(i)->procIndex())
+            <<
+                (
+                    CGAL::indexedVertexOps::uninitialised(c->vertex(i))
+                  ? " # This vertex is uninitialised!"
+                  : ""
+                )
+            << endl;
+
+        meshTools::writeOBJ(os, topoint(c->vertex(i)->point()));
+    }
+
+    os  << "f " << 1 + offset << " " << 3 + offset << " " << 2 + offset << nl
+        << "f " << 2 + offset << " " << 3 + offset << " " << 4 + offset << nl
+        << "f " << 1 + offset << " " << 4 + offset << " " << 3 + offset << nl
+        << "f " << 1 + offset << " " << 2 + offset << " " << 4 + offset << endl;
+
+//    os  << "# cicumcentre " << endl;
+
+//    meshTools::writeOBJ(os, c->dual());
+
+//    os  << "l " << 1 + offset << " " << 5 + offset << endl;
+}
+
+
+template<typename Triangulation>
+Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
+(
+    const Triangulation& t
+)
+{
+    tmp<pointField> tpts(new pointField(t.number_of_vertices(), point::max));
+    pointField& pts = tpts();
+
+    label nVert = 0;
+
+    for
+    (
+        typename Triangulation::Finite_vertices_iterator vit =
+            t.finite_vertices_begin();
+        vit != t.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalOrBoundaryPoint())
+        {
+            pts[nVert++] = topoint(vit->point());
+        }
+    }
+
+    return tpts;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files

@@ -1,5 +1,7 @@
 #include CGAL_FILES
 
+DelaunayMeshTools/DelaunayMeshTools.C
+
 conformalVoronoiMesh/indexedVertex/indexedVertexEnum.C
 conformalVoronoiMesh/indexedCell/indexedCellEnum.C
 
@@ -8,7 +10,10 @@ conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
 conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
 conformalVoronoiMesh/conformalVoronoiMeshIO.C
 conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
-conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
+
+conformalVoronoiMesh/featurePointConformer/pointFeatureEdgesTypes.C
+conformalVoronoiMesh/featurePointConformer/featurePointConformer.C
+conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C
 
 cvControls/cvControls.C
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/PrintTable/PrintTable.C

@@ -143,7 +143,7 @@ void Foam::PrintTable<KeyType, DataType>::print
         os  << nl << indent << tab << "# " << title_.c_str() << endl;
 
         os.width(largestKeyLength);
-        os  << indent << "# Proc No";
+        os  << indent << "# Proc";
 
         forAll(procData, procI)
         {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C

@@ -29,6 +29,7 @@ License
 #include "zeroGradientFvPatchFields.H"
 #include "Time.H"
 #include "Random.H"
+#include "pointConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1043,22 +1044,6 @@ Foam::backgroundMeshDecomposition::distribute
 }
 
 
-Foam::autoPtr<Foam::mapDistribute>
-Foam::backgroundMeshDecomposition::distributePoints
-(
-    List<point>& points
-) const
-{
-    labelList toProc(processorPosition(points));
-
-    autoPtr<mapDistribute> map(buildMap(toProc));
-
-    map().distribute(points);
-
-    return map;
-}
-
-
 bool Foam::backgroundMeshDecomposition::positionOnThisProcessor
 (
     const point& pt
@@ -1128,141 +1113,6 @@ Foam::pointIndexHit Foam::backgroundMeshDecomposition::findLineAny
 }
 
 
-Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
-(
-    const List<point>& pts
-) const
-{
-    DynamicList<label> toCandidateProc;
-    DynamicList<point> testPoints;
-    labelList ptBlockStart(pts.size(), -1);
-    labelList ptBlockSize(pts.size(), -1);
-
-    label nTotalCandidates = 0;
-
-    forAll(pts, pI)
-    {
-        const point& pt = pts[pI];
-
-        label nCandidates = 0;
-
-        forAll(allBackgroundMeshBounds_, procI)
-        {
-            if (allBackgroundMeshBounds_[procI].contains(pt))
-            {
-                toCandidateProc.append(procI);
-                testPoints.append(pt);
-
-                nCandidates++;
-            }
-        }
-
-        ptBlockStart[pI] = nTotalCandidates;
-        ptBlockSize[pI] = nCandidates;
-
-        nTotalCandidates += nCandidates;
-    }
-
-    // Needed for reverseDistribute
-    label preDistributionToCandidateProcSize = toCandidateProc.size();
-
-    autoPtr<mapDistribute> map(buildMap(toCandidateProc));
-
-    map().distribute(testPoints);
-
-    List<bool> pointOnCandidate(testPoints.size(), false);
-
-    // Test candidate points on candidate processors
-    forAll(testPoints, tPI)
-    {
-        pointOnCandidate[tPI] = positionOnThisProcessor(testPoints[tPI]);
-    }
-
-    map().reverseDistribute
-    (
-        preDistributionToCandidateProcSize,
-        pointOnCandidate
-    );
-
-    labelList ptProc(pts.size(), -1);
-
-    DynamicList<label> failedPointIndices;
-    DynamicList<point> failedPoints;
-
-    forAll(pts, pI)
-    {
-        // Extract the sub list of results for this point
-
-        SubList<bool> ptProcResults
-        (
-            pointOnCandidate,
-            ptBlockSize[pI],
-            ptBlockStart[pI]
-        );
-
-        forAll(ptProcResults, pPRI)
-        {
-            if (ptProcResults[pPRI])
-            {
-                ptProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
-
-                break;
-            }
-        }
-
-        if (ptProc[pI] < 0)
-        {
-            if (!globalBackgroundBounds_.contains(pts[pI]))
-            {
-                FatalErrorIn
-                (
-                    "Foam::labelList"
-                    "Foam::backgroundMeshDecomposition::processorPosition"
-                    "("
-                        "const List<point>&"
-                    ") const"
-                )
-                    << "The position " << pts[pI]
-                    << " is not in any part of the background mesh "
-                    << globalBackgroundBounds_ << endl
-                    << "A background mesh with a wider margin around "
-                    << "the geometry may help."
-                    << exit(FatalError);
-            }
-
-            if (debug)
-            {
-                WarningIn
-                (
-                    "Foam::labelList"
-                    "Foam::backgroundMeshDecomposition::processorPosition"
-                    "("
-                        "const List<point>&"
-                    ") const"
-                )   << "The position " << pts[pI]
-                    << " was not found in the background mesh "
-                    << globalBackgroundBounds_ << ", finding nearest."
-                    << endl;
-            }
-
-            failedPointIndices.append(pI);
-            failedPoints.append(pts[pI]);
-        }
-    }
-
-    labelList ptNearestProc(processorNearestPosition(failedPoints));
-
-    forAll(failedPoints, fPI)
-    {
-        label pI = failedPointIndices[fPI];
-
-        ptProc[pI] = ptNearestProc[fPI];
-    }
-
-    return ptProc;
-}
-
-
 Foam::labelList Foam::backgroundMeshDecomposition::processorNearestPosition
 (
     const List<point>& pts

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H

@@ -72,6 +72,7 @@ SourceFiles
 #include "indexedOctree.H"
 #include "treeDataPrimitivePatch.H"
 #include "volumeType.H"
+#include "CGALTriangulation3Ddefs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -224,7 +225,8 @@ public:
         );
 
         //- Distribute supplied the points to the appropriate processor
-        autoPtr<mapDistribute> distributePoints(List<point>& points) const;
+        template<typename PointType>
+        autoPtr<mapDistribute> distributePoints(List<PointType>& points) const;
 
         //- Is the given position inside the domain of this decomposition
         bool positionOnThisProcessor(const point& pt) const;
@@ -261,7 +263,8 @@ public:
         ) const;
 
         //- What processor is the given position on?
-        labelList processorPosition(const List<point>& pts) const;
+        template<typename PointType>
+        labelList processorPosition(const List<PointType>& pts) const;
 
         //- What is the nearest processor to the given position?
         labelList processorNearestPosition(const List<point>& pts) const;
@@ -334,6 +337,10 @@ public:
 
 #include "backgroundMeshDecompositionI.H"
 
+#ifdef NoRepository
+#   include "backgroundMeshDecompositionTemplates.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionTemplates.C

@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "backgroundMeshDecomposition.H"
+#include "pointConversion.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<typename PointType>
+Foam::autoPtr<Foam::mapDistribute>
+Foam::backgroundMeshDecomposition::distributePoints
+(
+    List<PointType>& points
+) const
+{
+    labelList toProc(processorPosition(points));
+
+    autoPtr<mapDistribute> map(buildMap(toProc));
+
+    map().distribute(points);
+
+    return map;
+}
+
+
+template<typename PointType>
+Foam::labelList Foam::backgroundMeshDecomposition::processorPosition
+(
+    const List<PointType>& pts
+) const
+{
+    DynamicList<label> toCandidateProc;
+    DynamicList<point> testPoints;
+    labelList ptBlockStart(pts.size(), -1);
+    labelList ptBlockSize(pts.size(), -1);
+
+    label nTotalCandidates = 0;
+
+    forAll(pts, pI)
+    {
+        const pointFromPoint pt = topoint(pts[pI]);
+
+        label nCandidates = 0;
+
+        forAll(allBackgroundMeshBounds_, procI)
+        {
+            if (allBackgroundMeshBounds_[procI].contains(pt))
+            {
+                toCandidateProc.append(procI);
+                testPoints.append(pt);
+
+                nCandidates++;
+            }
+        }
+
+        ptBlockStart[pI] = nTotalCandidates;
+        ptBlockSize[pI] = nCandidates;
+
+        nTotalCandidates += nCandidates;
+    }
+
+    // Needed for reverseDistribute
+    label preDistributionToCandidateProcSize = toCandidateProc.size();
+
+    autoPtr<mapDistribute> map(buildMap(toCandidateProc));
+
+    map().distribute(testPoints);
+
+    List<bool> pointOnCandidate(testPoints.size(), false);
+
+    // Test candidate points on candidate processors
+    forAll(testPoints, tPI)
+    {
+        pointOnCandidate[tPI] = positionOnThisProcessor(testPoints[tPI]);
+    }
+
+    map().reverseDistribute
+    (
+        preDistributionToCandidateProcSize,
+        pointOnCandidate
+    );
+
+    labelList ptProc(pts.size(), -1);
+
+    DynamicList<label> failedPointIndices;
+    DynamicList<point> failedPoints;
+
+    forAll(pts, pI)
+    {
+        // Extract the sub list of results for this point
+
+        SubList<bool> ptProcResults
+        (
+            pointOnCandidate,
+            ptBlockSize[pI],
+            ptBlockStart[pI]
+        );
+
+        forAll(ptProcResults, pPRI)
+        {
+            if (ptProcResults[pPRI])
+            {
+                ptProc[pI] = toCandidateProc[ptBlockStart[pI] + pPRI];
+
+                break;
+            }
+        }
+
+        if (ptProc[pI] < 0)
+        {
+            const pointFromPoint pt = topoint(pts[pI]);
+
+            if (!globalBackgroundBounds_.contains(pt))
+            {
+                FatalErrorIn
+                (
+                    "Foam::labelList"
+                    "Foam::backgroundMeshDecomposition::processorPosition"
+                    "("
+                        "const List<point>&"
+                    ") const"
+                )
+                    << "The position " << pt
+                    << " is not in any part of the background mesh "
+                    << globalBackgroundBounds_ << endl
+                    << "A background mesh with a wider margin around "
+                    << "the geometry may help."
+                    << exit(FatalError);
+            }
+
+            if (debug)
+            {
+                WarningIn
+                (
+                    "Foam::labelList"
+                    "Foam::backgroundMeshDecomposition::processorPosition"
+                    "("
+                        "const List<point>&"
+                    ") const"
+                )   << "The position " << pt
+                    << " was not found in the background mesh "
+                    << globalBackgroundBounds_ << ", finding nearest."
+                    << endl;
+            }
+
+            failedPointIndices.append(pI);
+            failedPoints.append(pt);
+        }
+    }
+
+    labelList ptNearestProc(processorNearestPosition(failedPoints));
+
+    forAll(failedPoints, fPI)
+    {
+        label pI = failedPointIndices[fPI];
+
+        ptProc[pI] = ptNearestProc[fPI];
+    }
+
+    return ptProc;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControl/cellShapeControl.C

@@ -30,6 +30,7 @@ License
 #include "cellSizeAndAlignmentControl.H"
 #include "searchableSurfaceControl.H"
 #include "cellSizeFunction.H"
+#include "indexedVertexOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -152,7 +153,7 @@ Foam::scalar Foam::cellShapeControl::cellSize(const point& pt) const
         {
             for (label pI = 0; pI < 4; ++pI)
             {
-                if (!ch->vertex(pI)->uninitialised())
+                if (!CGAL::indexedVertexOps::uninitialised(ch->vertex(pI)))
                 {
                     size = ch->vertex(pI)->targetCellSize();
                     return size;
@@ -276,7 +277,7 @@ void Foam::cellShapeControl::cellSizeAndAlignment
         {
             for (label pI = 0; pI < 4; ++pI)
             {
-                if (!ch->vertex(pI)->uninitialised())
+                if (!CGAL::indexedVertexOps::uninitialised(ch->vertex(pI)))
                 {
                     size = ch->vertex(pI)->targetCellSize();
                     alignment = ch->vertex(pI)->alignment();

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C

@@ -425,7 +425,6 @@ Foam::cellShapeControlMesh::cellShapeControlMesh(const Time& runTime)
              && alignments.size() == this->vertexCount()
             )
             {
-                label count = 0;
                 for
                 (
                     Finite_vertices_iterator vit = finite_vertices_begin();
@@ -433,9 +432,8 @@ Foam::cellShapeControlMesh::cellShapeControlMesh(const Time& runTime)
                     ++vit
                 )
                 {
-                    vit->targetCellSize() = sizes[count];
+                    vit->targetCellSize() = sizes[vit->index()];
-                    vit->alignment() = alignments[count];
+                    vit->alignment() = alignments[vit->index()];
-                    count++;
                 }
             }
             else
@@ -453,155 +451,6 @@ Foam::cellShapeControlMesh::cellShapeControlMesh(const Time& runTime)
 }
 
 
-//Foam::triangulatedMesh::triangulatedMesh
-//(
-//    const Time& runTime,
-//    const fileName& pointsFile,
-//    const fileName& sizesFile,
-//    const fileName& alignmentsFile,
-//    const scalar& defaultCellSize
-//)
-//:
-//    defaultCellSize_(defaultCellSize)
-//{
-//    Info<< "    Reading points from file     : " << pointsFile << endl;
-//
-//    pointIOField points
-//    (
-//        IOobject
-//        (
-//            pointsFile,
-//            runTime.constant(),
-//            runTime,
-//            IOobject::MUST_READ,
-//            IOobject::NO_WRITE
-//        )
-//    );
-//
-//    Info<< "    Reading sizes from file      : " << sizesFile << endl;
-//
-//    scalarIOField sizes
-//    (
-//        IOobject
-//        (
-//            sizesFile,
-//            runTime.constant(),
-//            runTime,
-//            IOobject::MUST_READ,
-//            IOobject::NO_WRITE
-//        )
-//    );
-//
-//    Info<< "    Reading alignments from file : " << alignmentsFile << endl;
-//
-//    tensorIOField alignments
-//    (
-//        IOobject
-//        (
-//            alignmentsFile,
-//            runTime.constant(),
-//            runTime,
-//            IOobject::MUST_READ,
-//            IOobject::NO_WRITE
-//        )
-//    );
-//
-//    Info<< "    Number of points : " << points.size() << endl;
-//    Info<< "    Minimum size     : " << min(sizes) << endl;
-//    Info<< "    Average size     : " << average(sizes) << endl;
-//    Info<< "    Maximum size     : " << max(sizes) << endl;
-//
-//    forAll(points, pI)
-//    {
-//        size_t nVert = number_of_vertices();
-//
-//        Vertex_handle v = insert
-//        (
-//            Point(points[pI].x(), points[pI].y(), points[pI].z())
-//        );
-//
-//        if (number_of_vertices() == nVert)
-//        {
-//            Info<< "    Failed to insert point : " << points[pI] << endl;
-//        }
-//
-//        v->targetCellSize() = sizes[pI];
-//
-//        const tensor& alignment = alignments[pI];
-//
-//
-//
-//        v->alignment() = alignment;
-//    }
-//
-////    scalar factor = 1.0;
-////    label maxIteration = 1;
-////
-////    for (label iteration = 0; iteration < maxIteration; ++iteration)
-////    {
-////        Info<< "Iteration : " << iteration << endl;
-////
-////        label nRefined = refineTriangulation(factor);
-////
-////        Info<< "    Number of cells refined in refinement iteration : "
-////            << nRefined << nl << endl;
-////
-////        if (nRefined <= 0 && iteration != 0)
-////        {
-////            break;
-////        }
-////
-////        factor *= 1.5;
-////    }
-//
-//    //writeRefinementTriangulation();
-//}
-
-
-//Foam::triangulatedMesh::triangulatedMesh
-//(
-//    const Time& runTime,
-//    const DynamicList<Foam::point>& points,
-//    const DynamicList<scalar>& sizes,
-//    const DynamicList<tensor>& alignments,
-//    const scalar& defaultCellSize
-//)
-//:
-//    defaultCellSize_(defaultCellSize)
-//{
-//    forAll(points, pI)
-//    {
-//        size_t nVert = number_of_vertices();
-//
-//        Vertex_handle v = insert
-//        (
-//            Point(points[pI].x(), points[pI].y(), points[pI].z())
-//        );
-//
-//        if (number_of_vertices() == nVert)
-//        {
-//            Info<< "Failed to insert point : " << points[pI] << endl;
-//        }
-//
-//        v->targetCellSize() = sizes[pI];
-//
-//        v->alignment() = alignments[pI];
-//    }
-//
-//    //writeRefinementTriangulation();
-//
-//    Info<< nl << "Refinement triangulation information: " << endl;
-//    Info<< "    Number of vertices: " << label(number_of_vertices()) << endl;
-//    Info<< "    Number of cells   : "
-//        << label(number_of_finite_cells()) << endl;
-//    Info<< "    Number of faces   : "
-//        << label(number_of_finite_facets()) << endl;
-//    Info<< "    Number of edges   : "
-//        << label(number_of_finite_edges()) << endl;
-//    Info<< "    Dimensionality    : " << label(dimension()) << nl << endl;
-//}
-
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cellShapeControlMesh::~cellShapeControlMesh()
@@ -619,7 +468,7 @@ void Foam::cellShapeControlMesh::barycentricCoords
 {
     // Use the previous cell handle as a hint on where to start searching
     // Giving a hint causes strange errors...
-    ch = locate(toPoint<Point>(pt));
+    ch = locate(toPoint(pt));
 
     if (dimension() > 2 && !is_infinite(ch))
     {
@@ -666,14 +515,6 @@ void Foam::cellShapeControlMesh::distribute
     const backgroundMeshDecomposition& decomposition
 )
 {
-    if (!Pstream::parRun())
-    {
-        return;
-    }
-
-    autoPtr<mapDistribute> mapDist =
-        DistributedDelaunayMesh<CellSizeDelaunay>::distribute(decomposition);
-
     DynamicList<Foam::point> points(number_of_vertices());
     DynamicList<scalar> sizes(number_of_vertices());
     DynamicList<tensor> alignments(number_of_vertices());
@@ -711,7 +552,13 @@ void Foam::cellShapeControlMesh::distribute
         }
     }
 
-    mapDist().distribute(points);
+    autoPtr<mapDistribute> mapDist =
+        DistributedDelaunayMesh<CellSizeDelaunay>::distribute
+        (
+            decomposition,
+            points
+        );
+
     mapDist().distribute(sizes);
     mapDist().distribute(alignments);
 
@@ -735,7 +582,7 @@ void Foam::cellShapeControlMesh::distribute
         (
             Vb
             (
-                toPoint<Point>(points[pI]),
+                toPoint(points[pI]),
                 -1,
                 Vb::vtInternal,
                 Pstream::myProcNo()

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C

@@ -328,7 +328,7 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
 
     cellSizeFunctions_.reorder(invertedFunctionPriorities);
 
-    Info<< nl << "There are " << cellSizeFunctions_.size()
+    Info<< nl << indent << "There are " << cellSizeFunctions_.size()
         << " region control functions" << endl;
 }
 
@@ -487,24 +487,25 @@ void Foam::searchableSurfaceControl::cellSizeFunctionVertices
 
     const scalar nearFeatDistSqrCoeff = 1e-8;
 
+    pointField ptField(1, vector::min);
+    scalarField distField(1, nearFeatDistSqrCoeff);
+    List<pointIndexHit> infoList(1, pointIndexHit());
+
+    vectorField normals(1);
+    labelList region(1, -1);
+
     forAll(points, pI)
     {
         // Is the point in the extendedFeatureEdgeMesh? If so get the
         // point normal, otherwise get the surface normal from
         // searchableSurface
-
+        ptField[0] = points[pI];
-        pointField ptField(1, points[pI]);
-        scalarField distField(1, nearFeatDistSqrCoeff);
-        List<pointIndexHit> infoList(1, pointIndexHit());
 
         searchableSurface_.findNearest(ptField, distField, infoList);
 
         if (infoList[0].hit())
         {
-            vectorField normals(1);
             searchableSurface_.getNormal(infoList, normals);
-
-            labelList region(1, -1);
             searchableSurface_.getRegion(infoList, region);
 
             const cellSizeFunction& sizeFunc =

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C

@@ -284,14 +284,9 @@ void Foam::controlMeshRefinement::initialMeshPopulation
         {
             vertices[vI] = Vb(pts[vI], Vb::vtInternalNearBoundary);
 
-//            Info<< "Find size of vertex " << vI << endl;
-
             label maxPriority = -1;
             scalar size = sizeControls_.cellSize(pts[vI], maxPriority);
 
-//            Info<< "    Size = " << size << ", priority = " << maxPriority
-//                << endl;
-
             if (maxPriority > controlFunction.maxPriority())
             {
                 vertices[vI].targetCellSize() = max
@@ -763,7 +758,7 @@ Foam::label Foam::controlMeshRefinement::refineMesh
             (
                 Vb
                 (
-                    toPoint<cellShapeControlMesh::Point>(pt),
+                    toPoint(pt),
                     Vb::vtInternal
                 )
             );

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C

@@ -58,9 +58,14 @@ Foam::fieldFromFile::fieldFromFile
         surface,
         defaultCellSize
     ),
+    coeffsDict_(cellSizeCalcTypeDict.subDict(typeName + "Coeffs")),
     fileName_
     (
         cellSizeCalcTypeDict.subDict(typeName + "Coeffs").lookup("fieldFile")
+    ),
+    cellSizeMultipleCoeff_
+    (
+        coeffsDict_.lookupOrDefault<scalar>("cellSizeMultipleCoeff", 1)
     )
 {}
 
@@ -90,6 +95,8 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::fieldFromFile::load()
         )
     );
 
+    pointCellSize() *= cellSizeMultipleCoeff_;
+
     return pointCellSize;
 }
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.H

@@ -58,10 +58,16 @@ private:
 
     // Private data
 
+        //- Dictionary of coefficients for automatic cell sizing
+        const dictionary& coeffsDict_;
+
         //- Name of the triSurfaceScalarField file to load in. Must be in
         //  constant/triSurface
         const word fileName_;
 
+        //- Multiply all the point sizes by this value
+        const scalar cellSizeMultipleCoeff_;
+
 
 public:
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C

@@ -109,7 +109,7 @@ Foam::scalar Foam::nonUniformField::interpolate
     const face& faceHitByPt = surfaceTriMesh_.triSurface::operator[](index);
 
     const pointField& pts = surfaceTriMesh_.points();
-    const Map<label>& pMap = surfaceTriMesh_.meshPointMap();
+//    const Map<label>& pMap = surfaceTriMesh_.meshPointMap();
 
     triPointRef tri
     (
@@ -122,9 +122,12 @@ Foam::scalar Foam::nonUniformField::interpolate
 
     tri.barycentric(pt, bary);
 
-    return pointCellSize_[pMap[faceHitByPt[0]]]*bary[0]
+//    return pointCellSize_[pMap[faceHitByPt[0]]]*bary[0]
-         + pointCellSize_[pMap[faceHitByPt[1]]]*bary[1]
+//         + pointCellSize_[pMap[faceHitByPt[1]]]*bary[1]
-         + pointCellSize_[pMap[faceHitByPt[2]]]*bary[2];
+//         + pointCellSize_[pMap[faceHitByPt[2]]]*bary[2];
+    return pointCellSize_[faceHitByPt[0]]*bary[0]
+         + pointCellSize_[faceHitByPt[1]]*bary[1]
+         + pointCellSize_[faceHitByPt[2]]*bary[2];
 }
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H

@@ -38,12 +38,14 @@ Description
 #ifdef CGAL_INEXACT
 
     // Fast kernel using a double as the storage type but the triangulation may
-    // fail. Adding robust circumcentre traits
+    // fail. Adding robust circumcentre traits.
     #include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
     typedef CGAL::Exact_predicates_inexact_constructions_kernel baseK;
+
 //    #include <CGAL/Robust_circumcenter_traits_3.h>
-    #include <CGAL/Robust_circumcenter_filtered_traits_3.h>
 //    typedef CGAL::Robust_circumcenter_traits_3<baseK>           K;
+
+    #include <CGAL/Robust_circumcenter_filtered_traits_3.h>
     typedef CGAL::Robust_circumcenter_filtered_traits_3<baseK>  K;
 
 #else
@@ -55,6 +57,7 @@ Description
 
 #endif
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H

@@ -53,6 +53,7 @@ typedef CGAL::Triangulation_data_structure_3<Vb, Cb>           Tds;
 typedef CGAL::Delaunay_triangulation_3<K, Tds, CompactLocator> Delaunay;
 typedef CGAL::Delaunay_triangulation_3<K, Tds, FastLocator>    CellSizeDelaunay;
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -34,6 +34,8 @@ License
 #include "controlMeshRefinement.H"
 #include "smoothAlignmentSolver.H"
 #include "OBJstream.H"
+#include "indexedVertexOps.H"
+#include "DelaunayMeshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,8 +44,6 @@ namespace Foam
 defineTypeNameAndDebug(conformalVoronoiMesh, 0);
 }
 
-const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3;
-
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
@@ -124,12 +124,13 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
             decomposition_().distributePoints(transferPoints)
         );
 
+        transferPoints.clear();
+
         map().distribute(points);
     }
 
     label nVert = number_of_vertices();
 
-    // using the range insert (faster than inserting points one by one)
     insert(points.begin(), points.end());
 
     label nInserted(number_of_vertices() - nVert);
@@ -152,7 +153,7 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
         ++vit
     )
     {
-        if (vit->uninitialised())
+        if (CGAL::indexedVertexOps::uninitialised(vit))
         {
             vit->index() = getNewVertexIndex();
             vit->type() = Vb::vtInternal;
@@ -169,24 +170,7 @@ void Foam::conformalVoronoiMesh::insertPoints
 {
     if (Pstream::parRun() && distribute)
     {
-        const label preDistributionSize = vertices.size();
+        decomposition_().distributePoints(vertices);
-
-        List<Foam::point> pts(preDistributionSize);
-
-        forAll(vertices, vI)
-        {
-            const Foam::point& pt = topoint(vertices[vI].point());
-
-            pts[vI] = pt;
-        }
-
-        // Distribute points to their appropriate processor
-        autoPtr<mapDistribute> map
-        (
-            decomposition_().distributePoints(pts)
-        );
-
-        map().distribute(vertices);
 
         forAll(vertices, vI)
         {
@@ -196,12 +180,7 @@ void Foam::conformalVoronoiMesh::insertPoints
 
     label preReinsertionSize(number_of_vertices());
 
-    rangeInsertWithInfo
+    this->DelaunayMesh<Delaunay>::insertPoints(vertices);
-    (
-        vertices.begin(),
-        vertices.end(),
-        false
-    );
 
     const label nReinserted = returnReduce
     (
@@ -288,7 +267,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 
     if (foamyHexMeshControls().objOutput() && fName != fileName::null)
     {
-        writePoints(fName, pts);
+        DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
     }
 }
 
@@ -328,7 +307,7 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
 
     if (foamyHexMeshControls().objOutput() && fName != fileName::null)
     {
-        writePoints(fName, pts);
+        DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
     }
 }
 
@@ -381,9 +360,10 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
 
     if (foamyHexMeshControls().objOutput())
     {
-        writePoints
+        DelaunayMeshTools::writeOBJ
         (
-            "initialPoints.obj",
+            time().path()/"initialPoints.obj",
+            *this,
             Foam::indexedVertexEnum::vtInternal
         );
     }
@@ -397,9 +377,6 @@ void Foam::conformalVoronoiMesh::distribute()
         return;
     }
 
-    autoPtr<mapDistribute> mapDist =
-        DistributedDelaunayMesh<Delaunay>::distribute(decomposition_());
-
     DynamicList<Foam::point> points(number_of_vertices());
     DynamicList<Foam::indexedVertexEnum::vertexType> types
     (
@@ -424,7 +401,9 @@ void Foam::conformalVoronoiMesh::distribute()
         }
     }
 
-    mapDist().distribute(points);
+    autoPtr<mapDistribute> mapDist =
+        DistributedDelaunayMesh<Delaunay>::distribute(decomposition_(), points);
+
     mapDist().distribute(types);
     mapDist().distribute(sizes);
     mapDist().distribute(alignments);
@@ -444,7 +423,7 @@ void Foam::conformalVoronoiMesh::distribute()
         (
             Vb
             (
-                toPoint<Point>(points[pI]),
+                toPoint(points[pI]),
                 -1,
                 types[pI],
                 Pstream::myProcNo()
@@ -549,7 +528,7 @@ void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
         << (cellSizeMesh.is_valid() ? "valid" : "not valid!" )
         << endl;
 
-    if (!Pstream::parRun() && foamyHexMeshControls().objOutput())
+    if (foamyHexMeshControls().writeCellShapeControlMesh())
     {
         //cellSizeMesh.writeTriangulation();
         cellSizeMesh.write();
@@ -583,7 +562,7 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
     {
         if (vit->internalOrBoundaryPoint())
         {
-            pointFromPoint pt = topoint(vit->point());
+            const pointFromPoint pt = topoint(vit->point());
 
             cellShapeControls().cellSizeAndAlignment
             (
@@ -625,18 +604,19 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
             Vertex_handle vA = c->vertex(eit->second);
             Vertex_handle vB = c->vertex(eit->third);
 
-            drawDelaunayCell(Pout, cc1);
+//            DelaunayMeshTools::drawDelaunayCell(Pout, cc1);
-            drawDelaunayCell(Pout, cc2);
+//            DelaunayMeshTools::drawDelaunayCell(Pout, cc2);
 
-            FatalErrorIn("Foam::conformalVoronoiMesh::buildDualFace")
+            WarningIn("Foam::conformalVoronoiMesh::buildDualFace")
                 << "Dual face uses circumcenter defined by a "
                 << "Delaunay tetrahedron with no internal "
                 << "or boundary points.  Defining Delaunay edge ends: "
                 << vA->info() << " "
                 << vB->info() << nl
-                << exit(FatalError);
+                <<endl;//<< exit(FatalError);
         }
-
+        else
+        {
             label cc1I = cc1->cellIndex();
             label cc2I = cc2->cellIndex();
 
@@ -649,6 +629,7 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
 
                 verticesOnFace.append(cc1I);
             }
+        }
 
         cc1++;
         cc2++;
@@ -838,6 +819,21 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         allGeometry_,
         foamyHexMeshDict.subDict("surfaceConformation")
     ),
+    decomposition_
+    (
+        Pstream::parRun()
+      ? new backgroundMeshDecomposition
+        (
+            runTime_,
+            rndGen_,
+            geometryToConformTo_,
+            foamyHexMeshControls().foamyHexMeshDict().subDict
+            (
+                "backgroundMeshDecomposition"
+            )
+        )
+      : NULL
+    ),
     cellShapeControl_
     (
         runTime_,
@@ -846,8 +842,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         geometryToConformTo_
     ),
     limitBounds_(),
-    featureVertices_(),
+    ftPtConformer_(*this),
-    featurePointLocations_(),
     edgeLocationTreePtr_(),
     surfacePtLocationTreePtr_(),
     surfaceConformationVertices_(),
@@ -856,7 +851,11 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         initialPointsMethod::New
         (
             foamyHexMeshDict.subDict("initialPoints"),
-            *this
+            runTime_,
+            rndGen_,
+            geometryToConformTo_,
+            cellShapeControl_,
+            decomposition_
         )
     ),
     relaxationModel_
@@ -873,8 +872,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         (
             foamyHexMeshDict.subDict("motionControl")
         )
-    ),
+    )
-    decomposition_()
 {}
 
 
@@ -893,28 +891,11 @@ void Foam::conformalVoronoiMesh::initialiseForMotion()
         geometryToConformTo_.writeFeatureObj("foamyHexMesh");
     }
 
-    if (Pstream::parRun())
-    {
-        decomposition_.reset
-        (
-            new backgroundMeshDecomposition
-            (
-                runTime_,
-                rndGen_,
-                geometryToConformTo_,
-                foamyHexMeshControls().foamyHexMeshDict().subDict
-                (
-                    "backgroundMeshDecomposition"
-                )
-            )
-        );
-    }
-
     buildCellSizeAndAlignmentMesh();
 
     insertInitialPoints();
 
-    insertFeaturePoints();
+    insertFeaturePoints(true);
 
     setVertexSizeAndAlignment();
 
@@ -953,9 +934,10 @@ void Foam::conformalVoronoiMesh::initialiseForMotion()
 
     if (foamyHexMeshControls().objOutput())
     {
-        writePoints
+        DelaunayMeshTools::writeOBJ
         (
-            "internalPoints_" + runTime_.timeName() + ".obj",
+            time().path()/"internalPoints_" + time().timeName() + ".obj",
+            *this,
             Foam::indexedVertexEnum::vtUnassigned,
             Foam::indexedVertexEnum::vtExternalFeaturePoint
         );
@@ -1159,10 +1141,12 @@ void Foam::conformalVoronoiMesh::move()
                      && pointToBeRetained[vB->index()] == true
                     )
                     {
-                        pointsToInsert.append
+                        const Foam::point pt(0.5*(dVA + dVB));
-                        (
+
-                            toPoint<Point>(0.5*(dVA + dVB))
+                        if (internalPointIsInside(pt))
-                        );
+                        {
+                            pointsToInsert.append(toPoint(pt));
+                        }
                     }
                 }
 
@@ -1203,7 +1187,8 @@ void Foam::conformalVoronoiMesh::move()
                   > foamyHexMeshControls().cosAlignmentAcceptanceAngle()
                 )
                 {
-                    scalar targetCellSize = averageCellSize(vA, vB);
+                    scalar targetCellSize =
+                        CGAL::indexedVertexOps::averageCellSize(vA, vB);
 
                     scalar targetFaceArea = sqr(targetCellSize);
 
@@ -1292,12 +1277,28 @@ void Foam::conformalVoronoiMesh::move()
                             )
                         )
                         {
-                            const Foam::point& newPt = 0.5*(dVA + dVB);
+                            const Foam::point newPt(0.5*(dVA + dVB));
 
-                            if (positionOnThisProc(newPt))
-                            {
                             // Prevent insertions spanning surfaces
-                                pointsToInsert.append(toPoint<Point>(newPt));
+                            if (internalPointIsInside(newPt))
+                            {
+                                if (Pstream::parRun())
+                                {
+                                    if
+                                    (
+                                        decomposition().positionOnThisProcessor
+                                        (
+                                            newPt
+                                        )
+                                    )
+                                    {
+                                        pointsToInsert.append(toPoint(newPt));
+                                    }
+                                }
+                                else
+                                {
+                                    pointsToInsert.append(toPoint(newPt));
+                                }
                             }
                         }
                     }
@@ -1333,10 +1334,12 @@ void Foam::conformalVoronoiMesh::move()
                              && pointToBeRetained[vB->index()] == true
                             )
                             {
-                                pointsToInsert.append
+                                const Foam::point pt(0.5*(dVA + dVB));
-                                (
+
-                                    toPoint<Point>(0.5*(dVA + dVB))
+                                if (internalPointIsInside(pt))
-                                );
+                                {
+                                    pointsToInsert.append(toPoint(pt));
+                                }
                             }
                         }
 
@@ -1431,6 +1434,9 @@ void Foam::conformalVoronoiMesh::move()
 
     Info<< "Sum displacements" << endl;
 
+    label nPointsToRetain = 0;
+    label nPointsToRemove = 0;
+
     for
     (
         Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -1451,21 +1457,33 @@ void Foam::conformalVoronoiMesh::move()
                 // 14/1/2011.
                 // Only necessary if using an exact constructions kernel
                 // (extended precision)
-
+                Foam::point pt
-                pointsToInsert.append
-                (
-                    toPoint<Point>
                 (
                     topoint(vit->point())
                   + displacementAccumulator[vit->index()]
-                    )
                 );
+
+                if (internalPointIsInside(pt))
+                {
+                    pointsToInsert.append(toPoint(pt));
+                    nPointsToRemove++;
+                }
+
+                nPointsToRetain++;
             }
         }
     }
 
     pointsToInsert.shrink();
 
+    Info<< returnReduce
+           (
+               nPointsToRetain - nPointsToRemove,
+               sumOp<label>()
+           )
+        << " internal points are outside the domain. "
+        << "They will not be inserted." << endl;
+
     // Save displacements to file.
     if (foamyHexMeshControls().objOutput() && time().outputTime())
     {
@@ -1500,73 +1518,12 @@ void Foam::conformalVoronoiMesh::move()
 
     timeCheck("Displacement calculated");
 
-    Info<< nl << "Inserting displaced tessellation" << endl;
+    Info<< nl<< "Inserting displaced tessellation" << endl;
 
-    reinsertFeaturePoints(true);
+    insertFeaturePoints(true);
 
     insertInternalPoints(pointsToInsert, true);
 
-    {
-        // Remove internal points that have been inserted outside the surface.
-        label internalPtIsOutside = 0;
-
-        autoPtr<OBJstream> str;
-
-        if (foamyHexMeshControls().objOutput() && time().outputTime())
-        {
-            str.set
-            (
-                new OBJstream(time().path()/"internalPointsOutsideDomain.obj")
-            );
-        }
-
-        DynamicList<Vertex_handle> pointsToRemove;
-
-        for
-        (
-            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-            vit != finite_vertices_end();
-            ++vit
-        )
-        {
-            if
-            (
-                (vit->internalPoint() || vit->internalBoundaryPoint())
-             //&& !vit->referred()
-            )
-            {
-                const Foam::point& pt = topoint(vit->point());
-
-                bool inside = geometryToConformTo_.inside(pt);
-
-                if
-                (
-                    !inside
-                 || !geometryToConformTo_.globalBounds().contains(pt)
-                )
-                {
-                    if
-                    (
-                        foamyHexMeshControls().objOutput()
-                     && time().outputTime()
-                    )
-                    {
-                        str().write(topoint(vit->point()));
-                    }
-
-                    pointsToRemove.append(vit);
-                    internalPtIsOutside++;
-                }
-            }
-        }
-
-        remove(pointsToRemove.begin(), pointsToRemove.end());
-
-        Info<< "    " << returnReduce(internalPtIsOutside, sumOp<label>())
-            << " internal points were inserted outside the domain. "
-            << "They have been removed." << endl;
-    }
-
     // Fix points that have not been significantly displaced
 //    for
 //    (
@@ -1614,16 +1571,21 @@ void Foam::conformalVoronoiMesh::move()
 
     if (foamyHexMeshControls().objOutput())
     {
-        writePoints
+        DelaunayMeshTools::writeOBJ
         (
-            "internalPoints_" + time().timeName() + ".obj",
+            time().path()/"internalPoints_" + time().timeName() + ".obj",
+            *this,
             Foam::indexedVertexEnum::vtInternal
         );
-    }
 
-    if (foamyHexMeshControls().objOutput() && time().outputTime())
+        if (reconformToSurface())
         {
-        writeBoundaryPoints("boundaryPoints_" + time().timeName() + ".obj");
+            DelaunayMeshTools::writeBoundaryPoints
+            (
+                time().path()/"boundaryPoints_" + time().timeName() + ".obj",
+                *this
+            );
+        }
     }
 
     timeCheck("After conformToSurface");
@@ -1646,64 +1608,6 @@ void Foam::conformalVoronoiMesh::move()
 }
 
 
-bool Foam::conformalVoronoiMesh::positionOnThisProc
-(
-    const Foam::point& pt
-) const
-{
-    if (Pstream::parRun())
-    {
-        return decomposition_().positionOnThisProcessor(pt);
-    }
-
-    return true;
-}
-
-
-Foam::boolList Foam::conformalVoronoiMesh::positionOnThisProc
-(
-    const Foam::List<Foam::point>& pts
-) const
-{
-    if (Pstream::parRun())
-    {
-        return decomposition_().positionOnThisProcessor(pts);
-    }
-
-    return boolList(pts.size(), true);
-}
-
-
-Foam::labelList Foam::conformalVoronoiMesh::positionProc
-(
-    const Foam::List<Foam::point>& pts
-) const
-{
-    if (!Pstream::parRun())
-    {
-        return labelList(pts.size(), -1);
-    }
-
-    return decomposition_().processorPosition(pts);
-}
-
-
-Foam::List<Foam::List<Foam::pointIndexHit> >
-Foam::conformalVoronoiMesh::intersectsProc
-(
-    const List<Foam::point>& starts,
-    const List<Foam::point>& ends
-) const
-{
-    if (!Pstream::parRun())
-    {
-        return List<List<pointIndexHit> >(starts.size());
-    }
-
-    return decomposition_().intersectsProcessors(starts, ends, false);
-}
-
-
 //Foam::labelListList Foam::conformalVoronoiMesh::overlapsProc
 //(
 //    const List<Foam::point>& centres,
@@ -1776,7 +1680,7 @@ void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
                 << "    quality = " << quality << nl
                 << "    dual    = " << topoint(cit->dual()) << endl;
 
-            drawDelaunayCell(str, cit, badCells++);
+            DelaunayMeshTools::drawDelaunayCell(str, cit, badCells++);
 
             FixedList<PointExact, 4> cellVerticesExact(PointExact(0,0,0));
             forAll(cellVerticesExact, vI)

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -74,9 +74,9 @@ SourceFiles
 #include "globalIndex.H"
 #include "pointFeatureEdgesTypes.H"
 #include "pointConversion.H"
-#include "Tuple2.H"
+#include "Pair.H"
 #include "DistributedDelaunayMesh.H"
-#include "tensorIOField.H"
+#include "featurePointConformer.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -103,12 +103,6 @@ class conformalVoronoiMesh
 {
 public:
 
-    enum reconformationMode
-    {
-        rmOff,   // Do not rebuild the surface conformation
-        rmOn
-    };
-
     typedef Delaunay::Vertex_handle    Vertex_handle;
     typedef Delaunay::Cell_handle      Cell_handle;
     typedef Delaunay::Edge             Edge;
@@ -126,8 +120,6 @@ private:
 
     // Static data
 
-        static const scalar tolParallel;
-
         static const scalar searchConeAngle;
 
         static const scalar searchAngleOppositeSurface;
@@ -151,19 +143,17 @@ private:
         //- The surfaces to conform to
         conformationSurfaces geometryToConformTo_;
 
+        //- Background mesh decomposition, only available in parallel.
+        autoPtr<backgroundMeshDecomposition> decomposition_;
+
         //- The cell shape control object
         cellShapeControl cellShapeControl_;
 
         //- Limiting bound box before infinity begins
         treeBoundBox limitBounds_;
 
-        //- Store the feature constraining points to be reinserted after a
+        //-
-        //  triangulation clear. Maintained with relative types and indices.
+        featurePointConformer ftPtConformer_;
-        List<Vb> featureVertices_;
-
-        //- Storing the locations of all of the features to be conformed to.
-        //  Single pointField required by the featurePointTree.
-        pointField featurePointLocations_;
 
         //- Search tree for edge point locations
         mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
@@ -190,9 +180,6 @@ private:
         //-  Face area weight function.  Runtime selectable.
         autoPtr<faceAreaWeightModel> faceAreaWeightModel_;
 
-        //- Background mesh decomposition, only available in parallel.
-        autoPtr<backgroundMeshDecomposition> decomposition_;
-
 
     // Private Member Functions
 
@@ -203,14 +190,6 @@ private:
         //  to be on a surface.
         inline scalar targetCellSize(const Foam::point& pt) const;
 
-        //- Return the target cell size from that stored on a pair of
-        //  Delaunay vertices, using a mean function.
-        inline scalar averageCellSize
-        (
-            const Vertex_handle& vA,
-            const Vertex_handle& vB
-        ) const;
-
         //- Return the target cell size from that stored on a pair of
         //  Delaunay vertices, including the possibility that one of
         //  them is not an internalOrBoundaryPoint, and so will not
@@ -228,36 +207,6 @@ private:
             const Delaunay::Finite_facets_iterator& fit
         ) const;
 
-        //- Return the local point pair separation at the given location
-        inline scalar pointPairDistance(const Foam::point& pt) const;
-
-        //- Return the local mixed feature point placement distance
-        inline scalar mixedFeaturePointDistance(const Foam::point& pt) const;
-
-        //- Return the square of the local feature point exclusion distance
-        inline scalar featurePointExclusionDistanceSqr
-        (
-            const Foam::point& pt
-        ) const;
-
-        //- Return the square of the local feature edge exclusion distance
-        inline scalar featureEdgeExclusionDistanceSqr
-        (
-            const Foam::point& pt
-        ) const;
-
-        //- Return the square of the local surface point exclusion distance
-        inline scalar surfacePtExclusionDistanceSqr
-        (
-            const Foam::point& pt
-        ) const;
-
-        //- Return the square of the local surface search distance
-        inline scalar surfaceSearchDistanceSqr(const Foam::point& pt) const;
-
-        //- Return the local maximum surface protrusion distance
-        inline scalar maxSurfaceProtrusion(const Foam::point& pt) const;
-
         //- Insert Delaunay vertices using the CGAL range insertion method,
         //  optionally check processor occupancy and distribute to other
         //  processors
@@ -281,7 +230,7 @@ private:
             const Foam::point& surfPt,
             const vector& n,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         inline Foam::point perturbPoint(const Foam::point& pt) const;
 
@@ -293,7 +242,10 @@ private:
             const Foam::point& surfPt,
             const vector& n,
             DynamicList<Vb>& pts
-        );
+        ) const;
+
+        //- Check internal point is completely inside the meshable region
+        inline bool internalPointIsInside(const Foam::point& pt) const;
 
         inline bool isPointPair
         (
@@ -318,20 +270,12 @@ private:
             const fileName fName = fileName::null
         );
 
-        //- Call the appropriate function to conform to an edge
-        void createEdgePointGroup
-        (
-            const extendedFeatureEdgeMesh& feMesh,
-            const pointIndexHit& edHit,
-            DynamicList<Vb>& pts
-        );
-
         void createEdgePointGroupByCirculating
         (
             const extendedFeatureEdgeMesh& feMesh,
             const pointIndexHit& edHit,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         bool meshableRegion
         (
@@ -354,7 +298,7 @@ private:
             const extendedFeatureEdgeMesh& feMesh,
             const pointIndexHit& edHit,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         //- Create points to conform to an internal edge
         void createInternalEdgePointGroup
@@ -362,7 +306,7 @@ private:
             const extendedFeatureEdgeMesh& feMesh,
             const pointIndexHit& edHit,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         //- Create points to conform to a flat edge
         void createFlatEdgePointGroup
@@ -370,7 +314,7 @@ private:
             const extendedFeatureEdgeMesh& feMesh,
             const pointIndexHit& edHit,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         //- Create points to conform to an open edge
         void createOpenEdgePointGroup
@@ -378,7 +322,7 @@ private:
             const extendedFeatureEdgeMesh& feMesh,
             const pointIndexHit& edHit,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         //- Create points to conform to multiply connected edge
         void createMultipleEdgePointGroup
@@ -386,85 +330,10 @@ private:
             const extendedFeatureEdgeMesh& feMesh,
             const pointIndexHit& edHit,
             DynamicList<Vb>& pts
-        );
+        ) const;
 
         //- Determine and insert point groups at the feature points
-        void insertFeaturePoints();
+        void insertFeaturePoints(bool distribute = false);
-
-        //- Create point groups at mixed feature points
-        void createMixedFeaturePoints(DynamicList<Vb>& pts);
-
-        void addMasterAndSlavePoints
-        (
-            const DynamicList<Foam::point>& masterPoints,
-            const DynamicList<indexedVertexEnum::vertexType>&masterPointsTypes,
-            const Map<DynamicList<autoPtr<plane> > >& masterPointReflections,
-            DynamicList<Vb>& pts,
-            const label ptI
-        ) const;
-
-        label getSign(const extendedFeatureEdgeMesh::edgeStatus eStatus) const;
-
-        void createMasterAndSlavePoints
-        (
-            const extendedFeatureEdgeMesh& feMesh,
-            const label ptI,
-            DynamicList<Vb>& pts
-        ) const;
-
-        void createFeaturePoints(DynamicList<Vb>& pts);
-
-        vector sharedFaceNormal
-        (
-            const extendedFeatureEdgeMesh& feMesh,
-            const label edgeI,
-            const label nextEdgeI
-        ) const;
-
-        List<Foam::point> reflectPointInPlanes
-        (
-            const Foam::point p,
-            const DynamicList<autoPtr<plane> >& planes
-        ) const;
-
-        Foam::point reflectPointInPlane
-        (
-            const Foam::point p,
-            const plane& planeN
-        ) const;
-
-
-        //- Fill the pointFeatureEdgesType struct with the types of feature
-        //  edges that are attached to the point.
-        List<extendedFeatureEdgeMesh::edgeStatus> calcPointFeatureEdgesTypes
-        (
-            const extendedFeatureEdgeMesh& feMesh,
-            const labelList& pEds,
-            pointFeatureEdgesTypes& pFEdgeTypes
-        ) const;
-
-        //- Create feature point groups if a specialisation exists for the
-        //  structure
-        bool createSpecialisedFeaturePoint
-        (
-            const extendedFeatureEdgeMesh& feMesh,
-            const labelList& pEds,
-            const pointFeatureEdgesTypes& pFEdgeTypes,
-            const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
-            const label ptI,
-            DynamicList<Vb>& pts
-        );
-
-        //- Store the locations of all of the features to be conformed to
-        void constructFeaturePointLocations();
-
-        List<pointIndexHit> findSurfacePtLocationsNearFeaturePoint
-        (
-            const Foam::point& featurePoint
-        ) const;
-
-        //- Reinsert stored feature point defining points
-        void reinsertFeaturePoints(bool distribute = false);
 
         //- Check if a location is in exclusion range around a feature point
         bool nearFeaturePt(const Foam::point& pt) const;
@@ -473,10 +342,6 @@ private:
         //  initialPointsMethod
         void insertInitialPoints();
 
-        //- Calculate the worst load balance
-        template<class Triangulation>
-        scalar calculateLoadUnbalance(const Triangulation& mesh) const;
-
         //- In parallel redistribute the backgroundMeshDecomposition and
         //  vertices to balance the number of vertices on each processor.
         //  Returns true if the background mesh changes as this removes all
@@ -534,7 +399,7 @@ private:
 
         //- Decision making function for when to rebuild the surface
         //  conformation
-        reconformationMode reconformationControl() const;
+        bool reconformToSurface() const;
 
         //- Determines geometrically whether a vertex is close to a surface
         //  This is an optimisation
@@ -850,9 +715,6 @@ private:
         //- Create the cell centres to use for the mesh
         void createCellCentres(pointField& cellCentres) const;
 
-        //- Extract all points in vertex-index order
-        tmp<pointField> allPoints() const;
-
         //- Sort the faces, owner and neighbour lists into
         //  upper-triangular order.  For internal faces only, use
         //  before adding patch faces
@@ -974,22 +836,6 @@ public:
         //  surface as required
         void move();
 
-        //- Check if the point is in the domain handled by this processor
-        bool positionOnThisProc(const Foam::point& pt) const;
-
-        //- Check if the point is in the domain handled by this processor
-        boolList positionOnThisProc(const Foam::List<Foam::point>& pts) const;
-
-        //- Which processor's domain handles this point
-        labelList positionProc(const Foam::List<Foam::point>& pts) const;
-
-        //- Which other processors does each line segment intersect
-        List<List<pointIndexHit> > intersectsProc
-        (
-            const List<Foam::point>& starts,
-            const List<Foam::point>& ends
-        ) const;
-
 //        //- Which other processors does each sphere overlap
 //        labelListList overlapsProc
 //        (
@@ -1022,6 +868,50 @@ public:
             inline const cvControls& foamyHexMeshControls() const;
 
 
+        // Query
+
+            //- Return the local point pair separation at the given location
+            inline scalar pointPairDistance(const Foam::point& pt) const;
+
+            //- Return the local mixed feature point placement distance
+            inline scalar mixedFeaturePointDistance
+            (
+                const Foam::point& pt
+            ) const;
+
+            //- Return the square of the local feature point exclusion distance
+            inline scalar featurePointExclusionDistanceSqr
+            (
+                const Foam::point& pt
+            ) const;
+
+            //- Return the square of the local feature edge exclusion distance
+            inline scalar featureEdgeExclusionDistanceSqr
+            (
+                const Foam::point& pt
+            ) const;
+
+            //- Return the square of the local surface point exclusion distance
+            inline scalar surfacePtExclusionDistanceSqr
+            (
+                const Foam::point& pt
+            ) const;
+
+            //- Return the square of the local surface search distance
+            inline scalar surfaceSearchDistanceSqr(const Foam::point& pt) const;
+
+            //- Return the local maximum surface protrusion distance
+            inline scalar maxSurfaceProtrusion(const Foam::point& pt) const;
+
+            //- Call the appropriate function to conform to an edge
+            void createEdgePointGroup
+            (
+                const extendedFeatureEdgeMesh& feMesh,
+                const pointIndexHit& edHit,
+                DynamicList<Vb>& pts
+            ) const;
+
+
         // Write
 
             //- Write the elapsedCpuTime and memory usage, with an optional
@@ -1038,53 +928,6 @@ public:
                 const string& description = string::null
             ) const;
 
-            //- Write the Delaunay cell
-            void drawDelaunayCell
-            (
-                Ostream& os,
-                const Cell_handle& c,
-                label offset = 0
-            ) const;
-
-            //- Write Delaunay points in the range between (and including)
-            //  type startPointType and endPointType to .obj file
-            void writePoints
-            (
-                const fileName& fName,
-                const Foam::indexedVertexEnum::vertexType startPointType,
-                const Foam::indexedVertexEnum::vertexType endPointType
-            ) const;
-
-            //- Write Delaunay points of type pointType to .obj file
-            void writePoints
-            (
-                const fileName& fName,
-                const Foam::indexedVertexEnum::vertexType pointType
-            ) const;
-
-            void writeFixedPoints(const fileName& fName) const;
-
-            //- Write the boundary Delaunay points to .obj file
-            void writeBoundaryPoints(const fileName& fName) const;
-
-            //- Write list of points to file
-            void writePoints
-            (
-                const fileName& fName,
-                const List<Foam::point>& points
-            ) const;
-
-            //- Write list of points to file
-            void writePoints
-            (
-                const fileName& fName,
-                const List<Vb>& points
-            ) const;
-
-            //- Write the internal Delaunay vertices of the tessellation as a
-            //  pointField that may be used to restart the meshing process
-            void writeInternalDelaunayVertices(const fileName& instance) const;
-
             //- Prepare data and call writeMesh for polyMesh and
             //  tetDualMesh
             void writeMesh(const fileName& instance);
@@ -1105,14 +948,6 @@ public:
                 const PackedBoolList& boundaryFacesToRemove
             ) const;
 
-            //- Write points and faces as .obj file
-            void writeObjMesh
-            (
-                const pointField& points,
-                const faceList& faces,
-                const fileName& fName
-            ) const;
-
             //- Calculate and write a field of the target cell size,
             //  target cell volume, actual cell volume and equivalent
             //  actual cell size (cbrt(actual cell volume)).
@@ -1126,12 +961,6 @@ public:
             //- Find the cellSet of the boundary cells which have points that
             //  protrude out of the surface beyond a tolerance.
             labelHashSet findRemainingProtrusionSet(const polyMesh& mesh) const;
-
-            void writeProcessorInterface
-            (
-                const fileName& fName,
-                const faceList& faces
-            ) const;
 };
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -29,6 +29,8 @@ License
 #include "polyMeshGeometry.H"
 #include "indexedCellChecks.H"
 #include "OBJstream.H"
+#include "indexedCellOps.H"
+#include "DelaunayMeshTools.H"
 
 #include "CGAL/Exact_predicates_exact_constructions_kernel.h"
 #include "CGAL/Gmpq.h"
@@ -598,7 +600,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 
     // deferredCollapseFaceSet(owner, neighbour, deferredCollapseFaces);
 
-    cellCentres = allPoints();
+    cellCentres = DelaunayMeshTools::allPoints(*this);
 
     cellToDelaunayVertex = removeUnusedCells(owner, neighbour);
 
@@ -1227,7 +1229,7 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
     );
 
     //createCellCentres(cellCentres);
-    cellCentres = allPoints();
+    cellCentres = DelaunayMeshTools::allPoints(*this);
 
     labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
     cellCentres = pointField(cellCentres, cellToDelaunayVertex);
@@ -2443,7 +2445,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                          && vc2->neighbor(cI)->hasConstrainedPoint()
                         )
                         {
-                            drawDelaunayCell
+                            DelaunayMeshTools::drawDelaunayCell
                             (
                                 cellStr,
                                 vc2->neighbor(cI),
@@ -2837,7 +2839,12 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                       + name(neighbour)
                       + "_interface.obj";
 
-                    writeProcessorInterface(fName, procPatchFaces);
+                    DelaunayMeshTools::writeProcessorInterface
+                    (
+                        time().path()/fName,
+                        *this,
+                        procPatchFaces
+                    );
                 }
             }
         }
@@ -2871,30 +2878,6 @@ void Foam::conformalVoronoiMesh::createCellCentres
 }
 
 
-Foam::tmp<Foam::pointField> Foam::conformalVoronoiMesh::allPoints() const
-{
-    tmp<pointField> tpts(new pointField(number_of_vertices(), point::max));
-    pointField& pts = tpts();
-
-    label nVert = 0;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->internalOrBoundaryPoint())
-        {
-            pts[nVert++] = topoint(vit->point());
-        }
-    }
-
-    return tpts;
-}
-
-
 void Foam::conformalVoronoiMesh::sortFaces
 (
     faceList& faces,

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -55,7 +55,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
         cit->cellIndex() = Cb::ctUnassigned;
     }
 
-    if (reconformationControl() == rmOff)
+    if (!reconformToSurface())
     {
         // Reinsert stored surface conformation
         reinsertSurfaceConformation();
@@ -87,8 +87,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 }
 
 
-Foam::conformalVoronoiMesh::reconformationMode
+bool Foam::conformalVoronoiMesh::reconformToSurface() const
-Foam::conformalVoronoiMesh::reconformationControl() const
 {
     if
     (
@@ -96,10 +95,10 @@ Foam::conformalVoronoiMesh::reconformationControl() const
       % foamyHexMeshControls().surfaceConformationRebuildFrequency() == 0
     )
     {
-        return rmOn;
+        return true;
     }
 
-    return rmOff;
+    return false;
 }
 
 
@@ -220,20 +219,10 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
 {
     timeCheck("Start buildSurfaceConformation");
 
-    if (reconformationControl() == rmOff)
+    Info<< nl
-    {
+        << "Rebuilding surface conformation for more iterations"
-        WarningIn("buildSurfaceConformation()")
-            << "reconformationMode rmNone specified, not building conformation"
         << endl;
 
-        return;
-    }
-    else
-    {
-        Info<< nl << "Rebuilding surface conformation for more iterations"
-            << endl;
-    }
-
     existingEdgeLocations_.clearStorage();
     existingSurfacePtLocations_.clearStorage();
 
@@ -1297,14 +1286,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
         if
         (
             is_infinite(c1) || is_infinite(c2)
-         || !c1->hasInternalPoint() || !c2->hasInternalPoint()
+         || (!c1->hasInternalPoint() && !c2->hasInternalPoint())
          || !c1->real() || !c2->real()
         )
         {
             continue;
         }
 
-        Foam::point edgeMid = 0.5*(c1->dual() + c2->dual());
+//        Foam::point endPt = 0.5*(c1->dual() + c2->dual());
+        Foam::point endPt = c1->dual();
+
+        if
+        (
+            magSqr(vert - c1->dual())
+          < magSqr(vert - c2->dual())
+        )
+        {
+            endPt = c2->dual();
+        }
 
         pointIndexHit surfHit;
         label hitSurface;
@@ -1312,7 +1311,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
         geometryToConformTo_.findSurfaceAnyIntersection
         (
             vert,
-            edgeMid,
+            endPt,
             surfHit,
             hitSurface
         );
@@ -1330,7 +1329,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
             const vector& n = norm[0];
 
             const scalar normalProtrusionDistance =
-                (edgeMid - surfHit.hitPoint()) & n;
+                (endPt - surfHit.hitPoint()) & n;
 
             if (normalProtrusionDistance > maxProtrusionDistance)
             {
@@ -1347,7 +1346,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
     if
     (
         surfHitLargest.hit()
-     && !positionOnThisProc(surfHitLargest.hitPoint())
+     && (
+            Pstream::parRun()
+         && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
+        )
     )
     {
         // A protrusion was identified, but not penetrating on this processor,
@@ -1390,14 +1392,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
         if
         (
             is_infinite(c1) || is_infinite(c2)
-         || !c1->hasInternalPoint() || !c2->hasInternalPoint()
+         || (!c1->hasInternalPoint() && !c2->hasInternalPoint())
          || !c1->real() || !c2->real()
         )
         {
             continue;
         }
 
-        Foam::point edgeMid = 0.5*(c1->dual() + c2->dual());
+//        Foam::point endPt = 0.5*(c1->dual() + c2->dual());
+        Foam::point endPt = c1->dual();
+
+        if
+        (
+            magSqr(vert - c1->dual())
+          < magSqr(vert - c2->dual())
+        )
+        {
+            endPt = c2->dual();
+        }
 
         pointIndexHit surfHit;
         label hitSurface;
@@ -1405,7 +1417,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
         geometryToConformTo_.findSurfaceAnyIntersection
         (
             vert,
-            edgeMid,
+            endPt,
             surfHit,
             hitSurface
         );
@@ -1422,8 +1434,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
 
             const vector& n = norm[0];
 
-            scalar normalIncursionDistance =
+            scalar normalIncursionDistance = (endPt - surfHit.hitPoint()) & n;
-                (edgeMid - surfHit.hitPoint()) & n;
 
             if (normalIncursionDistance < minIncursionDistance)
             {
@@ -1445,7 +1456,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
     if
     (
         surfHitLargest.hit()
-     && !positionOnThisProc(surfHitLargest.hitPoint())
+     && (
+            Pstream::parRun()
+         && !decomposition().positionOnThisProcessor(surfHitLargest.hitPoint())
+        )
     )
     {
         // A protrusion was identified, but not penetrating on this processor,
@@ -1468,7 +1482,11 @@ void Foam::conformalVoronoiMesh::reportProcessorOccupancy()
     {
         if (vit->real())
         {
-            if (!positionOnThisProc(topoint(vit->point())))
+            if
+            (
+                Pstream::parRun()
+             && !decomposition().positionOnThisProcessor(topoint(vit->point()))
+            )
             {
                 Pout<< topoint(vit->point()) << " is not on this processor "
                     << endl;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C

@@ -28,6 +28,7 @@ License
 #include "triangle.H"
 #include "tetrahedron.H"
 #include "const_circulator.H"
+#include "DelaunayMeshTools.H"
 
 using namespace Foam::vectorTools;
 
@@ -38,7 +39,7 @@ void Foam::conformalVoronoiMesh::createEdgePointGroup
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+) const
 {
     if (foamyHexMeshControls().circulateEdges())
     {
@@ -164,7 +165,7 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+)  const
 {
     typedef Foam::indexedVertexEnum::vertexType         vertexType;
     typedef extendedFeatureEdgeMesh::sideVolumeType     sideVolumeType;
@@ -413,6 +414,7 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
 
                 // Remove the first reflection plane if we are no longer
                 // circulating
+
                 masterPointReflectionsPrev.erase(initialRegion);
                 masterPointReflectionsNext.erase(circ.nRotations());
             }
@@ -480,7 +482,7 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+) const
 {
     const Foam::point& edgePt = edHit.hitPoint();
 
@@ -581,7 +583,7 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+) const
 {
     const Foam::point& edgePt = edHit.hitPoint();
 
@@ -757,7 +759,7 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+) const
 {
     const Foam::point& edgePt = edHit.hitPoint();
 
@@ -789,7 +791,7 @@ void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+) const
 {
     // Assume it is a baffle and insert flat edge point pairs
     const Foam::point& edgePt = edHit.hitPoint();
@@ -836,7 +838,7 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
     const extendedFeatureEdgeMesh& feMesh,
     const pointIndexHit& edHit,
     DynamicList<Vb>& pts
-)
+) const
 {
 //    Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
 
@@ -867,618 +869,27 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
 }
 
 
-void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
+void Foam::conformalVoronoiMesh::insertFeaturePoints(bool distribute)
 {
-    Info<< nl << "Reinserting stored feature points" << endl;
+    Info<< nl
+        << "Inserting feature points" << endl;
 
-    insertPoints(featureVertices_, distribute);
+    const label preFeaturePointSize(number_of_vertices());
-}
 
-
+    if (Pstream::parRun() && distribute)
-void Foam::conformalVoronoiMesh::createMixedFeaturePoints
-(
-    DynamicList<Vb>& pts
-)
-{
-    if (foamyHexMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
     {
-        Info<< nl << "Skipping specialised handling for mixed feature points"
+        ftPtConformer_.distribute(decomposition());
-            << endl;
-        return;
     }
 
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
-    (
-        geometryToConformTo_.features()
-    );
 
-    forAll(feMeshes, i)
+    // Insert the created points directly as already distributed.
-    {
+    this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices);
-        const extendedFeatureEdgeMesh& feMesh = feMeshes[i];
-        const labelListList& pointsEdges = feMesh.pointEdges();
-        const pointField& points = feMesh.points();
-
-        for
-        (
-            label ptI = feMesh.mixedStart();
-            ptI < feMesh.nonFeatureStart();
-            ptI++
-        )
-        {
-            const Foam::point& featPt = points[ptI];
-
-            if (!positionOnThisProc(featPt))
-            {
-                continue;
-            }
-
-            const labelList& pEds = pointsEdges[ptI];
-
-            pointFeatureEdgesTypes pFEdgeTypes(ptI);
-
-            const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat
-                = calcPointFeatureEdgesTypes(feMesh, pEds, pFEdgeTypes);
-
-            bool specialisedSuccess = false;
-
-            if (foamyHexMeshControls().specialiseFeaturePoints())
-            {
-                specialisedSuccess = createSpecialisedFeaturePoint
-                (
-                    feMesh, pEds, pFEdgeTypes, allEdStat, ptI, pts
-                );
-            }
-
-            if (!specialisedSuccess && foamyHexMeshControls().edgeAiming())
-            {
-                // Specialisations available for some mixed feature points.  For
-                // non-specialised feature points, inserting mixed internal and
-                // external edge groups at feature point.
-
-                // Skipping unsupported mixed feature point types
-//                bool skipEdge = false;
-//
-//                forAll(pEds, e)
-//                {
-//                    const label edgeI = pEds[e];
-//
-//                    const extendedFeatureEdgeMesh::edgeStatus edStatus
-//                        = feMesh.getEdgeStatus(edgeI);
-//
-//                    if
-//                    (
-//                        edStatus == extendedFeatureEdgeMesh::OPEN
-//                     || edStatus == extendedFeatureEdgeMesh::MULTIPLE
-//                    )
-//                    {
-//                        Info<< "Edge type " << edStatus
-//                            << " found for mixed feature point " << ptI
-//                            << ". Not supported."
-//                            << endl;
-//
-//                        skipEdge = true;
-//                    }
-//                }
-//
-//                if (skipEdge)
-//                {
-//                    Info<< "Skipping point " << ptI << nl << endl;
-//
-//                    continue;
-//                }
-
-//                createFeaturePoints(feMesh, ptI, pts, types);
-
-                const Foam::point pt = points[ptI];
-
-                const scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
-
-                forAll(pEds, e)
-                {
-                    const label edgeI = pEds[e];
-
-                    const Foam::point edgePt =
-                        pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
-
-                    const pointIndexHit edgeHit(true, edgePt, edgeI);
-
-                    createEdgePointGroup(feMesh, edgeHit, pts);
-                }
-            }
-        }
-    }
-}
-//
-//
-//void Foam::conformalVoronoiMesh::createFeaturePoints
-//(
-//    DynamicList<Foam::point>& pts,
-//    DynamicList<label>& indices,
-//    DynamicList<label>& types
-//)
-//{
-//    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-//    (
-//        geometryToConformTo_.features()
-//    );
-//
-//    forAll(feMeshes, i)
-//    {
-//        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-//
-//        for
-//        (
-//            label ptI = feMesh.convexStart();
-//            ptI < feMesh.mixedStart();
-//            ++ptI
-//        )
-//        {
-//            const Foam::point& featPt = feMesh.points()[ptI];
-//
-//            if (!positionOnThisProc(featPt))
-//            {
-//                continue;
-//            }
-//
-//            const scalar searchRadiusSqr = 5.0*targetCellSize(featPt);
-//
-//            labelList indices = surfacePtLocationTreePtr_().findSphere
-//            (
-//                featPt,
-//                searchRadiusSqr
-//            );
-//
-//            pointField nearestSurfacePoints(indices.size());
-//
-//            forAll(indices, pI)
-//            {
-//                nearestSurfacePoints[pI] =
-//                    surfaceConformationVertices_[indices[pI]];
-//            }
-//
-//            forAll(feMesh.)
-//
-//            // Now find the nearest points within the edge cones.
-//
-//            // Calculate preliminary surface point locations
-//
-//
-//        }
-//    }
-//}
-
-
-void Foam::conformalVoronoiMesh::insertFeaturePoints()
-{
-    Info<< nl << "Conforming to feature points" << endl;
-
-    Info<< "    Circulating edges is: "
-        << foamyHexMeshControls().circulateEdges().asText() << nl
-        << "    Guarding feature points is: "
-        << foamyHexMeshControls().guardFeaturePoints().asText() << nl
-        << "    Snapping to feature points is: "
-        << foamyHexMeshControls().snapFeaturePoints().asText() << nl
-        << "    Specialising feature points is: "
-        << foamyHexMeshControls().specialiseFeaturePoints().asText() << endl;
-
-
-    DynamicList<Vb> pts;
-
-    const label preFeaturePointSize = number_of_vertices();
-
-    createFeaturePoints(pts);
-
-    createMixedFeaturePoints(pts);
-
-    // Points added using the createEdgePointGroup function will be labelled as
-    // internal/external feature edges. Relabel them as feature points,
-    // otherwise they are inserted as both feature points and surface points.
-    forAll(pts, pI)
-    {
-        Vb& pt = pts[pI];
-
-        //if (pt.featureEdgePoint())
-        {
-            if (pt.internalBoundaryPoint())
-            {
-                pt.type() = Vb::vtInternalFeaturePoint;
-            }
-            else if (pt.externalBoundaryPoint())
-            {
-                pt.type() = Vb::vtExternalFeaturePoint;
-            }
-        }
-    }
-
-    // Insert the created points, distributing to the appropriate processor
-    insertPoints(pts, true);
-
-    if (foamyHexMeshControls().objOutput())
-    {
-        writePoints("featureVertices.obj", pts);
-    }
 
     label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
     reduce(nFeatureVertices, sumOp<label>());
 
     Info<< "    Inserted " << nFeatureVertices << " feature vertices" << endl;
-
-    featureVertices_.clear();
-    featureVertices_.setSize(pts.size());
-
-    forAll(pts, pI)
-    {
-        featureVertices_[pI] = pts[pI];
-    }
-
-    constructFeaturePointLocations();
-}
-
-
-void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
-{
-    DynamicList<Foam::point> ftPtLocs;
-
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-    (
-        geometryToConformTo_.features()
-    );
-
-    forAll(feMeshes, i)
-    {
-        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-
-        if (foamyHexMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
-        {
-            // Ignoring mixed feature points
-            for
-            (
-                label ptI = feMesh.convexStart();
-                ptI < feMesh.mixedStart();
-                ptI++
-            )
-            {
-                ftPtLocs.append(feMesh.points()[ptI]);
-            }
-        }
-        else
-        {
-            for
-            (
-                label ptI = feMesh.convexStart();
-                ptI < feMesh.nonFeatureStart();
-                ptI++
-            )
-            {
-                ftPtLocs.append(feMesh.points()[ptI]);
-            }
-        }
-    }
-
-    featurePointLocations_.transfer(ftPtLocs);
-}
-
-
-Foam::List<Foam::pointIndexHit>
-Foam::conformalVoronoiMesh::findSurfacePtLocationsNearFeaturePoint
-(
-    const Foam::point& featurePoint
-) const
-{
-    DynamicList<pointIndexHit> dynPointList;
-
-    const scalar searchRadiusSqr = 3*targetCellSize(featurePoint);
-
-    labelList surfacePtList = surfacePtLocationTreePtr_().findSphere
-    (
-       featurePoint,
-       searchRadiusSqr
-    );
-
-    forAll(surfacePtList, elemI)
-    {
-       label index = surfacePtList[elemI];
-
-       const Foam::point& p
-           = surfacePtLocationTreePtr_().shapes().shapePoints()[index];
-
-       pointIndexHit nearHit(true, p, index);
-
-       dynPointList.append(nearHit);
-    }
-
-    return dynPointList.shrink();
-}
-
-
-void Foam::conformalVoronoiMesh::addMasterAndSlavePoints
-(
-    const DynamicList<Foam::point>& masterPoints,
-    const DynamicList<Foam::indexedVertexEnum::vertexType>& masterPointsTypes,
-    const Map<DynamicList<autoPtr<plane> > >& masterPointReflections,
-    DynamicList<Vb>& pts,
-    const label ptI
-) const
-{
-    typedef DynamicList<autoPtr<plane> >        planeDynList;
-    typedef Foam::indexedVertexEnum::vertexType vertexType;
-
-    forAll(masterPoints, pI)
-    {
-        // Append master to the list of points
-
-        const Foam::point& masterPt = masterPoints[pI];
-        const vertexType masterType = masterPointsTypes[pI];
-
-//        Info<< "    Master = " << masterPt << endl;
-
-        pts.append
-        (
-            Vb
-            (
-                masterPt,
-                vertexCount() + pts.size(),
-                masterType,
-                Pstream::myProcNo()
-            )
-        );
-
-//        const label masterIndex = pts[pts.size() - 1].index();
-
-        //meshTools::writeOBJ(strMasters, masterPt);
-
-        const planeDynList& masterPointPlanes = masterPointReflections[pI];
-
-        forAll(masterPointPlanes, planeI)
-        {
-            // Reflect master points in the planes and insert the slave points
-
-            const plane& reflPlane = masterPointPlanes[planeI]();
-
-            const Foam::point slavePt = reflPlane.mirror(masterPt);
-
-//            Info<< "        Slave " << planeI << " = " << slavePt << endl;
-
-            const vertexType slaveType =
-            (
-                masterType == Vb::vtInternalFeaturePoint
-              ? Vb::vtExternalFeaturePoint // true
-              : Vb::vtInternalFeaturePoint // false
-            );
-
-            pts.append
-            (
-                Vb
-                (
-                    slavePt,
-                    vertexCount() + pts.size(),
-                    slaveType,
-                    Pstream::myProcNo()
-                )
-            );
-
-
-            //meshTools::writeOBJ(strSlaves, slavePt);
-        }
-    }
-}
-
-
-Foam::label Foam::conformalVoronoiMesh::getSign
-(
-    const extendedFeatureEdgeMesh::edgeStatus eStatus
-) const
-{
-    if (eStatus == extendedFeatureEdgeMesh::EXTERNAL)
-    {
-        return -1;
-    }
-    else if (eStatus == extendedFeatureEdgeMesh::INTERNAL)
-    {
-        return 1;
-    }
-
-    return 0;
-}
-
-
-void Foam::conformalVoronoiMesh::createMasterAndSlavePoints
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const label ptI,
-    DynamicList<Vb>& pts
-) const
-{
-    typedef DynamicList<autoPtr<plane> >                planeDynList;
-    typedef Foam::indexedVertexEnum::vertexType         vertexType;
-    typedef Foam::extendedFeatureEdgeMesh::edgeStatus   edgeStatus;
-
-    const Foam::point& featPt = feMesh.points()[ptI];
-
-    if (!positionOnThisProc(featPt) || geometryToConformTo_.outside(featPt))
-    {
-        return;
-    }
-
-    const scalar ppDist = pointPairDistance(featPt);
-
-    // Maintain a list of master points and the planes to relect them in
-    DynamicList<Foam::point> masterPoints;
-    DynamicList<vertexType> masterPointsTypes;
-    Map<planeDynList> masterPointReflections;
-
-    const labelList& featPtEdges = feMesh.featurePointEdges()[ptI];
-
-    pointFeatureEdgesTypes pointEdgeTypes(ptI);
-
-    const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat
-        = calcPointFeatureEdgesTypes(feMesh, featPtEdges, pointEdgeTypes);
-
-//    Info<< nl << featPt << "  " << pointEdgeTypes;
-
-    const_circulator<labelList> circ(featPtEdges);
-
-    // Loop around the edges of the feature point
-    if (circ.size()) do
-    {
-//        const edgeStatus eStatusPrev = feMesh.getEdgeStatus(circ.prev());
-        const edgeStatus eStatusCurr = feMesh.getEdgeStatus(circ());
-//        const edgeStatus eStatusNext = feMesh.getEdgeStatus(circ.next());
-
-//        Info<< "    Prev = "
-//            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusPrev]
-//            << " Curr = "
-//            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusCurr]
-////            << " Next = "
-////            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusNext]
-//            << endl;
-
-        // Get the direction in which to move the point in relation to the
-        // feature point
-        label sign = getSign(eStatusCurr);
-
-        const vector n = sharedFaceNormal(feMesh, circ(), circ.next());
-
-        const vector pointMotionDirection = sign*0.5*ppDist*n;
-
-//        Info<< "    Shared face normal      = " << n << endl;
-//        Info<< "    Direction to move point = " << pointMotionDirection
-//            << endl;
-
-        if (masterPoints.empty())
-        {
-            // Initialise with the first master point
-            Foam::point pt = featPt + pointMotionDirection;
-
-            planeDynList firstPlane;
-            firstPlane.append(autoPtr<plane>(new plane(featPt, n)));
-
-            masterPoints.append(pt);
-
-            masterPointsTypes.append
-            (
-                sign == 1
-              ? Vb::vtExternalFeaturePoint // true
-              : Vb::vtInternalFeaturePoint // false
-            );
-
-            //Info<< "    " << " " << firstPlane << endl;
-
-//            const Foam::point reflectedPoint = reflectPointInPlane
-//            (
-//                masterPoints.last(),
-//                firstPlane.last()()
-//            );
-
-            masterPointReflections.insert
-            (
-                masterPoints.size() - 1,
-                firstPlane
-            );
-        }
-//        else if
-//        (
-//            eStatusPrev == extendedFeatureEdgeMesh::INTERNAL
-//         && eStatusCurr == extendedFeatureEdgeMesh::EXTERNAL
-//        )
-//        {
-//            // Insert a new master point.
-//            Foam::point pt = featPt + pointMotionDirection;
-//
-//            planeDynList firstPlane;
-//            firstPlane.append(autoPtr<plane>(new plane(featPt, n)));
-//
-//            masterPoints.append(pt);
-//
-//            masterPointsTypes.append
-//            (
-//                sign == 1
-//              ? Vb::vtExternalFeaturePoint // true
-//              : Vb::vtInternalFeaturePoint // false
-//            );
-//
-//            masterPointReflections.insert
-//            (
-//                masterPoints.size() - 1,
-//                firstPlane
-//            );
-//        }
-//        else if
-//        (
-//            eStatusPrev == extendedFeatureEdgeMesh::EXTERNAL
-//         && eStatusCurr == extendedFeatureEdgeMesh::INTERNAL
-//        )
-//        {
-//
-//        }
-        else
-        {
-            // Just add this face contribution to the latest master point
-
-            masterPoints.last() += pointMotionDirection;
-
-            masterPointReflections[masterPoints.size() - 1].append
-            (
-                autoPtr<plane>(new plane(featPt, n))
-            );
-        }
-
-    } while (circ.circulate(CirculatorBase::CLOCKWISE));
-
-    addMasterAndSlavePoints
-    (
-        masterPoints,
-        masterPointsTypes,
-        masterPointReflections,
-        pts,
-        ptI
-    );
-}
-
-
-void Foam::conformalVoronoiMesh::createFeaturePoints(DynamicList<Vb>& pts)
-{
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-    (
-        geometryToConformTo_.features()
-    );
-
-    forAll(feMeshes, i)
-    {
-        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-
-        for
-        (
-            label ptI = feMesh.convexStart();
-            ptI < feMesh.mixedStart();
-//            ptI < feMesh.nonFeatureStart();
-            ptI++
-        )
-        {
-            createMasterAndSlavePoints(feMesh, ptI, pts);
-        }
-
-        if (foamyHexMeshControls().guardFeaturePoints())
-        {
-            for
-            (
-                //label ptI = feMesh.convexStart();
-                label ptI = feMesh.mixedStart();
-                ptI < feMesh.nonFeatureStart();
-                ptI++
-            )
-            {
-                pts.append
-                (
-                    Vb
-                    (
-                        feMesh.points()[ptI],
-                        Vb::vtConstrained
-                    )
-                );
-            }
-        }
-    }
 }
 
 
@@ -1654,71 +1065,3 @@ void Foam::conformalVoronoiMesh::createFeaturePoints(DynamicList<Vb>& pts)
 //        }
 //    }
 //}
-
-
-Foam::vector Foam::conformalVoronoiMesh::sharedFaceNormal
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const label edgeI,
-    const label nextEdgeI
-) const
-{
-    const labelList& edgeInormals = feMesh.edgeNormals()[edgeI];
-    const labelList& nextEdgeInormals = feMesh.edgeNormals()[nextEdgeI];
-
-    const vector& A1 = feMesh.normals()[edgeInormals[0]];
-    const vector& A2 = feMesh.normals()[edgeInormals[1]];
-
-    const vector& B1 = feMesh.normals()[nextEdgeInormals[0]];
-    const vector& B2 = feMesh.normals()[nextEdgeInormals[1]];
-
-//    Info<< "        A1 = " << A1 << endl;
-//    Info<< "        A2 = " << A2 << endl;
-//    Info<< "        B1 = " << B1 << endl;
-//    Info<< "        B2 = " << B2 << endl;
-
-    const scalar A1B1 = mag((A1 & B1) - 1.0);
-    const scalar A1B2 = mag((A1 & B2) - 1.0);
-    const scalar A2B1 = mag((A2 & B1) - 1.0);
-    const scalar A2B2 = mag((A2 & B2) - 1.0);
-
-//    Info<< "      A1B1 = " << A1B1 << endl;
-//    Info<< "      A1B2 = " << A1B2 << endl;
-//    Info<< "      A2B1 = " << A2B1 << endl;
-//    Info<< "      A2B2 = " << A2B2 << endl;
-
-    if (A1B1 < A1B2 && A1B1 < A2B1 && A1B1 < A2B2)
-    {
-        return 0.5*(A1 + B1);
-    }
-    else if (A1B2 < A1B1 && A1B2 < A2B1 && A1B2 < A2B2)
-    {
-        return 0.5*(A1 + B2);
-    }
-    else if (A2B1 < A1B1 && A2B1 < A1B2 && A2B1 < A2B2)
-    {
-        return 0.5*(A2 + B1);
-    }
-    else
-    {
-        return 0.5*(A2 + B2);
-    }
-}
-
-
-Foam::List<Foam::point> Foam::conformalVoronoiMesh::reflectPointInPlanes
-(
-    const Foam::point p,
-    const DynamicList<autoPtr<plane> >& planes
-) const
-{
-    List<Foam::point> reflectedPoints(planes.size());
-
-    forAll(planes, planeI)
-    {
-        reflectedPoints[planeI] = planes[planeI]().mirror(p);
-    }
-
-    return reflectedPoints;
-}
-

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H

@@ -23,6 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "indexedVertexOps.H"
+#include "indexedCellOps.H"
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 inline Foam::scalar Foam::conformalVoronoiMesh::defaultCellSize() const
@@ -40,25 +43,6 @@ inline Foam::scalar Foam::conformalVoronoiMesh::targetCellSize
 }
 
 
-inline Foam::scalar Foam::conformalVoronoiMesh::averageCellSize
-(
-    const Vertex_handle& vA,
-    const Vertex_handle& vB
-) const
-{
-    // Arithmetic mean
-    // return 0.5*(vA->targetCellSize() + vB->targetCellSize());
-
-    // Geometric mean
-    return sqrt(vA->targetCellSize()*vB->targetCellSize());
-
-    // Harmonic mean
-    // return
-    //     2.0*(vA->targetCellSize()*vB->targetCellSize())
-    //    /(vA->targetCellSize() + vB->targetCellSize());
-}
-
-
 inline Foam::scalar Foam::conformalVoronoiMesh::averageAnyCellSize
 (
     const Vertex_handle& vA,
@@ -90,7 +74,7 @@ inline Foam::scalar Foam::conformalVoronoiMesh::averageAnyCellSize
         return vB->targetCellSize();
     }
 
-    return averageCellSize(vA, vB);
+    return CGAL::indexedVertexOps::averageCellSize(vA, vB);
 }
 
 
@@ -242,7 +226,7 @@ inline void Foam::conformalVoronoiMesh::createPointPair
     const Foam::point& surfPt,
     const vector& n,
     DynamicList<Vb>& pts
-)
+) const
 {
     vector ppDistn = ppDist*n;
 
@@ -313,7 +297,7 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
     const Foam::point& surfPt,
     const vector& n,
     DynamicList<Vb>& pts
-)
+) const
 {
     vector ppDistn = ppDist*n;
 
@@ -341,6 +325,24 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
 }
 
 
+inline bool Foam::conformalVoronoiMesh::internalPointIsInside
+(
+    const Foam::point& pt
+) const
+{
+    if
+    (
+        !geometryToConformTo_.inside(pt)
+     || !geometryToConformTo_.globalBounds().contains(pt)
+    )
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
 inline bool Foam::conformalVoronoiMesh::isPointPair
 (
     const Vertex_handle& vA,
@@ -507,9 +509,9 @@ inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
     const int oppositeVertex = fit->second;
     const Cell_handle c2(c1->neighbor(oppositeVertex));
 
-    FixedList<label, 4> c1Procs(c1->processorsAttached());
+    FixedList<label, 4> c1Procs(CGAL::indexedCellOps::processorsAttached(c1));
 
-    FixedList<label, 4> c2Procs(c2->processorsAttached());
+    FixedList<label, 4> c2Procs(CGAL::indexedCellOps::processorsAttached(c2));
 
     forAll(c1Procs, aPI)
     {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -33,6 +33,8 @@ License
 #include "polyTopoChange.H"
 #include "PrintTable.H"
 #include "pointMesh.H"
+#include "indexedVertexOps.H"
+#include "DelaunayMeshTools.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -91,272 +93,9 @@ void Foam::conformalVoronoiMesh::timeCheck
 }
 
 
-void Foam::conformalVoronoiMesh::drawDelaunayCell
-(
-    Ostream& os,
-    const Cell_handle& c,
-    label offset
-) const
-{
-    // Supply offset as tet number
-    offset *= 4;
-
-    os  << "# cell index: " << label(c->cellIndex())
-        << " INT_MIN = " << INT_MIN
-        << endl;
-
-    os  << "# circumradius "
-        << mag(c->dual() - topoint(c->vertex(0)->point()))
-        << endl;
-
-    for (int i = 0; i < 4; i++)
-    {
-        os  << "# index / type / procIndex: "
-            << label(c->vertex(i)->index()) << " "
-            << label(c->vertex(i)->type()) << " "
-            << label(c->vertex(i)->procIndex())
-            << (is_infinite(c->vertex(i)) ? " # This vertex is infinite!" : "")
-            <<
-                (
-                    c->vertex(i)->uninitialised()
-                  ? " # This vertex is uninitialised!"
-                  : ""
-                )
-            << endl;
-
-        meshTools::writeOBJ(os, topoint(c->vertex(i)->point()));
-    }
-
-    os  << "f " << 1 + offset << " " << 3 + offset << " " << 2 + offset << nl
-        << "f " << 2 + offset << " " << 3 + offset << " " << 4 + offset << nl
-        << "f " << 1 + offset << " " << 4 + offset << " " << 3 + offset << nl
-        << "f " << 1 + offset << " " << 2 + offset << " " << 4 + offset << endl;
-
-//    os  << "# cicumcentre " << endl;
-
-//    meshTools::writeOBJ(os, c->dual());
-
-//    os  << "l " << 1 + offset << " " << 5 + offset << endl;
-}
-
-
-void Foam::conformalVoronoiMesh::writePoints
-(
-    const fileName& fName,
-    const Foam::indexedVertexEnum::vertexType startPointType,
-    const Foam::indexedVertexEnum::vertexType endPointType
-) const
-{
-    OFstream str(runTime_.path()/fName);
-
-    Pout<< nl << "Writing points of types:" << nl;
-
-    forAllConstIter
-    (
-        HashTable<int>,
-        Foam::indexedVertexEnum::vertexTypeNames_,
-        iter
-    )
-    {
-        if (iter() >= startPointType && iter() <= endPointType)
-        {
-            Pout<< "    " << iter.key() << nl;
-        }
-    }
-
-    Pout<< "to " << str.name() << endl;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->type() >= startPointType && vit->type() <= endPointType)
-        {
-            meshTools::writeOBJ(str, topoint(vit->point()));
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::writePoints
-(
-    const fileName& fName,
-    const Foam::indexedVertexEnum::vertexType pointType
-) const
-{
-    writePoints(fName, pointType, pointType);
-}
-
-
-void Foam::conformalVoronoiMesh::writeFixedPoints
-(
-    const fileName& fName
-) const
-{
-    OFstream str(runTime_.path()/fName);
-
-    Pout<< nl << "Writing fixed points to " << str.name() << endl;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->fixed())
-        {
-            meshTools::writeOBJ(str, topoint(vit->point()));
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::writeBoundaryPoints
-(
-    const fileName& fName
-) const
-{
-    OFstream str(runTime_.path()/fName);
-
-    Pout<< nl << "Writing boundary points to " << str.name() << endl;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (!vit->internalPoint())
-        {
-            meshTools::writeOBJ(str, topoint(vit->point()));
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::writePoints
-(
-    const fileName& fName,
-    const List<Foam::point>& points
-) const
-{
-    if (points.size())
-    {
-        OFstream str(runTime_.path()/fName);
-
-        Pout<< nl << "Writing " << points.size() << " points from pointList to "
-            << str.name() << endl;
-
-        forAll(points, p)
-        {
-            meshTools::writeOBJ(str, points[p]);
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::writePoints
-(
-    const fileName& fName,
-    const List<Vb>& points
-) const
-{
-    if (points.size())
-    {
-        OFstream str(runTime_.path()/fName);
-
-        Pout<< nl << "Writing " << points.size() << " points from pointList to "
-            << str.name() << endl;
-
-        forAll(points, p)
-        {
-            meshTools::writeOBJ(str, topoint(points[p].point()));
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::writeProcessorInterface
-(
-    const fileName& fName,
-    const faceList& faces
-) const
-{
-    OFstream str(runTime_.path()/fName);
-
-    pointField points(number_of_finite_cells(), point::max);
-
-    for
-    (
-        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if (!cit->hasFarPoint() && !is_infinite(cit))
-        {
-            points[cit->cellIndex()] = cit->dual();
-        }
-    }
-
-    meshTools::writeOBJ(str, faces, points);
-}
-
-
-void Foam::conformalVoronoiMesh::writeInternalDelaunayVertices
-(
-    const fileName& instance
-) const
-{
-    pointField internalDelaunayVertices(number_of_vertices());
-
-    label vertI = 0;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->internalPoint())
-        {
-            internalDelaunayVertices[vertI++] = topoint(vit->point());
-        }
-    }
-
-    internalDelaunayVertices.setSize(vertI);
-
-    pointIOField internalDVs
-    (
-        IOobject
-        (
-            "internalDelaunayVertices",
-            instance,
-            runTime_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        internalDelaunayVertices
-    );
-
-    Info<< nl
-        << "Writing " << internalDVs.name()
-        << " to " << internalDVs.instance()
-        << endl;
-
-    internalDVs.write();
-}
-
-
 void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 {
-    writeInternalDelaunayVertices(instance);
+    DelaunayMeshTools::writeInternalDelaunayVertices(instance, *this);
 
     // Per cell the Delaunay vertex
     labelList cellToDelaunayVertex;
@@ -935,6 +674,8 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
 
     pBufs.finishedSends();
 
+    Info<< incrIndent << indent << "    Face ordering initialised..." << endl;
+
     // Receive and calculate ordering
     bool anyChanged = false;
 
@@ -992,6 +733,8 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
         }
     }
 
+    Info<< incrIndent << indent << "    Faces matched." << endl;
+
     reduce(anyChanged, orOp<bool>());
 
     if (anyChanged)
@@ -1049,7 +792,12 @@ void Foam::conformalVoronoiMesh::writeMesh
 {
     if (foamyHexMeshControls().objOutput())
     {
-        writeObjMesh(points, faces, word(meshName + ".obj"));
+        DelaunayMeshTools::writeObjMesh
+        (
+            time().path()/word(meshName + ".obj"),
+            points,
+            faces
+        );
     }
 
     const label nInternalFaces = readLabel(patchDicts[0].lookup("startFace"));
@@ -1404,38 +1152,6 @@ void Foam::conformalVoronoiMesh::writeMesh
 }
 
 
-void Foam::conformalVoronoiMesh::writeObjMesh
-(
-    const pointField& points,
-    const faceList& faces,
-    const fileName& fName
-) const
-{
-    OFstream str(runTime_.path()/fName);
-
-    Pout<< nl << "Writing points and faces to " << str.name() << endl;
-
-    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::conformalVoronoiMesh::writeCellSizes
 (
     const fvMesh& mesh

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C

@@ -27,47 +27,6 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Triangulation>
-Foam::scalar Foam::conformalVoronoiMesh::calculateLoadUnbalance
-(
-    const Triangulation& mesh
-) const
-{
-    label nRealVertices = 0;
-
-    for
-    (
-        typename Triangulation::Finite_vertices_iterator vit
-            = mesh.finite_vertices_begin();
-        vit != mesh.finite_vertices_end();
-        ++vit
-    )
-    {
-        // Only store real vertices that are not feature vertices
-        if (vit->real() && !vit->featurePoint())
-        {
-            nRealVertices++;
-        }
-    }
-
-    scalar globalNRealVertices = returnReduce
-    (
-        nRealVertices,
-        sumOp<label>()
-    );
-
-    scalar unbalance = returnReduce
-    (
-        mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
-        maxOp<scalar>()
-    );
-
-    Info<< "    Processor unbalance " << unbalance << endl;
-
-    return unbalance;
-}
-
-
 template<class Triangulation>
 bool Foam::conformalVoronoiMesh::distributeBackground(const Triangulation& mesh)
 {
@@ -86,7 +45,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground(const Triangulation& mesh)
 
     while (true)
     {
-        scalar maxLoadUnbalance = calculateLoadUnbalance(mesh);
+        scalar maxLoadUnbalance = mesh.calculateLoadUnbalance();
 
         if
         (

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C

@@ -0,0 +1,604 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "featurePointConformer.H"
+#include "cvControls.H"
+#include "conformationSurfaces.H"
+#include "conformalVoronoiMesh.H"
+#include "cellShapeControl.H"
+#include "DelaunayMeshTools.H"
+#include "const_circulator.H"
+#include "backgroundMeshDecomposition.H"
+#include "autoPtr.H"
+#include "mapDistribute.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(featurePointConformer, 0);
+}
+
+const Foam::scalar Foam::featurePointConformer::tolParallel = 1e-3;
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+Foam::vector Foam::featurePointConformer::sharedFaceNormal
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const label edgeI,
+    const label nextEdgeI
+) const
+{
+    const labelList& edgeInormals = feMesh.edgeNormals()[edgeI];
+    const labelList& nextEdgeInormals = feMesh.edgeNormals()[nextEdgeI];
+
+    const vector& A1 = feMesh.normals()[edgeInormals[0]];
+    const vector& A2 = feMesh.normals()[edgeInormals[1]];
+
+    const vector& B1 = feMesh.normals()[nextEdgeInormals[0]];
+    const vector& B2 = feMesh.normals()[nextEdgeInormals[1]];
+
+//    Info<< "        A1 = " << A1 << endl;
+//    Info<< "        A2 = " << A2 << endl;
+//    Info<< "        B1 = " << B1 << endl;
+//    Info<< "        B2 = " << B2 << endl;
+
+    const scalar A1B1 = mag((A1 & B1) - 1.0);
+    const scalar A1B2 = mag((A1 & B2) - 1.0);
+    const scalar A2B1 = mag((A2 & B1) - 1.0);
+    const scalar A2B2 = mag((A2 & B2) - 1.0);
+
+//    Info<< "      A1B1 = " << A1B1 << endl;
+//    Info<< "      A1B2 = " << A1B2 << endl;
+//    Info<< "      A2B1 = " << A2B1 << endl;
+//    Info<< "      A2B2 = " << A2B2 << endl;
+
+    if (A1B1 < A1B2 && A1B1 < A2B1 && A1B1 < A2B2)
+    {
+        return 0.5*(A1 + B1);
+    }
+    else if (A1B2 < A1B1 && A1B2 < A2B1 && A1B2 < A2B2)
+    {
+        return 0.5*(A1 + B2);
+    }
+    else if (A2B1 < A1B1 && A2B1 < A1B2 && A2B1 < A2B2)
+    {
+        return 0.5*(A2 + B1);
+    }
+    else
+    {
+        return 0.5*(A2 + B2);
+    }
+}
+
+
+Foam::label Foam::featurePointConformer::getSign
+(
+    const extendedFeatureEdgeMesh::edgeStatus eStatus
+) const
+{
+    if (eStatus == extendedFeatureEdgeMesh::EXTERNAL)
+    {
+        return -1;
+    }
+    else if (eStatus == extendedFeatureEdgeMesh::INTERNAL)
+    {
+        return 1;
+    }
+
+    return 0;
+}
+
+
+void Foam::featurePointConformer::addMasterAndSlavePoints
+(
+    const DynamicList<Foam::point>& masterPoints,
+    const DynamicList<Foam::indexedVertexEnum::vertexType>& masterPointsTypes,
+    const Map<DynamicList<autoPtr<plane> > >& masterPointReflections,
+    DynamicList<Vb>& pts,
+    const label ptI
+) const
+{
+    typedef DynamicList<autoPtr<plane> >        planeDynList;
+    typedef Foam::indexedVertexEnum::vertexType vertexType;
+
+    forAll(masterPoints, pI)
+    {
+        // Append master to the list of points
+
+        const Foam::point& masterPt = masterPoints[pI];
+        const vertexType masterType = masterPointsTypes[pI];
+
+//        Info<< "    Master = " << masterPt << endl;
+
+        pts.append
+        (
+            Vb
+            (
+                masterPt,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                masterType,
+                Pstream::myProcNo()
+            )
+        );
+
+//        const label masterIndex = pts[pts.size() - 1].index();
+
+        //meshTools::writeOBJ(strMasters, masterPt);
+
+        const planeDynList& masterPointPlanes = masterPointReflections[pI];
+
+        forAll(masterPointPlanes, planeI)
+        {
+            // Reflect master points in the planes and insert the slave points
+
+            const plane& reflPlane = masterPointPlanes[planeI]();
+
+            const Foam::point slavePt = reflPlane.mirror(masterPt);
+
+//            Info<< "        Slave " << planeI << " = " << slavePt << endl;
+
+            const vertexType slaveType =
+            (
+                masterType == Vb::vtInternalFeaturePoint
+              ? Vb::vtExternalFeaturePoint // true
+              : Vb::vtInternalFeaturePoint // false
+            );
+
+            pts.append
+            (
+                Vb
+                (
+                    slavePt,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    slaveType,
+                    Pstream::myProcNo()
+                )
+            );
+
+
+            //meshTools::writeOBJ(strSlaves, slavePt);
+        }
+    }
+}
+
+
+void Foam::featurePointConformer::createMasterAndSlavePoints
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const label ptI,
+    DynamicList<Vb>& pts
+) const
+{
+    typedef DynamicList<autoPtr<plane> >          planeDynList;
+    typedef indexedVertexEnum::vertexType         vertexType;
+    typedef extendedFeatureEdgeMesh::edgeStatus   edgeStatus;
+
+    const Foam::point& featPt = feMesh.points()[ptI];
+
+    if
+    (
+        (
+            Pstream::parRun()
+         && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt)
+        )
+     || geometryToConformTo_.outside(featPt)
+    )
+    {
+        return;
+    }
+
+    const scalar ppDist = foamyHexMesh_.pointPairDistance(featPt);
+
+    // Maintain a list of master points and the planes to relect them in
+    DynamicList<Foam::point> masterPoints;
+    DynamicList<vertexType> masterPointsTypes;
+    Map<planeDynList> masterPointReflections;
+
+    const labelList& featPtEdges = feMesh.featurePointEdges()[ptI];
+
+    pointFeatureEdgesTypes pointEdgeTypes(feMesh, ptI);
+
+    const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat =
+        pointEdgeTypes.calcPointFeatureEdgesTypes();
+
+//    Info<< nl << featPt << "  " << pointEdgeTypes;
+
+    const_circulator<labelList> circ(featPtEdges);
+
+    // Loop around the edges of the feature point
+    if (circ.size()) do
+    {
+//        const edgeStatus eStatusPrev = feMesh.getEdgeStatus(circ.prev());
+        const edgeStatus eStatusCurr = feMesh.getEdgeStatus(circ());
+//        const edgeStatus eStatusNext = feMesh.getEdgeStatus(circ.next());
+
+//        Info<< "    Prev = "
+//            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusPrev]
+//            << " Curr = "
+//            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusCurr]
+////            << " Next = "
+////            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusNext]
+//            << endl;
+
+        // Get the direction in which to move the point in relation to the
+        // feature point
+        label sign = getSign(eStatusCurr);
+
+        const vector n = sharedFaceNormal(feMesh, circ(), circ.next());
+
+        const vector pointMotionDirection = sign*0.5*ppDist*n;
+
+//        Info<< "    Shared face normal      = " << n << endl;
+//        Info<< "    Direction to move point = " << pointMotionDirection
+//            << endl;
+
+        if (masterPoints.empty())
+        {
+            // Initialise with the first master point
+            Foam::point pt = featPt + pointMotionDirection;
+
+            planeDynList firstPlane;
+            firstPlane.append(autoPtr<plane>(new plane(featPt, n)));
+
+            masterPoints.append(pt);
+
+            masterPointsTypes.append
+            (
+                sign == 1
+              ? Vb::vtExternalFeaturePoint // true
+              : Vb::vtInternalFeaturePoint // false
+            );
+
+            //Info<< "    " << " " << firstPlane << endl;
+
+//            const Foam::point reflectedPoint = reflectPointInPlane
+//            (
+//                masterPoints.last(),
+//                firstPlane.last()()
+//            );
+
+            masterPointReflections.insert
+            (
+                masterPoints.size() - 1,
+                firstPlane
+            );
+        }
+//        else if
+//        (
+//            eStatusPrev == extendedFeatureEdgeMesh::INTERNAL
+//         && eStatusCurr == extendedFeatureEdgeMesh::EXTERNAL
+//        )
+//        {
+//            // Insert a new master point.
+//            Foam::point pt = featPt + pointMotionDirection;
+//
+//            planeDynList firstPlane;
+//            firstPlane.append(autoPtr<plane>(new plane(featPt, n)));
+//
+//            masterPoints.append(pt);
+//
+//            masterPointsTypes.append
+//            (
+//                sign == 1
+//              ? Vb::vtExternalFeaturePoint // true
+//              : Vb::vtInternalFeaturePoint // false
+//            );
+//
+//            masterPointReflections.insert
+//            (
+//                masterPoints.size() - 1,
+//                firstPlane
+//            );
+//        }
+//        else if
+//        (
+//            eStatusPrev == extendedFeatureEdgeMesh::EXTERNAL
+//         && eStatusCurr == extendedFeatureEdgeMesh::INTERNAL
+//        )
+//        {
+//
+//        }
+        else
+        {
+            // Just add this face contribution to the latest master point
+
+            masterPoints.last() += pointMotionDirection;
+
+            masterPointReflections[masterPoints.size() - 1].append
+            (
+                autoPtr<plane>(new plane(featPt, n))
+            );
+        }
+
+    } while (circ.circulate(CirculatorBase::CLOCKWISE));
+
+    addMasterAndSlavePoints
+    (
+        masterPoints,
+        masterPointsTypes,
+        masterPointReflections,
+        pts,
+        ptI
+    );
+}
+
+
+void Foam::featurePointConformer::createMixedFeaturePoints
+(
+    DynamicList<Vb>& pts
+) const
+{
+    if (foamyHexMeshControls_.mixedFeaturePointPPDistanceCoeff() < 0)
+    {
+        Info<< nl
+            << "Skipping specialised handling for mixed feature points" << endl;
+
+        return;
+    }
+
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh = feMeshes[i];
+        const labelListList& pointsEdges = feMesh.pointEdges();
+        const pointField& points = feMesh.points();
+
+        for
+        (
+            label ptI = feMesh.mixedStart();
+            ptI < feMesh.nonFeatureStart();
+            ptI++
+        )
+        {
+            const Foam::point& featPt = points[ptI];
+
+            if
+            (
+                Pstream::parRun()
+             && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt))
+            {
+                continue;
+            }
+
+            const labelList& pEds = pointsEdges[ptI];
+
+            pointFeatureEdgesTypes pFEdgeTypes(feMesh, ptI);
+
+            const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat =
+                pFEdgeTypes.calcPointFeatureEdgesTypes();
+
+            bool specialisedSuccess = false;
+
+            if (foamyHexMeshControls_.specialiseFeaturePoints())
+            {
+                specialisedSuccess =
+                    createSpecialisedFeaturePoint
+                    (
+                        feMesh, pEds, pFEdgeTypes, allEdStat, ptI, pts
+                    );
+            }
+
+            if (!specialisedSuccess && foamyHexMeshControls_.edgeAiming())
+            {
+                // Specialisations available for some mixed feature points.  For
+                // non-specialised feature points, inserting mixed internal and
+                // external edge groups at feature point.
+
+                // Skipping unsupported mixed feature point types
+//                bool skipEdge = false;
+//
+//                forAll(pEds, e)
+//                {
+//                    const label edgeI = pEds[e];
+//
+//                    const extendedFeatureEdgeMesh::edgeStatus edStatus
+//                        = feMesh.getEdgeStatus(edgeI);
+//
+//                    if
+//                    (
+//                        edStatus == extendedFeatureEdgeMesh::OPEN
+//                     || edStatus == extendedFeatureEdgeMesh::MULTIPLE
+//                    )
+//                    {
+//                        Info<< "Edge type " << edStatus
+//                            << " found for mixed feature point " << ptI
+//                            << ". Not supported."
+//                            << endl;
+//
+//                        skipEdge = true;
+//                    }
+//                }
+//
+//                if (skipEdge)
+//                {
+//                    Info<< "Skipping point " << ptI << nl << endl;
+//
+//                    continue;
+//                }
+
+//                createFeaturePoints(feMesh, ptI, pts, types);
+
+                const Foam::point& pt = points[ptI];
+
+                const scalar edgeGroupDistance =
+                    foamyHexMesh_.mixedFeaturePointDistance(pt);
+
+                forAll(pEds, e)
+                {
+                    const label edgeI = pEds[e];
+
+                    const Foam::point edgePt =
+                        pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
+
+                    const pointIndexHit edgeHit(true, edgePt, edgeI);
+
+                    foamyHexMesh_.createEdgePointGroup(feMesh, edgeHit, pts);
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::featurePointConformer::createFeaturePoints(DynamicList<Vb>& pts)
+{
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+
+        for
+        (
+            label ptI = feMesh.convexStart();
+            ptI < feMesh.mixedStart();
+//            ptI < feMesh.nonFeatureStart();
+            ptI++
+        )
+        {
+            createMasterAndSlavePoints(feMesh, ptI, pts);
+        }
+
+        if (foamyHexMeshControls_.guardFeaturePoints())
+        {
+            for
+            (
+                //label ptI = feMesh.convexStart();
+                label ptI = feMesh.mixedStart();
+                ptI < feMesh.nonFeatureStart();
+                ptI++
+            )
+            {
+                pts.append
+                (
+                    Vb
+                    (
+                        feMesh.points()[ptI],
+                        Vb::vtConstrained
+                    )
+                );
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::featurePointConformer::featurePointConformer
+(
+    const conformalVoronoiMesh& foamyHexMesh
+)
+:
+    foamyHexMesh_(foamyHexMesh),
+    foamyHexMeshControls_(foamyHexMesh.foamyHexMeshControls()),
+    geometryToConformTo_(foamyHexMesh.geometryToConformTo()),
+    featurePointVertices_()
+{
+    Info<< nl
+        << "Conforming to feature points" << endl;
+
+    Info<< incrIndent
+        << indent << "Circulating edges is: "
+        << foamyHexMeshControls_.circulateEdges().asText() << nl
+        << indent << "Guarding feature points is: "
+        << foamyHexMeshControls_.guardFeaturePoints().asText() << nl
+        << indent << "Snapping to feature points is: "
+        << foamyHexMeshControls_.snapFeaturePoints().asText() << nl
+        << indent << "Specialising feature points is: "
+        << foamyHexMeshControls_.specialiseFeaturePoints().asText()
+        << decrIndent
+        << endl;
+
+    DynamicList<Vb> pts;
+
+    createFeaturePoints(pts);
+
+    createMixedFeaturePoints(pts);
+
+    // Points added using the createEdgePointGroup function will be labelled as
+    // internal/external feature edges. Relabel them as feature points,
+    // otherwise they are inserted as both feature points and surface points.
+    forAll(pts, pI)
+    {
+        Vb& pt = pts[pI];
+
+        //if (pt.featureEdgePoint())
+        {
+            if (pt.internalBoundaryPoint())
+            {
+                pt.type() = Vb::vtInternalFeaturePoint;
+            }
+            else if (pt.externalBoundaryPoint())
+            {
+                pt.type() = Vb::vtExternalFeaturePoint;
+            }
+        }
+    }
+
+    if (foamyHexMeshControls_.objOutput())
+    {
+        DelaunayMeshTools::writeOBJ("featureVertices.obj", pts);
+    }
+
+    featurePointVertices_.transfer(pts);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::featurePointConformer::~featurePointConformer()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::featurePointConformer::distribute
+(
+    const backgroundMeshDecomposition& decomposition
+)
+{
+    // Distribute points to their appropriate processor
+    decomposition.distributePoints(featurePointVertices_);
+
+    // Update the processor indices of the points to the new processor number
+    forAll(featurePointVertices_, vI)
+    {
+        featurePointVertices_[vI].procIndex() = Pstream::myProcNo();
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::featurePointConformer
+
+Description
+    The Delaunay vertices required to conform to a feature point can be
+    determined upon initialisation because the feature points are fixed and
+    do not change throughout the meshing process.
+
+SourceFiles
+    featurePointConformerI.H
+    featurePointConformer.C
+    featurePointConformerSpecialisations.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef featurePointConformer_H
+#define featurePointConformer_H
+
+#include "CGALTriangulation3Ddefs.H"
+#include "vector.H"
+#include "DynamicList.H"
+#include "List.H"
+#include "extendedFeatureEdgeMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class conformalVoronoiMesh;
+class cvControls;
+class conformationSurfaces;
+class pointFeatureEdgesTypes;
+class backgroundMeshDecomposition;
+
+/*---------------------------------------------------------------------------*\
+                   Class featurePointConformer Declaration
+\*---------------------------------------------------------------------------*/
+
+class featurePointConformer
+{
+    // Static data
+
+        //- Tolerance within which two lines are said to be parallel.
+        static const scalar tolParallel;
+
+
+    // Private data
+
+        //- Reference to the mesher.
+        const conformalVoronoiMesh& foamyHexMesh_;
+
+        //- Reference to the mesher controls.
+        const cvControls& foamyHexMeshControls_;
+
+        //- Reference to the conformation surfaces.
+        const conformationSurfaces& geometryToConformTo_;
+
+        //- Store the feature constraining points, to be reinserted after a
+        //  triangulation clear.
+        List<Vb> featurePointVertices_;
+
+
+    // Private Member Functions
+
+        //- Calculate the shared face normal between two edges geometrically.
+        vector sharedFaceNormal
+        (
+            const extendedFeatureEdgeMesh& feMesh,
+            const label edgeI,
+            const label nextEdgeI
+        ) const;
+
+        label getSign(const extendedFeatureEdgeMesh::edgeStatus eStatus) const;
+
+        bool createSpecialisedFeaturePoint
+        (
+            const extendedFeatureEdgeMesh& feMesh,
+            const labelList& pEds,
+            const pointFeatureEdgesTypes& pFEdgesTypes,
+            const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
+            const label ptI,
+            DynamicList<Vb>& pts
+        ) const;
+
+
+        void addMasterAndSlavePoints
+        (
+            const DynamicList<point>& masterPoints,
+            const DynamicList<indexedVertexEnum::vertexType>& masterPointsTypes,
+            const Map<DynamicList<autoPtr<plane> > >& masterPointReflections,
+            DynamicList<Vb>& pts,
+            const label ptI
+        ) const;
+
+        //- Helper function for conforming to feature points
+        void createMasterAndSlavePoints
+        (
+            const extendedFeatureEdgeMesh& feMesh,
+            const label ptI,
+            DynamicList<Vb>& pts
+        ) const;
+
+        void createMixedFeaturePoints(DynamicList<Vb>& pts) const;
+
+        //- Create the points that will conform to the feature
+        void createFeaturePoints(DynamicList<Vb>& pts);
+
+        //- Disallow default bitwise copy construct
+        featurePointConformer(const featurePointConformer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const featurePointConformer&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("featurePointConformer");
+
+
+    // Constructors
+
+        //- Construct from components
+        explicit featurePointConformer
+        (
+            const conformalVoronoiMesh& foamyHexMesh
+        );
+
+
+    //- Destructor
+    ~featurePointConformer();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the feature point vertices for insertion into the
+            //  triangulation.
+            inline const List<Vb>& featurePointVertices() const;
+
+
+        // Edit
+
+            //- Distribute the feature point vertices according to the
+            //  supplied background mesh
+            void distribute(const backgroundMeshDecomposition& decomposition);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "featurePointConformerI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H

@@ -0,0 +1,32 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+const Foam::List<Vb>& Foam::featurePointConformer::featurePointVertices() const
+{
+    return featurePointVertices_;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C

@@ -23,39 +23,17 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "conformalVoronoiMesh.H"
+#include "featurePointConformer.H"
 #include "vectorTools.H"
+#include "pointFeatureEdgesTypes.H"
+#include "conformalVoronoiMesh.H"
+#include "pointConversion.H"
 
 using namespace Foam::vectorTools;
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-Foam::List<Foam::extendedFeatureEdgeMesh::edgeStatus>
+bool Foam::featurePointConformer::createSpecialisedFeaturePoint
-Foam::conformalVoronoiMesh::calcPointFeatureEdgesTypes
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const labelList& pEds,
-    pointFeatureEdgesTypes& pFEdgeTypes
-) const
-{
-    List<extendedFeatureEdgeMesh::edgeStatus> allEdStat(pEds.size());
-
-    forAll(pEds, i)
-    {
-        label edgeI = pEds[i];
-
-        extendedFeatureEdgeMesh::edgeStatus& eS = allEdStat[i];
-
-        eS = feMesh.getEdgeStatus(edgeI);
-
-        pFEdgeTypes(eS)++;
-    }
-
-    return allEdStat;
-}
-
-
-bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
 (
     const extendedFeatureEdgeMesh& feMesh,
     const labelList& pEds,
@@ -63,7 +41,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
     const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
     const label ptI,
     DynamicList<Vb>& pts
-)
+) const
 {
     if
     (
@@ -81,20 +59,24 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
      && pEds.size() == 3
     )
     {
-        Info<< "nExternal == 2 && nInternal == 1" << endl;
+        if (debug) Info<< "nExternal == 2 && nInternal == 1" << endl;
 
         const Foam::point& featPt = feMesh.points()[ptI];
 
-        if (!positionOnThisProc(featPt))
+        if
+        (
+            Pstream::parRun()
+         && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt)
+        )
         {
             return false;
         }
 
-        label nVert = number_of_vertices();
+        label nVert = foamyHexMesh_.number_of_vertices();
 
         const label initialNumOfPoints = pts.size();
 
-        const scalar ppDist = pointPairDistance(featPt);
+        const scalar ppDist = foamyHexMesh_.pointPairDistance(featPt);
 
         const vectorField& normals = feMesh.normals();
 
@@ -204,8 +186,8 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
         Foam::point externalPtD;
         Foam::point externalPtE;
 
-        vector convexEdgePlaneCNormal(0,0,0);
+        vector convexEdgePlaneCNormal(vector::zero);
-        vector convexEdgePlaneDNormal(0,0,0);
+        vector convexEdgePlaneDNormal(vector::zero);
 
         const labelList& concaveEdgeNormals = edgeNormals[concaveEdgeI];
         const labelList& convexEdgeANormals = edgeNormals[convexEdgesI[0]];
@@ -223,8 +205,11 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 const vector& convexNormal
                     = normals[convexEdgeANormals[edgeAnormalI]];
 
+                if (debug)
+                {
                     Info<< "Angle between vectors = "
                         << degAngleBetween(concaveNormal, convexNormal) << endl;
+                }
 
                 // Need a looser tolerance, because sometimes adjacent triangles
                 // on the same surface will be slightly out of alignment.
@@ -241,8 +226,11 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 const vector& convexNormal
                     = normals[convexEdgeBNormals[edgeBnormalI]];
 
+                if (debug)
+                {
                     Info<< "Angle between vectors = "
                         << degAngleBetween(concaveNormal, convexNormal) << endl;
+                }
 
                 // Need a looser tolerance, because sometimes adjacent triangles
                 // on the same surface will be slightly out of alignment.
@@ -349,7 +337,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
           + radAngleBetween(concaveEdgePlaneANormal, concaveEdgePlaneBNormal)
         );
 
-        if (totalAngle > foamyHexMeshControls().maxQuadAngle())
+        if (totalAngle > foamyHexMeshControls_.maxQuadAngle())
         {
             // Add additional mitreing points
             //scalar angleSign = 1.0;
@@ -419,21 +407,28 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
      && pFEdgesTypes[extendedFeatureEdgeMesh::INTERNAL] == 2
      && pEds.size() == 3
     )
+    {
+        if (debug)
         {
             Info<< "nExternal == 1 && nInternal == 2" << endl;
+        }
 
         const Foam::point& featPt = feMesh.points()[ptI];
 
-        if (!positionOnThisProc(featPt))
+        if
+        (
+            Pstream::parRun()
+         && !foamyHexMesh_.decomposition().positionOnThisProcessor(featPt)
+        )
         {
             return false;
         }
 
-        label nVert = number_of_vertices();
+        label nVert = foamyHexMesh_.number_of_vertices();
 
         const label initialNumOfPoints = pts.size();
 
-        const scalar ppDist = pointPairDistance(featPt);
+        const scalar ppDist = foamyHexMesh_.pointPairDistance(featPt);
 
         const vectorField& normals = feMesh.normals();
 
@@ -543,8 +538,8 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
         Foam::point externalPtD;
         Foam::point externalPtE;
 
-        vector concaveEdgePlaneCNormal(0,0,0);
+        vector concaveEdgePlaneCNormal(vector::zero);
-        vector concaveEdgePlaneDNormal(0,0,0);
+        vector concaveEdgePlaneDNormal(vector::zero);
 
         const labelList& convexEdgeNormals = edgeNormals[convexEdgeI];
         const labelList& concaveEdgeANormals = edgeNormals[concaveEdgesI[0]];
@@ -562,8 +557,11 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 const vector& concaveNormal
                     = normals[concaveEdgeANormals[edgeAnormalI]];
 
+                if (debug)
+                {
                     Info<< "Angle between vectors = "
                         << degAngleBetween(convexNormal, concaveNormal) << endl;
+                }
 
                 // Need a looser tolerance, because sometimes adjacent triangles
                 // on the same surface will be slightly out of alignment.
@@ -580,8 +578,11 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 const vector& concaveNormal
                     = normals[concaveEdgeBNormals[edgeBnormalI]];
 
+                if (debug)
+                {
                     Info<< "Angle between vectors = "
                         << degAngleBetween(convexNormal, concaveNormal) << endl;
+                }
 
                 // Need a looser tolerance, because sometimes adjacent triangles
                 // on the same surface will be slightly out of alignment.
@@ -691,7 +692,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
           + radAngleBetween(convexEdgePlaneANormal, convexEdgePlaneBNormal)
         );
 
-        if (totalAngle > foamyHexMeshControls().maxQuadAngle())
+        if (totalAngle > foamyHexMeshControls_.maxQuadAngle())
         {
             // Add additional mitreing points
             //scalar angleSign = 1.0;
@@ -751,6 +752,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 Info<< "Point " << ptI << " "
                     << indexedVertexEnum::vertexTypeNames_[pts[ptI].type()]
                     << " : ";
+
                 meshTools::writeOBJ(Info, topoint(pts[ptI].point()));
             }
         }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/pointFeatureEdgesTypes.C

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointFeatureEdgesTypes.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pointFeatureEdgesTypes::pointFeatureEdgesTypes
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const label pointLabel
+)
+:
+    HashTable<label, extendedFeatureEdgeMesh::edgeStatus>(),
+    feMesh_(feMesh),
+    pointLabel_(pointLabel)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::pointFeatureEdgesTypes::~pointFeatureEdgesTypes()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::List<Foam::extendedFeatureEdgeMesh::edgeStatus>
+Foam::pointFeatureEdgesTypes::calcPointFeatureEdgesTypes()
+{
+    const labelList& pEds = feMesh_.pointEdges()[pointLabel_];
+
+    List<extendedFeatureEdgeMesh::edgeStatus> allEdStat(pEds.size());
+
+    forAll(pEds, i)
+    {
+        label edgeI = pEds[i];
+
+        extendedFeatureEdgeMesh::edgeStatus& eS = allEdStat[i];
+
+        eS = feMesh_.getEdgeStatus(edgeI);
+
+        this->operator()(eS)++;
+    }
+
+    return allEdStat;
+}
+
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const pointFeatureEdgesTypes& p
+)
+{
+    os  << "Point = " << p.pointLabel_ << endl;
+
+    for
+    (
+        HashTable<label, extendedFeatureEdgeMesh::edgeStatus>
+            ::const_iterator iter = p.cbegin();
+        iter != p.cend();
+        ++iter
+    )
+    {
+        os  << "    "
+            << extendedFeatureEdgeMesh::edgeStatusNames_[iter.key()]
+            << " = "
+            << iter()
+            << endl;
+    }
+
+    return os;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/pointFeatureEdgesTypes.H

@@ -27,11 +27,18 @@ Class
 Description
     Holds information on the types of feature edges attached to feature points.
 
+SourceFiles
+    pointFeatureEdgesTypes.C
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef pointFeatureEdgesTypes_H
 #define pointFeatureEdgesTypes_H
 
+#include "HashTable.H"
+#include "extendedFeatureEdgeMesh.H"
+#include "List.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -46,38 +53,45 @@ class pointFeatureEdgesTypes
 :
     public HashTable<label, extendedFeatureEdgeMesh::edgeStatus>
 {
+    // Private data
+
+        //- Reference to the feature edge mesh
+        const extendedFeatureEdgeMesh& feMesh_;
+
+        //- label of the point
         label pointLabel_;
 
+
 public:
 
-    pointFeatureEdgesTypes(const label pointLabel)
+    // Constructors
-    :
-        HashTable<label, extendedFeatureEdgeMesh::edgeStatus>(),
-        pointLabel_(pointLabel)
-    {}
 
-
+        //- Construct from components
-    friend Ostream& operator<<(Ostream& os, const pointFeatureEdgesTypes& p)
+        pointFeatureEdgesTypes
-    {
-        os  << "Point = " << p.pointLabel_ << endl;
-
-        for
         (
-            HashTable<label, extendedFeatureEdgeMesh::edgeStatus>
+            const extendedFeatureEdgeMesh& feMesh,
-                ::const_iterator iter = p.cbegin();
+            const label pointLabel
-            iter != p.cend();
+        );
-            ++iter
-        )
-        {
-            os  << "    "
-                << extendedFeatureEdgeMesh::edgeStatusNames_[iter.key()]
-                << " = "
-                << iter()
-                << endl;
-        }
 
-        return os;
+
-    }
+    //- Destructor
+    ~pointFeatureEdgesTypes();
+
+
+    // Member Functions
+
+        //- Fill the pointFeatureEdgesType class with the types of feature
+        //  edges that are attached to the point.
+        List<extendedFeatureEdgeMesh::edgeStatus> calcPointFeatureEdgesTypes();
+
+
+    // Info
+
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const pointFeatureEdgesTypes& p
+        );
 };
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H

@@ -189,11 +189,6 @@ public:
         //- Does the Dual vertex form part of a processor patch
         inline bool parallelDualVertex() const;
 
-        //- Does the Dual vertex form part of a processor patch
-        inline Foam::label dualVertexMasterProc() const;
-
-        inline Foam::FixedList<Foam::label, 4> processorsAttached() const;
-
         //- Using the globalIndex object, return a list of four (sorted) global
         //  Delaunay vertex indices that uniquely identify this tet in parallel
         inline Foam::tetCell vertexGlobalIndices

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H

@@ -265,60 +265,6 @@ inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
 }
 
 
-template<class Gt, class Cb>
-inline Foam::label CGAL::indexedCell<Gt, Cb>::dualVertexMasterProc() const
-{
-    if (!parallelDualVertex())
-    {
-        return -1;
-    }
-
-    // The master processor is the lowest numbered of the four on this tet.
-
-    int masterProc = Foam::Pstream::nProcs() + 1;
-
-    for (int i = 0; i < 4; i++)
-    {
-        if (this->vertex(i)->referred())
-        {
-            masterProc = min(masterProc, this->vertex(i)->procIndex());
-        }
-        else
-        {
-            masterProc = min(masterProc, Foam::Pstream::myProcNo());
-        }
-    }
-
-    return masterProc;
-}
-
-
-template<class Gt, class Cb>
-inline Foam::FixedList<Foam::label, 4>
-CGAL::indexedCell<Gt, Cb>::processorsAttached() const
-{
-    if (!parallelDualVertex())
-    {
-        return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
-    }
-
-    Foam::FixedList<Foam::label, 4> procsAttached
-    (
-        Foam::Pstream::myProcNo()
-    );
-
-    for (int i = 0; i < 4; i++)
-    {
-        if (this->vertex(i)->referred())
-        {
-            procsAttached[i] = this->vertex(i)->procIndex();
-        }
-    }
-
-    return procsAttached;
-}
-
-
 template<class Gt, class Cb>
 inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
 (

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellOps.H

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    indexedCellOps
+
+Description
+
+SourceFiles
+    indexedCellOpsTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedCellOps_H
+#define indexedCellOps_H
+
+#include "label.H"
+#include "FixedList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                     Namespace indexedCellOps Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace indexedCellOps
+{
+
+//- Does the Dual vertex form part of a processor patch
+template<typename CellType>
+Foam::label dualVertexMasterProc(const CellType& c);
+
+template<typename CellType>
+Foam::FixedList<Foam::label, 4> processorsAttached(const CellType& c);
+
+} // End namespace indexedCellOps
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "indexedCellOpsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellOpsTemplates.C

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "indexedCellOps.H"
+#include "Pstream.H"
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<typename CellType>
+Foam::label CGAL::indexedCellOps::dualVertexMasterProc(const CellType& c)
+{
+    if (!c->parallelDualVertex())
+    {
+        return -1;
+    }
+
+    // The master processor is the lowest numbered of the four on this tet.
+
+    int masterProc = Foam::Pstream::nProcs() + 1;
+
+    for (Foam::label vI = 0; vI < 4; ++vI)
+    {
+        if (c->vertex(vI)->referred())
+        {
+            masterProc = min(masterProc, c->vertex(vI)->procIndex());
+        }
+        else
+        {
+            masterProc = min(masterProc, Foam::Pstream::myProcNo());
+        }
+    }
+
+    return masterProc;
+}
+
+
+template<typename CellType>
+Foam::FixedList<Foam::label, 4>
+CGAL::indexedCellOps::processorsAttached(const CellType& c)
+{
+    Foam::FixedList<Foam::label, 4> procsAttached(Foam::Pstream::myProcNo());
+
+    if (!c->parallelDualVertex())
+    {
+        return procsAttached;
+    }
+
+    for (Foam::label vI = 0; vI < 4; ++vI)
+    {
+        if (c->vertex(vI)->referred())
+        {
+            procsAttached[vI] = c->vertex(vI)->procIndex();
+        }
+    }
+
+    return procsAttached;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H

@@ -190,8 +190,6 @@ public:
 
         inline Foam::scalar targetCellSize() const;
 
-        inline bool uninitialised() const;
-
         //- Is point a far-point
         inline bool farPoint() const;
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H

@@ -211,13 +211,6 @@ inline Foam::scalar CGAL::indexedVertex<Gt, Vb>::targetCellSize() const
 }
 
 
-template<class Gt, class Vb>
-inline bool CGAL::indexedVertex<Gt, Vb>::uninitialised() const
-{
-    return type_ == vtUnassigned;
-}
-
-
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::farPoint() const
 {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexOps.H

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    indexedVertexOps
+
+Description
+
+SourceFiles
+    indexedVertexOpsTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedVertexOps_H
+#define indexedVertexOps_H
+
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+                    Namespace indexedVertexOps Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace indexedVertexOps
+{
+
+//- Return the target cell size from that stored on a pair of Delaunay vertices,
+//  using a mean function.
+template<typename VertexType>
+Foam::scalar averageCellSize(const VertexType& vA, const VertexType& vB);
+
+
+template<typename VertexType>
+inline bool uninitialised(const VertexType& v);
+
+} // End namespace indexedVertexOps
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "indexedVertexOpsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexOpsTemplates.C

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "indexedVertexOps.H"
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<typename VertexType>
+Foam::scalar CGAL::indexedVertexOps::averageCellSize
+(
+    const VertexType& vA,
+    const VertexType& vB
+)
+{
+    // Arithmetic mean
+    // return 0.5*(vA->targetCellSize() + vB->targetCellSize());
+
+    // Geometric mean
+    return sqrt(vA->targetCellSize()*vB->targetCellSize());
+
+    // Harmonic mean
+    // return
+    //     2.0*(vA->targetCellSize()*vB->targetCellSize())
+    //    /(vA->targetCellSize() + vB->targetCellSize());
+}
+
+
+template<typename VertexType>
+inline bool CGAL::indexedVertexOps::uninitialised(const VertexType& v)
+{
+    return v->type() == Foam::indexedVertexEnum::vtUnassigned;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointConversion.H

@@ -25,8 +25,7 @@ Class
     pointConversion
 
 Description
-
+    Conversion functions between point (Foam::) and Point (CGAL::)
-    Conversion functions between point (FOAM::) and Point (CGAL)
 
 \*---------------------------------------------------------------------------*/
 
@@ -61,7 +60,6 @@ namespace Foam
         return reinterpret_cast<pointFromPoint>(P);
     }
 
-    template<typename Point>
     inline PointFrompoint toPoint(const Foam::point& p)
     {
         return reinterpret_cast<PointFrompoint>(p);
@@ -80,14 +78,33 @@ namespace Foam
         );
     }
 
-    template<typename Point>
+    inline PointFrompoint toPoint(const Foam::point& p)
-    inline Point toPoint(const Foam::point& p)
     {
-        return Point(p.x(), p.y(), p.z());
+        return PointFrompoint(p.x(), p.y(), p.z());
     }
 
 #endif
 
+    //- Specialisation for indexedVertex.
+    template<>
+    inline pointFromPoint topoint<CGAL::indexedVertex<K> >
+    (
+        const CGAL::indexedVertex<K>& P
+    )
+    {
+        return topoint(P.point());
+    }
+
+    //- Specialisation for Foam::point. Used only as a dummy.
+    template<>
+    inline pointFromPoint topoint<Foam::point>
+    (
+        const Foam::point& P
+    )
+    {
+        return P;
+    }
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -597,10 +597,11 @@ bool Foam::conformationSurfaces::outside
 }
 
 
-Foam::Field<bool> Foam::conformationSurfaces::wellInside
+Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
 (
     const pointField& samplePts,
-    const scalarField& testDistSqr
+    const scalarField& testDistSqr,
+    const bool testForInside
 ) const
 {
     List<List<volumeType> > surfaceVolumeTests
@@ -622,43 +623,16 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 
         if (normalVolumeTypes_[regionI] != extendedFeatureEdgeMesh::BOTH)
         {
-//            if (surface.hasVolumeType())
-//            {
-//                List<List<pointIndexHit> > info;
-//
-//                // Count number of intersections
-//                surface.findLineAll
-//                (
-//                    samplePts,
-//                    pointField(samplePts.size(), locationInMesh()),
-//                    info
-//                );
-//
-//                forAll(info, ptI)
-//                {
-//                    if (info[ptI].size() % 2 == 0)
-//                    {
-//                        surfaceVolumeTests[s][ptI] = volumeType::INSIDE;
-//                    }
-//                    else
-//                    {
-//                        surfaceVolumeTests[s][ptI] = volumeType::OUTSIDE;
-//                    }
-//                }
-//            }
-//            else
-            {
             surface.getVolumeType(samplePts, surfaceVolumeTests[s]);
         }
     }
-    }
 
     // Compare the volumeType result for each point wrt to each surface with the
     // reference value and if the points are inside the surface by a given
     // distanceSquared
 
     // Assume that the point is wellInside until demonstrated otherwise.
-    Field<bool> insidePoint(samplePts.size(), true);
+    Field<bool> insideOutsidePoint(samplePts.size(), testForInside);
 
     //Check if the points are inside the surface by the given distance squared
 
@@ -682,7 +656,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
         {
             // If the point is within range of the surface, then it can't be
             // well (in|out)side
-            insidePoint[i] = false;
+            insideOutsidePoint[i] = false;
 
             continue;
         }
@@ -691,31 +665,12 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
         {
             const label regionI = regionOffset_[s];
 
-//            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
-
-//            if
-//            (
-//                !surface.hasVolumeType()
-//             //&& !surface.bounds().contains(samplePts[i])
-//            )
-//            {
-//                continue;
-//            }
-
-            // If one of the pattern tests is failed, then the point cannot be
-            // inside, therefore, if this is a testForInside = true call, the
-            // result is false.  If this is a testForInside = false call, then
-            // the result is true.
             if (normalVolumeTypes_[regionI] == extendedFeatureEdgeMesh::BOTH)
             {
                 continue;
             }
 
-//            Info<< surface.name() << " = "
-//                << volumeType::names[surfaceVolumeTests[s][i]] << endl;
-
             if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
-//            if (surfaceVolumeTests[s][i] != volumeType::INSIDE)
             {
                 if
                 (
@@ -723,7 +678,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
                  == extendedFeatureEdgeMesh::INSIDE
                 )
                 {
-                    insidePoint[i] = false;
+                    insideOutsidePoint[i] = !testForInside;
                     break;
                 }
             }
@@ -735,16 +690,62 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
                  == extendedFeatureEdgeMesh::OUTSIDE
                 )
                 {
-                    insidePoint[i] = false;
+                    insideOutsidePoint[i] = !testForInside;
                     break;
                 }
             }
+//            else
+//            {
+//                // Surface volume type is unknown
+//                Info<< "UNKNOWN" << endl;
+//                // Get nearest face normal
+//
+//                pointField sample(1, samplePts[i]);
+//                scalarField nearestDistSqr(1, GREAT);
+//                List<pointIndexHit> info;
+//                vectorField norms(1);
+//
+//                surface.findNearest(sample, nearestDistSqr, info);
+//                surface.getNormal(info, norms);
+//
+//                vector fN = norms[0];
+//                fN /= mag(fN);
+//
+//                vector hitDir = info[0].rawPoint() - samplePts[i];
+//                hitDir /= mag(hitDir);
+//
+//                if ((fN & hitDir) < 0)
+//                {
+//                    // Point is OUTSIDE
+//
+//                    if
+//                    (
+//                        normalVolumeTypes_[regionI]
+//                     == extendedFeatureEdgeMesh::OUTSIDE
+//                    )
+//                    {
+//                    }
+//                    else
+//                    {
+//                        insidePoint[i] = false;
+//                        break;
+//                    }
+//                }
+//            }
         }
     }
 
-    return insidePoint;
+    return insideOutsidePoint;
+}
 
-    //return wellInOutSide(samplePts, testDistSqr, true);
+
+Foam::Field<bool> Foam::conformationSurfaces::wellInside
+(
+    const pointField& samplePts,
+    const scalarField& testDistSqr
+) const
+{
+    return wellInOutSide(samplePts, testDistSqr, true);
 }
 
 
@@ -764,122 +765,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellOutside
     const scalarField& testDistSqr
 ) const
 {
-    // First check if it is outside any closed surface
+    return wellInOutSide(samplePts, testDistSqr, false);
-
-    List<List<volumeType> > surfaceVolumeTests
-    (
-        surfaces_.size(),
-        List<volumeType>
-        (
-            samplePts.size(),
-            volumeType::UNKNOWN
-        )
-    );
-
-    // Get lists for the volumeTypes for each sample wrt each surface
-    forAll(surfaces_, s)
-    {
-        const searchableSurface& surface(allGeometry_[surfaces_[s]]);
-
-        const label regionI = regionOffset_[s];
-
-        if (normalVolumeTypes_[regionI] != extendedFeatureEdgeMesh::BOTH)
-        {
-            surface.getVolumeType(samplePts, surfaceVolumeTests[s]);
-        }
-    }
-
-    // Compare the volumeType result for each point wrt to each surface with the
-    // reference value and if the points are inside the surface by a given
-    // distanceSquared
-
-    // Assume that the point is wellInside until demonstrated otherwise.
-    Field<bool> outsidePoint(samplePts.size(), false);
-
-    //Check if the points are inside the surface by the given distance squared
-
-    labelList hitSurfaces;
-    List<pointIndexHit> hitInfo;
-    searchableSurfacesQueries::findNearest
-    (
-        allGeometry_,
-        surfaces_,
-        samplePts,
-        testDistSqr,
-        hitSurfaces,
-        hitInfo
-    );
-
-    forAll(samplePts, i)
-    {
-        const pointIndexHit& pHit = hitInfo[i];
-
-        if (pHit.hit())
-        {
-            // If the point is within range of the surface, then it can't be
-            // well (in|out)side
-            outsidePoint[i] = false;
-            //continue;
-        }
-
-        forAll(surfaces_, s)
-        {
-            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
-
-//            if
-//            (
-//                !surface.hasVolumeType()
-//             //&& !surface.bounds().contains(samplePts[i])
-//            )
-//            {
-//                continue;
-//            }
-
-            const label regionI = regionOffset_[s];
-
-//            Info<< s << " " << surfaces_[s] << " " << surface.name() << " "
-//                << normalVolumeTypes_[regionI] << " "
-//                << surfaceVolumeTests[s][i] << endl;
-
-            // If one of the pattern tests is failed, then the point cannot be
-            // inside, therefore, if this is a testForInside = true call, the
-            // result is false.  If this is a testForInside = false call, then
-            // the result is true.
-            if (normalVolumeTypes_[regionI] == extendedFeatureEdgeMesh::BOTH)
-            {
-                continue;
-            }
-
-            if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
-            {
-                if
-                (
-                    normalVolumeTypes_[regionI]
-                 == extendedFeatureEdgeMesh::INSIDE
-                )
-                {
-                    outsidePoint[i] = true;
-                    break;
-                }
-            }
-            else if (surfaceVolumeTests[s][i] == volumeType::INSIDE)
-            {
-                if
-                (
-                    normalVolumeTypes_[regionI]
-                 == extendedFeatureEdgeMesh::OUTSIDE
-                )
-                {
-                    outsidePoint[i] = true;
-                    break;
-                }
-            }
-        }
-    }
-
-    return outsidePoint;
-
-    //return wellInOutSide(samplePts, testDistSqr, false);
 }
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H

@@ -198,12 +198,12 @@ public:
             //- Check if point is closer to the surfaces to conform to than
             //  testDistSqr, in which case return false, otherwise assess in or
             //  outside and return a result depending on the testForInside flag
-//            Field<bool> wellInOutSide
+            Field<bool> wellInOutSide
-//            (
+            (
-//                const pointField& samplePts,
+                const pointField& samplePts,
-//                const scalarField& testDistSqr,
+                const scalarField& testDistSqr,
-//                bool testForInside
+                bool testForInside
-//            ) const;
+            ) const;
 
             //- Check if point is inside surfaces to conform to by at least
             //  testDistSqr

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C

@@ -49,7 +49,7 @@ void Foam::autoDensity::writeOBJ
     fileName name
 ) const
 {
-    OFstream str(foamyHexMesh_.time().path()/name + ".obj");
+    OFstream str(time().path()/name + ".obj");
 
     Pout<< "Writing " << str.name() << endl;
 
@@ -74,11 +74,11 @@ bool Foam::autoDensity::combinedOverlaps(const treeBoundBox& box) const
     if (Pstream::parRun())
     {
         return
-            foamyHexMesh_.decomposition().overlapsThisProcessor(box)
+            decomposition().overlapsThisProcessor(box)
-         || foamyHexMesh_.geometryToConformTo().overlaps(box);
+         || geometryToConformTo().overlaps(box);
     }
 
-    return foamyHexMesh_.geometryToConformTo().overlaps(box);
+    return geometryToConformTo().overlaps(box);
 }
 
 
@@ -87,11 +87,11 @@ bool Foam::autoDensity::combinedInside(const point& p) const
     if (Pstream::parRun())
     {
         return
-            foamyHexMesh_.decomposition().positionOnThisProcessor(p)
+            decomposition().positionOnThisProcessor(p)
-         && foamyHexMesh_.geometryToConformTo().inside(p);
+         && geometryToConformTo().inside(p);
     }
 
-    return foamyHexMesh_.geometryToConformTo().inside(p);
+    return geometryToConformTo().inside(p);
 }
 
 
@@ -103,7 +103,7 @@ Foam::Field<bool> Foam::autoDensity::combinedWellInside
 {
     if (!Pstream::parRun())
     {
-        return foamyHexMesh_.geometryToConformTo().wellInside
+        return geometryToConformTo().wellInside
         (
             pts,
             minimumSurfaceDistanceCoeffSqr_*sqr(sizes)
@@ -117,7 +117,7 @@ Foam::Field<bool> Foam::autoDensity::combinedWellInside
 
     Field<bool> insideA
     (
-        foamyHexMesh_.geometryToConformTo().wellInside
+        geometryToConformTo().wellInside
         (
             pts,
             minimumSurfaceDistanceCoeffSqr_*sqr(sizes)
@@ -126,7 +126,7 @@ Foam::Field<bool> Foam::autoDensity::combinedWellInside
 
     Field<bool> insideB
     (
-        foamyHexMesh_.decomposition().positionOnThisProcessor(pts)
+        decomposition().positionOnThisProcessor(pts)
     );
 
     // inside = insideA && insideB;
@@ -162,14 +162,14 @@ bool Foam::autoDensity::combinedWellInside
 
     if (Pstream::parRun())
     {
-        inside = foamyHexMesh_.decomposition().positionOnThisProcessor(p);
+        inside = decomposition().positionOnThisProcessor(p);
     }
 
     // Perform AND operation between testing the surfaces and the previous
     // result, i.e the parallel result, or in serial, with true.
     inside =
         inside
-     && foamyHexMesh_.geometryToConformTo().wellInside
+     && geometryToConformTo().wellInside
         (
             p,
             minimumSurfaceDistanceCoeffSqr_*sqr(size)
@@ -179,7 +179,7 @@ bool Foam::autoDensity::combinedWellInside
 }
 
 
-void Foam::autoDensity::recurseAndFill
+Foam::label Foam::autoDensity::recurseAndFill
 (
     DynamicList<Vb::Point>& initialPoints,
     const treeBoundBox& bb,
@@ -187,24 +187,31 @@ void Foam::autoDensity::recurseAndFill
     word recursionName
 ) const
 {
+    label maxDepth = 0;
+
     for (direction i = 0; i < 8; i++)
     {
         treeBoundBox subBB = bb.subBbox(i);
 
         word newName = recursionName + "_" + Foam::name(i);
 
-        conformalVoronoiMesh::timeCheck(foamyHexMesh_.time(), newName, debug);
+        conformalVoronoiMesh::timeCheck(time(), newName, debug);
 
         if (combinedOverlaps(subBB))
         {
             if (levelLimit > 0)
             {
+                maxDepth =
+                    max
+                    (
+                        maxDepth,
                         recurseAndFill
                         (
                             initialPoints,
                             subBB,
                             levelLimit - 1,
                             newName
+                        )
                     );
             }
             else
@@ -222,12 +229,17 @@ void Foam::autoDensity::recurseAndFill
 
                 if (!fillBox(initialPoints, subBB, true))
                 {
+                    maxDepth =
+                        max
+                        (
+                            maxDepth,
                             recurseAndFill
                             (
                                 initialPoints,
                                 subBB,
                                 levelLimit - 1,
                                 newName
+                            )
                         );
                 }
             }
@@ -247,12 +259,17 @@ void Foam::autoDensity::recurseAndFill
 
             if (!fillBox(initialPoints, subBB, false))
             {
+                maxDepth =
+                    max
+                    (
+                        maxDepth,
                         recurseAndFill
                         (
                             initialPoints,
                             subBB,
                             levelLimit - 1,
                             newName
+                        )
                     );
             }
         }
@@ -268,6 +285,8 @@ void Foam::autoDensity::recurseAndFill
             }
         }
     }
+
+    return maxDepth + 1;
 }
 
 
@@ -278,9 +297,7 @@ bool Foam::autoDensity::fillBox
     bool overlapping
 ) const
 {
-    const conformationSurfaces& geometry(foamyHexMesh_.geometryToConformTo());
+    const conformationSurfaces& geometry = geometryToConformTo();
-
-    Random& rnd = foamyHexMesh_.rndGen();
 
     unsigned int initialSize = initialPoints.size();
 
@@ -333,15 +350,14 @@ bool Foam::autoDensity::fillBox
 
     if (!overlapping && !wellInside)
     {
-        // If this is an inside box then then it is possible to fill points very
+        // If this is an inside box then it is possible to fill points very
         // close to the boundary, to prevent this, check the corners and sides
         // of the box so ensure that they are "wellInside".  If not, set as an
         // overlapping box.
 
         pointField corners(bb.points());
 
-        scalarField cornerSizes =
+        scalarField cornerSizes = cellShapeControls().cellSize(corners);
-            foamyHexMesh_.cellShapeControls().cellSize(corners);
 
         Field<bool> insideCorners = combinedWellInside(corners, cornerSizes);
 
@@ -450,8 +466,7 @@ bool Foam::autoDensity::fillBox
                         );
                 }
 
-                lineSizes =
+                lineSizes = cellShapeControls().cellSize(linePoints);
-                    foamyHexMesh_.cellShapeControls().cellSize(linePoints);
 
                 Field<bool> insideLines = combinedWellInside
                 (
@@ -537,9 +552,12 @@ bool Foam::autoDensity::fillBox
                         min
                       + vector
                         (
-                            delta.x()*(i + 0.5 + 0.1*(rnd.scalar01() - 0.5)),
+                            delta.x()
-                            delta.y()*(j + 0.5 + 0.1*(rnd.scalar01() - 0.5)),
+                           *(i + 0.5 + 0.1*(rndGen().scalar01() - 0.5)),
-                            delta.z()*(k + 0.5 + 0.1*(rnd.scalar01() - 0.5))
+                            delta.y()
+                           *(j + 0.5 + 0.1*(rndGen().scalar01() - 0.5)),
+                            delta.z()
+                           *(k + 0.5 + 0.1*(rndGen().scalar01() - 0.5))
                         );
                 }
             }
@@ -550,10 +568,7 @@ bool Foam::autoDensity::fillBox
         // corner when only some these points are required.
         shuffle(samplePoints);
 
-        scalarField sampleSizes = foamyHexMesh_.cellShapeControls().cellSize
+        scalarField sampleSizes = cellShapeControls().cellSize(samplePoints);
-        (
-            samplePoints
-        );
 
         Field<bool> insidePoints = combinedWellInside
         (
@@ -632,7 +647,7 @@ bool Foam::autoDensity::fillBox
             {
                 trialPoints++;
 
-                point p = samplePoints[i];
+                const point& p = samplePoints[i];
 
                 scalar localSize = sampleSizes[i];
 
@@ -647,7 +662,7 @@ bool Foam::autoDensity::fillBox
                 // TODO - is there a lot of cost in the 1/density calc?  Could
                 // assess on
                 //    (1/maxDensity)/(1/localDensity) = minVolume/localVolume
-                if (localDensity/maxDensity > rnd.scalar01())
+                if (localDensity/maxDensity > rndGen().scalar01())
                 {
                     scalar localVolume = 1/localDensity;
 
@@ -660,7 +675,7 @@ bool Foam::autoDensity::fillBox
                         scalar addProbability =
                            (totalVolume - volumeAdded)/localVolume;
 
-                        scalar r = rnd.scalar01();
+                        scalar r = rndGen().scalar01();
 
                         if (debug)
                         {
@@ -714,9 +729,9 @@ bool Foam::autoDensity::fillBox
         {
             trialPoints++;
 
-            point p = min + cmptMultiply(span, rnd.vector01());
+            point p = min + cmptMultiply(span, rndGen().vector01());
 
-            scalar localSize = foamyHexMesh_.cellShapeControls().cellSize(p);
+            scalar localSize = cellShapeControls().cellSize(p);
 
             bool insidePoint = false;
 
@@ -770,7 +785,7 @@ bool Foam::autoDensity::fillBox
 
                 // Accept possible placements proportional to the relative local
                 // density
-                if (localDensity/maxDensity > rnd.scalar01())
+                if (localDensity/maxDensity > rndGen().scalar01())
                 {
                     scalar localVolume = 1/localDensity;
 
@@ -783,7 +798,7 @@ bool Foam::autoDensity::fillBox
                         scalar addProbability =
                             (totalVolume - volumeAdded)/localVolume;
 
-                        scalar r = rnd.scalar01();
+                        scalar r = rndGen().scalar01();
 
                         if (debug)
                         {
@@ -847,10 +862,23 @@ bool Foam::autoDensity::fillBox
 autoDensity::autoDensity
 (
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 :
-    initialPointsMethod(typeName, initialPointsDict, foamyHexMesh),
+    initialPointsMethod
+    (
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
     globalTrialPoints_(0),
     minCellSizeLimit_
     (
@@ -875,8 +903,7 @@ autoDensity::autoDensity
                 "const dictionary& initialPointsDict,"
                 "const conformalVoronoiMesh& foamyHexMesh"
             ")"
-        )
+        )   << "The maxSizeRatio must be greater than one to be sensible, "
-            << "The maxSizeRatio must be greater than one to be sensible, "
             << "setting to " << maxSizeRatio_
             << endl;
     }
@@ -893,14 +920,14 @@ List<Vb::Point> autoDensity::initialPoints() const
     // on whether this is a parallel run.
     if (Pstream::parRun())
     {
-        hierBB = foamyHexMesh_.decomposition().procBounds();
+        hierBB = decomposition().procBounds();
     }
     else
     {
         // Extend the global box to move it off large plane surfaces
-        hierBB = foamyHexMesh_.geometryToConformTo().globalBounds().extend
+        hierBB = geometryToConformTo().globalBounds().extend
         (
-            foamyHexMesh_.rndGen(),
+            rndGen(),
             1e-6
         );
     }
@@ -914,7 +941,7 @@ List<Vb::Point> autoDensity::initialPoints() const
         Pout<< "    Filling box " << hierBB << endl;
     }
 
-    recurseAndFill
+    label maxDepth = recurseAndFill
     (
         initialPoints,
         hierBB,
@@ -932,11 +959,16 @@ List<Vb::Point> autoDensity::initialPoints() const
         reduce(globalTrialPoints_, sumOp<label>());
     }
 
-    Info<< "        " << nInitialPoints << " points placed" << nl
+    Info<< incrIndent << incrIndent
-        << "        " << globalTrialPoints_ << " locations queried" << nl
+        << indent << nInitialPoints << " points placed" << nl
-        << "        "
+        << indent << globalTrialPoints_ << " locations queried" << nl
+        << indent
         << scalar(nInitialPoints)/scalar(max(globalTrialPoints_, 1))
-        << " success rate"
+        << " success rate" << nl
+        << indent
+        << returnReduce(maxDepth, maxOp<label>())
+        << " levels of recursion (maximum)"
+        << decrIndent << decrIndent
         << endl;
 
     return initialPoints;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.H

@@ -114,7 +114,7 @@ private:
         ) const;
 
         //- Descend into octants of the supplied bound box
-        void recurseAndFill
+        label recurseAndFill
         (
             DynamicList<Vb::Point>& initialPoints,
             const treeBoundBox& bb,
@@ -144,7 +144,11 @@ public:
         autoDensity
         (
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C

@@ -41,10 +41,23 @@ addToRunTimeSelectionTable(initialPointsMethod, bodyCentredCubic, dictionary);
 bodyCentredCubic::bodyCentredCubic
 (
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 :
-    initialPointsMethod(typeName, initialPointsDict, foamyHexMesh),
+    initialPointsMethod
+    (
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
     initialCellSize_(readScalar(detailsDict().lookup("initialCellSize"))),
     randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
     randomPerturbationCoeff_
@@ -64,11 +77,11 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
     // on whether this is a parallel run.
     if (Pstream::parRun())
     {
-        bb = foamyHexMesh_.decomposition().procBounds();
+        bb = decomposition().procBounds();
     }
     else
     {
-        bb = foamyHexMesh_.geometryToConformTo().globalBounds();
+        bb = geometryToConformTo().globalBounds();
     }
 
     scalar x0 = bb.min().x();
@@ -87,8 +100,6 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
 
     delta *= pow((1.0/2.0),-(1.0/3.0));
 
-    Random& rndGen = foamyHexMesh_.rndGen();
-
     scalar pert = randomPerturbationCoeff_*cmptMin(delta);
 
     DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
@@ -118,17 +129,14 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    pA.x() += pert*(rndGen.scalar01() - 0.5);
+                    pA.x() += pert*(rndGen().scalar01() - 0.5);
-                    pA.y() += pert*(rndGen.scalar01() - 0.5);
+                    pA.y() += pert*(rndGen().scalar01() - 0.5);
-                    pA.z() += pert*(rndGen.scalar01() - 0.5);
+                    pA.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
-                const backgroundMeshDecomposition& decomp =
-                    foamyHexMesh_.decomposition();
-
                 if (Pstream::parRun())
                 {
-                    if (decomp.positionOnThisProcessor(pA))
+                    if (decomposition().positionOnThisProcessor(pA))
                     {
                         // Add this point in parallel only if this position is
                         // on this processor.
@@ -142,14 +150,14 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    pB.x() += pert*(rndGen.scalar01() - 0.5);
+                    pB.x() += pert*(rndGen().scalar01() - 0.5);
-                    pB.y() += pert*(rndGen.scalar01() - 0.5);
+                    pB.y() += pert*(rndGen().scalar01() - 0.5);
-                    pB.z() += pert*(rndGen.scalar01() - 0.5);
+                    pB.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
                 if (Pstream::parRun())
                 {
-                    if (decomp.positionOnThisProcessor(pB))
+                    if (decomposition().positionOnThisProcessor(pB))
                     {
                         // Add this point in parallel only if this position is
                         // on this processor.
@@ -165,14 +173,11 @@ List<Vb::Point> bodyCentredCubic::initialPoints() const
             points.setSize(pI);
 
             Field<bool> insidePoints =
-                foamyHexMesh_.geometryToConformTo().wellInside
+                geometryToConformTo().wellInside
                 (
                     points,
                     minimumSurfaceDistanceCoeffSqr_
-                   *sqr
+                   *sqr(cellShapeControls().cellSize(points))
-                    (
-                        foamyHexMesh_.cellShapeControls().cellSize(points)
-                    )
                 );
 
             forAll(insidePoints, i)

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.H

@@ -77,7 +77,11 @@ public:
         bodyCentredCubic
         (
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C

@@ -41,10 +41,23 @@ addToRunTimeSelectionTable(initialPointsMethod, faceCentredCubic, dictionary);
 faceCentredCubic::faceCentredCubic
 (
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 :
-    initialPointsMethod(typeName, initialPointsDict, foamyHexMesh),
+    initialPointsMethod
+    (
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
     initialCellSize_(readScalar(detailsDict().lookup("initialCellSize"))),
     randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
     randomPerturbationCoeff_
@@ -64,11 +77,11 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
     // on whether this is a parallel run.
     if (Pstream::parRun())
     {
-        bb = foamyHexMesh_.decomposition().procBounds();
+        bb = decomposition().procBounds();
     }
     else
     {
-        bb = foamyHexMesh_.geometryToConformTo().globalBounds();
+        bb = geometryToConformTo().globalBounds();
     }
 
     scalar x0 = bb.min().x();
@@ -87,8 +100,6 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
     delta *= pow((1.0/4.0),-(1.0/3.0));
 
-    Random& rndGen = foamyHexMesh_.rndGen();
-
     scalar pert = randomPerturbationCoeff_*cmptMin(delta);
 
     DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
@@ -116,17 +127,14 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen.scalar01() - 0.5);
+                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen.scalar01() - 0.5);
+                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen.scalar01() - 0.5);
+                    p.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
-                const backgroundMeshDecomposition& decomp =
-                    foamyHexMesh_.decomposition();
-
                 if (Pstream::parRun())
                 {
-                    if (decomp.positionOnThisProcessor(p))
+                    if (decomposition().positionOnThisProcessor(p))
                     {
                         // Add this point in parallel only if this position is
                         // on this processor.
@@ -147,14 +155,14 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen.scalar01() - 0.5);
+                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen.scalar01() - 0.5);
+                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen.scalar01() - 0.5);
+                    p.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
                 if (Pstream::parRun())
                 {
-                    if (decomp.positionOnThisProcessor(p))
+                    if (decomposition().positionOnThisProcessor(p))
                     {
                         // Add this point in parallel only if this position is
                         // on this processor.
@@ -175,14 +183,14 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen.scalar01() - 0.5);
+                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen.scalar01() - 0.5);
+                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen.scalar01() - 0.5);
+                    p.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
                 if (Pstream::parRun())
                 {
-                    if (decomp.positionOnThisProcessor(p))
+                    if (decomposition().positionOnThisProcessor(p))
                     {
                         // Add this point in parallel only if this position is
                         // on this processor.
@@ -203,14 +211,14 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen.scalar01() - 0.5);
+                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen.scalar01() - 0.5);
+                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen.scalar01() - 0.5);
+                    p.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
                 if (Pstream::parRun())
                 {
-                    if (decomp.positionOnThisProcessor(p))
+                    if (decomposition().positionOnThisProcessor(p))
                     {
                         // Add this point in parallel only if this position is
                         // on this processor.
@@ -226,14 +234,11 @@ List<Vb::Point> faceCentredCubic::initialPoints() const
             points.setSize(pI);
 
             Field<bool> insidePoints =
-                foamyHexMesh_.geometryToConformTo().wellInside
+                geometryToConformTo().wellInside
                 (
                     points,
                     minimumSurfaceDistanceCoeffSqr_
-                   *sqr
+                   *sqr(cellShapeControls().cellSize(points))
-                    (
-                        foamyHexMesh_.cellShapeControls().cellSize(points)
-                    )
                 );
 
             forAll(insidePoints, i)

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.H

@@ -77,7 +77,11 @@ public:
         faceCentredCubic
         (
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.C

@@ -43,11 +43,19 @@ initialPointsMethod::initialPointsMethod
 (
     const word& type,
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 :
     dictionary(initialPointsDict),
-    foamyHexMesh_(foamyHexMesh),
+    runTime_(runTime),
+    rndGen_(rndGen),
+    geometryToConformTo_(geometryToConformTo),
+    cellShapeControls_(cellShapeControls),
+    decomposition_(decomposition),
     detailsDict_(subDict(type + "Coeffs")),
     minimumSurfaceDistanceCoeffSqr_
     (
@@ -68,7 +76,11 @@ initialPointsMethod::initialPointsMethod
 autoPtr<initialPointsMethod> initialPointsMethod::New
 (
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 {
     word initialPointsMethodTypeName
@@ -99,7 +111,15 @@ autoPtr<initialPointsMethod> initialPointsMethod::New
     return
         autoPtr<initialPointsMethod>
         (
-            cstrIter()(initialPointsDict, foamyHexMesh)
+            cstrIter()
+            (
+                initialPointsDict,
+                runTime,
+                rndGen,
+                geometryToConformTo,
+                cellShapeControls,
+                decomposition
+            )
         );
 }
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H

@@ -62,8 +62,15 @@ protected:
 
     // Protected data
 
-        //- Reference to the conformalVoronoiMesh holding this object
+        const Time& runTime_;
-        const conformalVoronoiMesh& foamyHexMesh_;
+
+        Random& rndGen_;
+
+        const conformationSurfaces& geometryToConformTo_;
+
+        const cellShapeControl& cellShapeControls_;
+
+        const autoPtr<backgroundMeshDecomposition>& decomposition_;
 
         //- Method details dictionary
         dictionary detailsDict_;
@@ -102,9 +109,20 @@ public:
             dictionary,
             (
                 const dictionary& initialPointsDict,
-                const conformalVoronoiMesh& foamyHexMesh
+                const Time& runTime,
+                Random& rndGen,
+                const conformationSurfaces& geometryToConformTo,
+                const cellShapeControl& cellShapeControls,
+                const autoPtr<backgroundMeshDecomposition>& decomposition
             ),
-            (initialPointsDict, foamyHexMesh)
+            (
+                initialPointsDict,
+                runTime,
+                rndGen,
+                geometryToConformTo,
+                cellShapeControls,
+                decomposition
+            )
         );
 
 
@@ -115,7 +133,11 @@ public:
         (
             const word& type,
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 
@@ -125,7 +147,11 @@ public:
         static autoPtr<initialPointsMethod> New
         (
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 
@@ -135,6 +161,33 @@ public:
 
     // Member Functions
 
+        // Access
+
+            const Time& time() const
+            {
+                return runTime_;
+            }
+
+            Random& rndGen() const
+            {
+                return rndGen_;
+            }
+
+            const conformationSurfaces& geometryToConformTo() const
+            {
+                return geometryToConformTo_;
+            }
+
+            const cellShapeControl& cellShapeControls() const
+            {
+                return cellShapeControls_;
+            }
+
+            const backgroundMeshDecomposition& decomposition() const
+            {
+                return decomposition_;
+            }
+
             //- Const access to the details dictionary
             const dictionary& detailsDict() const
             {
@@ -146,6 +199,9 @@ public:
                 return fixInitialPoints_;
             }
 
+
+        // Queries
+
             //- Return the initial points for the conformalVoronoiMesh
             virtual List<Vb::Point> initialPoints() const = 0;
 };

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C

@@ -41,12 +41,34 @@ addToRunTimeSelectionTable(initialPointsMethod, pointFile, dictionary);
 pointFile::pointFile
 (
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 :
-    initialPointsMethod(typeName, initialPointsDict, foamyHexMesh),
+    initialPointsMethod
-    pointFileName_(detailsDict().lookup("pointFile"))
+    (
-{}
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
+    pointFileName_(detailsDict().lookup("pointFile")),
+    insideOutsideCheck_(detailsDict().lookup("insideOutsideCheck")),
+    randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
+    randomPerturbationCoeff_
+    (
+        readScalar(detailsDict().lookup("randomPerturbationCoeff"))
+    )
+{
+    Info<< "    Inside/Outside check is " << insideOutsideCheck_.asText()
+        << endl;
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -58,8 +80,8 @@ List<Vb::Point> pointFile::initialPoints() const
         IOobject
         (
             pointFileName_.name(),
-            foamyHexMesh_.time().timeName(),
+            time().timeName(),
-            foamyHexMesh_.time(),
+            time(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
@@ -85,13 +107,13 @@ List<Vb::Point> pointFile::initialPoints() const
 
         if (!isParentFile)
         {
-            foamyHexMesh_.decomposition().distributePoints(points);
+            decomposition().distributePoints(points);
         }
         else
         {
             // Otherwise, this is assumed to be points covering the whole
             // domain, so filter the points to be only those on this processor
-            boolList procPt(foamyHexMesh_.positionOnThisProc(points));
+            boolList procPt(decomposition().positionOnThisProcessor(points));
 
             List<boolList> allProcPt(Pstream::nProcs());
 
@@ -126,15 +148,17 @@ List<Vb::Point> pointFile::initialPoints() const
         }
     }
 
-    Field<bool> insidePoints = foamyHexMesh_.geometryToConformTo().wellInside
+    Field<bool> insidePoints(points.size(), true);
+
+    if (insideOutsideCheck_)
+    {
+        insidePoints = geometryToConformTo().wellInside
         (
             points,
             minimumSurfaceDistanceCoeffSqr_
-       *sqr
+           *sqr(cellShapeControls().cellSize(points))
-        (
-            foamyHexMesh_.cellShapeControls().cellSize(points)
-        )
         );
+    }
 
     DynamicList<Vb::Point> initialPoints(insidePoints.size()/10);
 
@@ -142,7 +166,14 @@ List<Vb::Point> pointFile::initialPoints() const
     {
         if (insidePoints[i])
         {
-            const point& p(points[i]);
+            point& p = points[i];
+
+            if (randomiseInitialGrid_)
+            {
+                p.x() += randomPerturbationCoeff_*(rndGen().scalar01() - 0.5);
+                p.y() += randomPerturbationCoeff_*(rndGen().scalar01() - 0.5);
+                p.z() += randomPerturbationCoeff_*(rndGen().scalar01() - 0.5);
+            }
 
             initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
         }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "fileName.H"
 #include "pointIOField.H"
+#include "Switch.H"
 #include "initialPointsMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,6 +62,15 @@ private:
         //- The initial cell spacing
         fileName pointFileName_;
 
+        //- Check if inserted points are inside or outside
+        Switch insideOutsideCheck_;
+
+        //- Should the initial positions be randomised
+        Switch randomiseInitialGrid_;
+
+        //- Randomise the initial positions by fraction of the initialCellSize_
+        scalar randomPerturbationCoeff_;
+
 
 public:
 
@@ -73,7 +83,11 @@ public:
         pointFile
         (
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C

@@ -41,10 +41,23 @@ addToRunTimeSelectionTable(initialPointsMethod, uniformGrid, dictionary);
 uniformGrid::uniformGrid
 (
     const dictionary& initialPointsDict,
-    const conformalVoronoiMesh& foamyHexMesh
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
 )
 :
-    initialPointsMethod(typeName, initialPointsDict, foamyHexMesh),
+    initialPointsMethod
+    (
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
     initialCellSize_(readScalar(detailsDict().lookup("initialCellSize"))),
     randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
     randomPerturbationCoeff_
@@ -64,11 +77,11 @@ List<Vb::Point> uniformGrid::initialPoints() const
     // on whether this is a parallel run.
     if (Pstream::parRun())
     {
-        bb = foamyHexMesh_.decomposition().procBounds();
+        bb = decomposition().procBounds();
     }
     else
     {
-        bb = foamyHexMesh_.geometryToConformTo().globalBounds();
+        bb = geometryToConformTo().globalBounds();
     }
 
     scalar x0 = bb.min().x();
@@ -87,8 +100,6 @@ List<Vb::Point> uniformGrid::initialPoints() const
 
     delta *= pow((1.0),-(1.0/3.0));
 
-    Random& rndGen = foamyHexMesh_.rndGen();
-
     scalar pert = randomPerturbationCoeff_*cmptMin(delta);
 
     // Initialise points list
@@ -117,15 +128,15 @@ List<Vb::Point> uniformGrid::initialPoints() const
 
                 if (randomiseInitialGrid_)
                 {
-                    p.x() += pert*(rndGen.scalar01() - 0.5);
+                    p.x() += pert*(rndGen().scalar01() - 0.5);
-                    p.y() += pert*(rndGen.scalar01() - 0.5);
+                    p.y() += pert*(rndGen().scalar01() - 0.5);
-                    p.z() += pert*(rndGen.scalar01() - 0.5);
+                    p.z() += pert*(rndGen().scalar01() - 0.5);
                 }
 
                 if
                 (
                     Pstream::parRun()
-                 && !foamyHexMesh_.decomposition().positionOnThisProcessor(p)
+                 && !decomposition().positionOnThisProcessor(p)
                 )
                 {
                     // Skip this point if, in parallel, this position is not on
@@ -139,13 +150,13 @@ List<Vb::Point> uniformGrid::initialPoints() const
             points.setSize(pI);
 
             Field<bool> insidePoints =
-                foamyHexMesh_.geometryToConformTo().wellInside
+                geometryToConformTo().wellInside
                 (
                     points,
                     minimumSurfaceDistanceCoeffSqr_
                    *sqr
                     (
-                        foamyHexMesh_.cellShapeControls().cellSize(points)
+                        cellShapeControls().cellSize(points)
                     )
                 );
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.H

@@ -77,7 +77,11 @@ public:
         uniformGrid
         (
             const dictionary& initialPointsDict,
-            const conformalVoronoiMesh& foamyHexMesh
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
         );
 
 

applications/utilities/mesh/generation/foamyHexMesh/foamyHexMesh.C

@@ -30,6 +30,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
+#include "Time.H"
+#include "IOdictionary.H"
+#include "searchableSurfaces.H"
 #include "conformalVoronoiMesh.H"
 #include "vtkSetWriter.H"
 
@@ -137,9 +140,10 @@ int main(int argc, char *argv[])
 
             mesh.move();
 
-            Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            Info<< nl
+                << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                 << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-                << endl;
+                << nl << endl;
         }
     }
 

applications/utilities/mesh/generation/foamyHexMesh/vectorTools/vectorTools.H

@@ -27,8 +27,6 @@ Class
 Description
     Functions for analysing the relationships between vectors
 
-SourceFiles
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef vectorTools_H
@@ -49,7 +47,7 @@ namespace Foam
 //- Collection of functions for testing relationships between two vectors.
 namespace vectorTools
 {
-    //- Test if a and b are parallel: a.b = 1
+    //- Test if a and b are parallel: a^b = 0
     //  Uses the cross product, so the tolerance is proportional to
     //  the sine of the angle between a and b in radians
     template<typename T>
@@ -146,6 +144,7 @@ namespace vectorTools
 
 } // End namespace vectorTools
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -420,11 +420,24 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
 
     const refinementSurfaces& rf = surfacePtr();
 
-    Info<< setw(20) << "Region"
+    // Determine maximum region name length
+    label maxLen = 0;
+    forAll(rf.surfaces(), surfI)
+    {
+        label geomI = rf.surfaces()[surfI];
+        const wordList& regionNames = allGeometry.regionNames()[geomI];
+        forAll(regionNames, regionI)
+        {
+            maxLen = Foam::max(maxLen, label(regionNames[regionI].size()));
+        }
+    }
+
+
+    Info<< setw(maxLen) << "Region"
         << setw(10) << "Min Level"
         << setw(10) << "Max Level"
         << setw(10) << "Gap Level" << nl
-        << setw(20) << "------"
+        << setw(maxLen) << "------"
         << setw(10) << "---------"
         << setw(10) << "---------"
         << setw(10) << "---------" << endl;
@@ -441,7 +454,7 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
         {
             label globalI = rf.globalRegion(surfI, regionI);
 
-            Info<< setw(20) << regionNames[regionI]
+            Info<< setw(maxLen) << regionNames[regionI]
                 << setw(10) << rf.minLevel()[globalI]
                 << setw(10) << rf.maxLevel()[globalI]
                 << setw(10) << rf.gapLevel()[globalI] << endl;
@@ -723,111 +736,117 @@ int main(int argc, char *argv[])
 
     autoPtr<fvMesh> meshPtr;
 
-    if (surfaceSimplify)
+//    if (surfaceSimplify)
-    {
+//    {
-        IOdictionary foamyHexMeshDict
+//        IOdictionary foamyHexMeshDict
-        (
+//        (
-           IOobject
+//           IOobject
-           (
+//           (
-                "foamyHexMeshDict",
+//                "foamyHexMeshDict",
-                runTime.system(),
+//                runTime.system(),
-                runTime,
+//                runTime,
-                IOobject::MUST_READ_IF_MODIFIED,
+//                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE
+//                IOobject::NO_WRITE
-           )
+//           )
-        );
+//        );
-
+//
-        const dictionary& motionDict =
+//        const dictionary& motionDict =
-            foamyHexMeshDict.subDict("motionControl");
+//            foamyHexMeshDict.subDict("motionControl");
-
+//
-        const scalar defaultCellSize =
+//        const scalar defaultCellSize =
-            readScalar(motionDict.lookup("defaultCellSize"));
+//            readScalar(motionDict.lookup("defaultCellSize"));
-
+//
-        Info<< "Constructing single cell mesh from boundBox" << nl << endl;
+//        Info<< "Constructing single cell mesh from boundBox" << nl << endl;
-
+//
-        boundBox bb(args.optionRead<boundBox>("surfaceSimplify"));
+//        boundBox bb(args.optionRead<boundBox>("surfaceSimplify"));
-
+//
-        labelList owner(6, label(0));
+//        labelList owner(6, label(0));
-        labelList neighbour(0);
+//        labelList neighbour(0);
-
+//
-        const cellModel& hexa = *(cellModeller::lookup("hex"));
+//        const cellModel& hexa = *(cellModeller::lookup("hex"));
-        faceList faces = hexa.modelFaces();
+//        faceList faces = hexa.modelFaces();
-
+//
-        meshPtr.set
+//        meshPtr.set
-        (
+//        (
-            new fvMesh
+//            new fvMesh
-            (
+//            (
-                IOobject
+//                IOobject
-                (
+//                (
-                    fvMesh::defaultRegion,
+//                    fvMesh::defaultRegion,
-                    runTime.timeName(),
+//                    runTime.timeName(),
-                    runTime,
+//                    runTime,
-                    IOobject::NO_READ
+//                    IOobject::NO_READ
-                ),
+//                ),
-                xferMove<Field<vector> >(bb.points()()),
+//                xferMove<Field<vector> >(bb.points()()),
-                faces.xfer(),
+//                faces.xfer(),
-                owner.xfer(),
+//                owner.xfer(),
-                neighbour.xfer()
+//                neighbour.xfer()
-            )
+//            )
-        );
+//        );
-
+//
-        List<polyPatch*> patches(1);
+//        List<polyPatch*> patches(1);
-
+//
-        patches[0] = new wallPolyPatch
+//        patches[0] = new wallPolyPatch
-        (
+//        (
-            "boundary",
+//            "boundary",
-            6,
+//            6,
-            0,
+//            0,
-            0,
+//            0,
-            meshPtr().boundaryMesh(),
+//            meshPtr().boundaryMesh(),
-            emptyPolyPatch::typeName
+//            wallPolyPatch::typeName
-        );
+//        );
-
+//
-        meshPtr().addFvPatches(patches);
+//        meshPtr().addFvPatches(patches);
-
+//
-        const scalar initialCellSize = ::pow(meshPtr().V()[0], 1.0/3.0);
+//        const scalar initialCellSize = ::pow(meshPtr().V()[0], 1.0/3.0);
-        const label initialRefLevels =
+//        const label initialRefLevels =
-            ::log(initialCellSize/defaultCellSize)/::log(2);
+//            ::log(initialCellSize/defaultCellSize)/::log(2);
-
+//
-        Info<< "Default cell size = " << defaultCellSize << endl;
+//        Info<< "Default cell size = " << defaultCellSize << endl;
-        Info<< "Initial cell size = " << initialCellSize << endl;
+//        Info<< "Initial cell size = " << initialCellSize << endl;
-
+//
-        Info<< "Initial refinement levels = " << initialRefLevels << endl;
+//        Info<< "Initial refinement levels = " << initialRefLevels << endl;
-
+//
-        Info<< "Mesh starting size = " << meshPtr().nCells() << endl;
+//        Info<< "Mesh starting size = " << meshPtr().nCells() << endl;
-
+//
-        // meshCutter must be destroyed before writing the mesh otherwise it
+//        // meshCutter must be destroyed before writing the mesh otherwise it
-        // writes the cellLevel/pointLevel files
+//        // writes the cellLevel/pointLevel files
-        {
+//        {
-            hexRef8 meshCutter(meshPtr(), false);
+//            hexRef8 meshCutter(meshPtr(), false);
-
+//
-            for (label refineI = 0; refineI < initialRefLevels; ++refineI)
+//            for (label refineI = 0; refineI < initialRefLevels; ++refineI)
-            {
+//            {
-                // Mesh changing engine.
+//                // Mesh changing engine.
-                polyTopoChange meshMod(meshPtr(), true);
+//                polyTopoChange meshMod(meshPtr(), true);
-
+//
-                // Play refinement commands into mesh changer.
+//                // Play refinement commands into mesh changer.
-                meshCutter.setRefinement(identity(meshPtr().nCells()), meshMod);
+//                meshCutter.setRefinement
-
+//                (
-                // Create mesh (no inflation), return map from old to new mesh.
+//                    identity(meshPtr().nCells()),
-                autoPtr<mapPolyMesh> map = meshMod.changeMesh(meshPtr(), false);
+//                    meshMod
-
+//                );
-                // Update fields
+//
-                meshPtr().updateMesh(map);
+//                // Create mesh (no inflation), return map from old to new mesh
-
+//                autoPtr<mapPolyMesh> map =
-                // Delete mesh volumes.
+//                    meshMod.changeMesh(meshPtr(), false);
-                meshPtr().clearOut();
+//
-
+//                // Update fields
-                Info<< "Refinement Iteration " << refineI + 1
+//                meshPtr().updateMesh(map);
-                    << ", Mesh size = " << meshPtr().nCells() << endl;
+//
-            }
+//                // Delete mesh volumes.
-        }
+//                meshPtr().clearOut();
-
+//
-        Info<< "Mesh end size = " << meshPtr().nCells() << endl;
+//                Info<< "Refinement Iteration " << refineI + 1
-
+//                    << ", Mesh size = " << meshPtr().nCells() << endl;
-        meshPtr().write();
+//            }
-    }
+//        }
-    else
+//
+//        Info<< "Mesh end size = " << meshPtr().nCells() << endl;
+//
+//        Info<< "Create mesh" << endl;
+//        meshPtr().write();
+//    }
+//    else
     {
         Foam::Info
             << "Create mesh for time = "
@@ -1521,6 +1540,21 @@ int main(int argc, char *argv[])
             includePatches,
             outFileName
         );
+
+        pointIOField cellCentres
+        (
+            IOobject
+            (
+                "internalCellCentres",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh.cellCentres()
+        );
+
+        cellCentres.write();
     }
 
 

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C

@@ -206,7 +206,18 @@ Foam::functionEntries::codeStream::getFunction
             off_t masterSize = mySize;
             Pstream::scatter(masterSize);
 
+            if (debug)
+            {
+                Pout<< endl<< "on processor " << Pstream::myProcNo()
+                    << " have masterSize:" << masterSize
+                    << " and localSize:" << mySize
+                    << endl;
+            }
+
+
             if (mySize < masterSize)
+            {
+                if (debug)
                 {
                     Pout<< "Local file " << libPath
                         << " not of same size (" << mySize
@@ -214,6 +225,7 @@ Foam::functionEntries::codeStream::getFunction
                         << masterSize << "). Waiting for "
                         << regIOobject::fileModificationSkew
                         << " seconds." << endl;
+                }
                 Foam::sleep(regIOobject::fileModificationSkew);
 
                 // Recheck local size
@@ -237,6 +249,14 @@ Foam::functionEntries::codeStream::getFunction
                         << exit(FatalIOError);
                 }
             }
+
+            if (debug)
+            {
+                Pout<< endl<< "on processor " << Pstream::myProcNo()
+                    << " after waiting: have masterSize:" << masterSize
+                    << " and localSize:" << mySize
+                    << endl;
+            }
         }
 
         if (isA<IOdictionary>(topDict(parentDict)))
@@ -244,6 +264,12 @@ Foam::functionEntries::codeStream::getFunction
             // Cached access to dl libs. Guarantees clean up upon destruction
             // of Time.
             dlLibraryTable& dlLibs = libs(parentDict);
+
+            if (debug)
+            {
+                Pout<< "Opening cached dictionary:" << libPath << endl;
+            }
+
             if (!dlLibs.open(libPath, false))
             {
                 FatalIOErrorIn
@@ -261,10 +287,28 @@ Foam::functionEntries::codeStream::getFunction
         else
         {
             // Uncached opening of libPath
+            if (debug)
+            {
+                Pout<< "Opening uncached dictionary:" << libPath << endl;
+            }
             lib = dlOpen(libPath, true);
         }
     }
 
+    bool haveLib = lib;
+    reduce(haveLib, andOp<bool>());
+
+    if (!haveLib)
+    {
+        FatalIOErrorIn
+        (
+            "functionEntries::codeStream::execute(..)",
+            parentDict
+        )   << "Failed loading library " << libPath
+            << " on some processors."
+            << exit(FatalIOError);
+    }
+
 
     // Find the function handle in the library
     streamingFunctionType function =

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduAddressing.C

@@ -267,7 +267,13 @@ Foam::Tuple2<Foam::label, Foam::scalar> Foam::lduAddressing::band() const
     }
 
     label bandwidth = max(cellBandwidth);
-    scalar profile = sum(1.0*cellBandwidth);
+
+    // Do not use field algebra because of conversion label to scalar
+    scalar profile = 0.0;
+    forAll(cellBandwidth, cellI)
+    {
+        profile += 1.0*cellBandwidth[cellI];
+    }
 
     return Tuple2<label, scalar>(bandwidth, profile);
 }

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C

@@ -41,19 +41,6 @@ defineTypeNameAndDebug(polyBoundaryMesh, 0);
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::labelList Foam::polyBoundaryMesh::ident(const label len)
-{
-    labelList elems(len);
-    forAll(elems, elemI)
-    {
-        elems[elemI] = elemI;
-    }
-    return elems;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::polyBoundaryMesh::polyBoundaryMesh

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H

@@ -78,9 +78,6 @@ class polyBoundaryMesh
 
     // Private Member Functions
 
-        //- Create identity map
-        static labelList ident(const label len);
-
         //- Calculate the geometry for the patches (transformation tensors etc.)
         void calcGeometry();
 

src/OpenFOAM/primitives/triad/triad.C

@@ -66,6 +66,13 @@ const triad triad::min
     vector(-VGREAT, -VGREAT, -VGREAT)
 );
 
+const triad triad::I
+(
+    vector(1, 0, 0),
+    vector(0, 1, 0),
+    vector(0, 0, 1)
+);
+
 const triad triad::unset
 (
     vector(VGREAT, VGREAT, VGREAT),

src/OpenFOAM/primitives/triad/triad.H

@@ -100,6 +100,7 @@ public:
         static const triad one;
         static const triad max;
         static const triad min;
+        static const triad I;
         static const triad unset;
 
 

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1371,6 +1371,12 @@ void Foam::fvMeshSubset::setLargeCellSubset
 }
 
 
+bool Foam::fvMeshSubset::hasSubMesh() const
+{
+    return fvMeshSubsetPtr_.valid();
+}
+
+
 const fvMesh& Foam::fvMeshSubset::subMesh() const
 {
     checkCellSubset();

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -268,6 +268,9 @@ public:
                 return baseMesh_;
             }
 
+            //- Have subMesh?
+            bool hasSubMesh() const;
+
             //- Return reference to subset mesh
             const fvMesh& subMesh() const;
 

tutorials/mesh/foamyHexMesh/blob/Allclean

@@ -6,9 +6,20 @@ cd ${0%/*} || exit 1    # run from this directory
 
 cleanCase
 
+rm -rf constant/triSurface/blob.stl.gz > /dev/null 2>&1
+
 rm -rf constant/cellAlignments > /dev/null 2>&1
 rm -rf constant/targetCellSize > /dev/null 2>&1
 rm -rf constant/internalDelaunayVertices > /dev/null 2>&1
 
+rm -rf constant/backgroundMeshDecomposition/polyMesh/boundary > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/faces > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/neighbour > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/owner > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/points > /dev/null 2>&1
+
+rm -rf snapToSurface?.obj > /dev/null 2>&1
+rm -rf tetsToSnapTo.obj > /dev/null 2>&1
+
 
 # ----------------------------------------------------------------- end-of-file

tutorials/mesh/foamyHexMesh/flange/Allclean

@@ -4,11 +4,21 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
+rm -r constant/triSurface/flange.stl.gz > /dev/null 2>&1
+
 rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
 rm -r constant/ccx constant/ccy constant/ccz > /dev/null 2>&1
 rm -r constant/internalDelaunayVertices constant/targetCellSize > /dev/null 2>&1
 rm -r 0/ccx 0/ccy 0/ccz > /dev/null 2>&1
 
+rm -rf constant/backgroundMeshDecomposition/polyMesh/boundary > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/faces > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/neighbour > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/owner > /dev/null 2>&1
+rm -rf constant/backgroundMeshDecomposition/polyMesh/points > /dev/null 2>&1
+
+rm -r *.obj > /dev/null 2>&1
+
 cleanCase
 
 # ----------------------------------------------------------------- end-of-file

tutorials/mesh/foamyHexMesh/mixerVessel/Allclean

@@ -6,6 +6,23 @@ cd ${0%/*} || exit 1    # run from this directory
 
 rm -r constant/internalDelaunayVertices constant/targetCellSize > /dev/null 2>&1
 
+rm -r snapToSurface?.obj subsetBox.obj surf1.obj tetsToSnapTo.obj > /dev/null 2>&1
+
+rm -r constant/backgroundMeshDecomposition/polyMesh/boundary > /dev/null 2>&1
+rm -r constant/backgroundMeshDecomposition/polyMesh/faces > /dev/null 2>&1
+rm -r constant/backgroundMeshDecomposition/polyMesh/neighbour > /dev/null 2>&1
+rm -r constant/backgroundMeshDecomposition/polyMesh/owner > /dev/null 2>&1
+rm -r constant/backgroundMeshDecomposition/polyMesh/points > /dev/null 2>&1
+
+rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
+
+rm -r constant/triSurface/rawSurfaces/*.stl > /dev/null 2>&1
+
+rm -r constant/triSurface/*.stl > /dev/null 2>&1
+rm -r constant/triSurface/*.obj > /dev/null 2>&1
+rm -r constant/triSurface/*.eMesh > /dev/null 2>&1
+rm -r constant/triSurface/problemFaces > /dev/null 2>&1
+
 cleanCase
 
 # ----------------------------------------------------------------- end-of-file

tutorials/mesh/foamyHexMesh/simpleShapes/Allclean

@@ -6,11 +6,13 @@ cd ${0%/*} || exit 1    # run from this directory
 
 rm -r constant/extendedFeatureEdgeMesh constant/internalDelaunayVertices > /dev/null 2>&1
 rm constant/triSurface/*.eMesh > /dev/null 2>&1
-rm constant/triSurface/*_clean* > /dev/null 2>&1
+rm constant/triSurface/*_orient* > /dev/null 2>&1
+rm -r constant/polyMesh > /dev/null 2>&1
 rm -r constant/polyMesh > /dev/null 2>&1
-rm constant/targetCellSize > /dev/null 2>&1
 rm -r constant/tetDualMesh > /dev/null 2>&1
 
+rm -r snapToSurface?.obj tetsToSnapTo.obj > /dev/null 2>&1
+
 rm domain coneAndSphere > /dev/null 2>&1
 
 cleanCase