ENH: Add options to feature point handling

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

@@ -111,6 +111,8 @@ public:
 
     typedef Delaunay::Vertex_handle    Vertex_handle;
     typedef Delaunay::Cell_handle      Cell_handle;
+    typedef Delaunay::Edge             Edge;
+    typedef Delaunay::Facet            Facet;
     typedef Delaunay::Point            Point;
 
     typedef List<DynamicList<Pair<labelPair> > > labelPairPairDynListList;
@@ -371,6 +373,28 @@ private:
             DynamicList<Vb>& pts
         );
 
+        void createEdgePointGroupByCirculating
+        (
+            const extendedFeatureEdgeMesh& feMesh,
+            const pointIndexHit& edHit,
+            DynamicList<Vb>& pts
+        );
+
+        bool meshableRegion
+        (
+            const plane::side side,
+            const extendedFeatureEdgeMesh::sideVolumeType volType
+        ) const;
+
+        bool regionIsInside
+        (
+            const extendedFeatureEdgeMesh::sideVolumeType volTypeA,
+            const vector& normalA,
+            const extendedFeatureEdgeMesh::sideVolumeType volTypeB,
+            const vector& normalB,
+            const vector& masterPtVec
+        ) const;
+
         //- Create points to conform to an external edge
         void createExternalEdgePointGroup
         (
@@ -872,9 +896,14 @@ private:
             PtrList<dictionary>& patchDicts
         ) const;
 
+        vector calcSharedPatchNormal(Cell_handle c1, Cell_handle c2) const;
+
+        bool boundaryDualFace(Cell_handle c1, Cell_handle c2) const;
+
         //- Create all of the internal and boundary faces
         void createFacesOwnerNeighbourAndPatches
         (
+            const pointField& pts,
             faceList& faces,
             labelList& owner,
             labelList& neighbour,

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

@@ -536,6 +536,39 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 //        checkDuals();
 //
 //        Info<< nl << "Finished checks" << nl << endl;
+//    }
+
+//    OFstream str("attachedToFeature.obj");
+//    label offset = 0;
+//
+//    for
+//    (
+//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+//        vit != finite_vertices_end();
+//        ++vit
+//    )
+//    {
+//        if (vit->featurePoint())
+//        {
+//            std::list<Cell_handle> adjacentCells;
+//
+//            finite_incident_cells(vit, std::back_inserter(adjacentCells));
+//
+//            for
+//            (
+//                std::list<Cell_handle>::iterator acit = adjacentCells.begin();
+//                acit != adjacentCells.end();
+//                ++acit
+//            )
+//            {
+//                if ((*acit)->real())
+//                {
+//                    drawDelaunayCell(str, (*acit), offset);
+//                    offset++;
+////                    meshTools::writeOBJ(str, topoint((*acit)->dual()));
+//                }
+//            }
+//        }
 //    }
 
     setVertexSizeAndAlignment();
@@ -557,6 +590,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 
     createFacesOwnerNeighbourAndPatches
     (
+        points,
         faces,
         owner,
         neighbour,
@@ -1190,6 +1224,7 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
 
     createFacesOwnerNeighbourAndPatches
     (
+        pts,
         faces,
         owner,
         neighbour,
@@ -1703,9 +1738,53 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
     this->resetCellCount();
 
-    pts.setSize(number_of_finite_cells());
+    label nConstrainedVertices = 0;
+    if (foamyHexMeshControls().guardFeaturePoints())
+    {
+        for
+        (
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            ++vit
+        )
+        {
+            if (vit->constrained())
+            {
+                vit->index() = number_of_finite_cells() + nConstrainedVertices;
+                nConstrainedVertices++;
+            }
+        }
+    }
 
-    boundaryPts.setSize(number_of_finite_cells(), -1);
+    pts.setSize(number_of_finite_cells() + nConstrainedVertices);
+    boundaryPts.setSize
+    (
+        number_of_finite_cells() + nConstrainedVertices,
+        -1
+    );
+
+    if (foamyHexMeshControls().guardFeaturePoints())
+    {
+        nConstrainedVertices = 0;
+        for
+        (
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            ++vit
+        )
+        {
+            if (vit->constrained())
+            {
+                pts[number_of_finite_cells() + nConstrainedVertices] =
+                    topoint(vit->point());
+
+                boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
+                    1;
+
+                nConstrainedVertices++;
+            }
+        }
+    }
 
     for
     (
@@ -1780,28 +1859,35 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                 pts[cit->cellIndex()] = cit->dual();
             }
 
-            if (cit->featurePointDualVertex())
+            // Feature point snapping
+            if (foamyHexMeshControls().snapFeaturePoints())
             {
-                pointFromPoint dual = cit->dual();
+                if (cit->featurePointDualVertex())
-
-                pointIndexHit fpHit;
-                label featureHit;
-
-                // Find nearest feature point and compare
-                geometryToConformTo_.findFeaturePointNearest
-                (
-                    dual,
-                    sqr(targetCellSize(dual)),
-                    fpHit,
-                    featureHit
-                );
-
-                if (fpHit.hit())
                 {
-                    Info<< "Dual        = " << dual << nl
+                    pointFromPoint dual = cit->dual();
-                        << "    Nearest = " << fpHit.hitPoint() << endl;
 
-                    pts[cit->cellIndex()] = fpHit.hitPoint();
+                    pointIndexHit fpHit;
+                    label featureHit;
+
+                    // Find nearest feature point and compare
+                    geometryToConformTo_.findFeaturePointNearest
+                    (
+                        dual,
+                        sqr(targetCellSize(dual)),
+                        fpHit,
+                        featureHit
+                    );
+
+                    if (fpHit.hit())
+                    {
+                        if (debug)
+                        {
+                            Info<< "Dual        = " << dual << nl
+                                << "    Nearest = " << fpHit.hitPoint() << endl;
+                        }
+
+                        pts[cit->cellIndex()] = fpHit.hitPoint();
+                    }
                 }
             }
 
@@ -1823,9 +1909,9 @@ void Foam::conformalVoronoiMesh::indexDualVertices
         }
     }
 
-    pts.setSize(this->cellCount());
+    //pts.setSize(this->cellCount());
 
-    boundaryPts.setSize(this->cellCount());
+    //boundaryPts.setSize(this->cellCount());
 }
 
 
@@ -1873,9 +1959,13 @@ Foam::label Foam::conformalVoronoiMesh::createPatchInfo
                     "type",
                     wallPolyPatch::typeName
                 );
-        }
 
-        patchDicts.set(patchI, new dictionary());
+            patchDicts.set(patchI, new dictionary(patchInfo[patchI]));
+        }
+        else
+        {
+            patchDicts.set(patchI, new dictionary());
+        }
     }
 
     patchNames.setSize(patchNames.size() + 1);
@@ -1977,8 +2067,71 @@ Foam::label Foam::conformalVoronoiMesh::createPatchInfo
 }
 
 
+Foam::vector Foam::conformalVoronoiMesh::calcSharedPatchNormal
+(
+    Cell_handle c1,
+    Cell_handle c2
+) const
+{
+    List<Foam::point> patchEdge(2, point::max);
+
+    label count = 0;
+
+    // Get shared Facet
+    for (label cI = 0; cI < 4; ++cI)
+    {
+        if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
+        {
+            if (c1->vertex(cI)->internalBoundaryPoint())
+            {
+                patchEdge[0] = topoint(c1->vertex(cI)->point());
+            }
+            else
+            {
+                patchEdge[1] = topoint(c1->vertex(cI)->point());
+            }
+        }
+    }
+
+    Info<< "    " << patchEdge << endl;
+
+    return vector(patchEdge[1] - patchEdge[0]);
+}
+
+
+bool Foam::conformalVoronoiMesh::boundaryDualFace
+(
+    Cell_handle c1,
+    Cell_handle c2
+) const
+{
+    label nInternal = 0;
+    label nExternal = 0;
+
+    for (label cI = 0; cI < 4; ++cI)
+    {
+        if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
+        {
+            if (c1->vertex(cI)->internalBoundaryPoint())
+            {
+                nInternal++;
+            }
+            else if (c1->vertex(cI)->externalBoundaryPoint())
+            {
+                nExternal++;
+            }
+        }
+    }
+
+    Info<< "in = " << nInternal << " out = " << nExternal << endl;
+
+    return (nInternal == 1 && nExternal == 1);
+}
+
+
 void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 (
+    const pointField& pts,
     faceList& faces,
     labelList& owner,
     labelList& neighbour,
@@ -2023,11 +2176,244 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
     faces.setSize(number_of_finite_edges());
     owner.setSize(number_of_finite_edges());
     neighbour.setSize(number_of_finite_edges());
+    boundaryFacesToRemove.setSize(number_of_finite_edges(), false);
 
     labelPairPairDynListList procPatchSortingIndex(nPatches);
 
     label dualFaceI = 0;
 
+    if (foamyHexMeshControls().guardFeaturePoints())
+    {
+        OBJstream startCellStr("startingCell.obj");
+        OBJstream featurePointFacesStr("ftPtFaces.obj");
+        OBJstream featurePointDualsStr("ftPtDuals.obj");
+        OFstream cellStr("vertexCells.obj");
+
+        label vcount = 1;
+
+        for
+        (
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            ++vit
+        )
+        {
+            if (vit->constrained())
+            {
+                // Find a starting cell
+                std::list<Cell_handle> vertexCells;
+                finite_incident_cells(vit, std::back_inserter(vertexCells));
+
+                Cell_handle startCell;
+
+                for
+                (
+                    std::list<Cell_handle>::iterator vcit = vertexCells.begin();
+                    vcit != vertexCells.end();
+                    ++vcit
+                )
+                {
+                    if ((*vcit)->featurePointExternalCell())
+                    {
+                        startCell = *vcit;
+                    }
+
+                    if ((*vcit)->real())
+                    {
+                        featurePointDualsStr.write
+                        (
+                            linePointRef(topoint(vit->point()), (*vcit)->dual())
+                        );
+                    }
+                }
+
+                // Error if startCell is null
+                if (startCell == NULL)
+                {
+                    Pout<< "Start cell is null!" << endl;
+                }
+
+                // Need to pick a direction to walk in
+                Cell_handle vc1 = startCell;
+                Cell_handle vc2;
+
+                Info<< "c1 index = " << vc1->cellIndex() << " "
+                    << vc1->dual() << endl;
+
+                for (label cI = 0; cI < 4; ++cI)
+                {
+                    Info<< "c1 = " << cI << " "
+                        << vc1->neighbor(cI)->cellIndex() << " v = "
+                        << vc1->neighbor(cI)->dual() << endl;
+
+                    Info<< vc1->vertex(cI)->info();
+                }
+
+                Cell_handle nextCell;
+
+                for (label cI = 0; cI < 4; ++cI)
+                {
+                    if (vc1->vertex(cI)->externalBoundaryPoint())
+                    {
+                        vc2 = vc1->neighbor(cI);
+
+                        Info<< "    c2 is neighbor "
+                            << vc2->cellIndex()
+                            << " of c1" << endl;
+
+                        for (label cI = 0; cI < 4; ++cI)
+                        {
+                            Info<< "    c2 = " << cI << " "
+                                << vc2->neighbor(cI)->cellIndex() << " v = "
+                                << vc2->vertex(cI)->index() << endl;
+                        }
+
+                        face f(3);
+                        f[0] = vit->index();
+                        f[1] = vc1->cellIndex();
+                        f[2] = vc2->cellIndex();
+
+                        Info<< "f " << f << endl;
+                        forAll(f, pI)
+                        {
+                            Info<< "    " << pts[f[pI]] << endl;
+                        }
+
+                        vector correctNormal = calcSharedPatchNormal(vc1, vc2);
+                        correctNormal /= mag(correctNormal);
+
+                        Info<< "    cN " << correctNormal << endl;
+
+                        vector fN = f.normal(pts);
+
+                        if (mag(fN) < SMALL)
+                        {
+                            nextCell = vc2;
+                            continue;
+                        }
+
+                        fN /= mag(fN);
+                        Info<< "    fN " << fN << endl;
+
+                        if ((fN & correctNormal) > 0)
+                        {
+                            nextCell = vc2;
+                            break;
+                        }
+                    }
+                }
+
+                vc2 = nextCell;
+
+                label own = vit->index();
+                face f(3);
+                f[0] = own;
+
+                Info<< "Start walk from " << vc1->cellIndex()
+                    << " to " << vc2->cellIndex() << endl;
+
+                // Walk while not at start cell
+
+                label iter = 0;
+                do
+                {
+                    Info<< "     Walk from " << vc1->cellIndex()
+                        << " " << vc1->dual()
+                        << " to " << vc2->cellIndex()
+                        << " " << vc2->dual()
+                        << endl;
+
+                    startCellStr.write(linePointRef(vc1->dual(), vc2->dual()));
+
+                    // Get patch by getting face between cells and the two
+                    // points on the face that are not the feature vertex
+                    label patchIndex =
+                        geometryToConformTo_.findPatch
+                        (
+                            topoint(vit->point())
+                        );
+
+                    f[1] = vc1->cellIndex();
+                    f[2] = vc2->cellIndex();
+
+                    patchFaces[patchIndex].append(f);
+                    patchOwners[patchIndex].append(own);
+                    patchPPSlaves[patchIndex].append(own);
+
+                    // Find next cell
+                    Cell_handle nextCell;
+
+                    Info<< "    c1 vertices " << vc2->dual() << endl;
+                    for (label cI = 0; cI < 4; ++cI)
+                    {
+                        Info<< "        " << vc2->vertex(cI)->info();
+                    }
+                    Info<< "    c1 neighbour vertices " << endl;
+                    for (label cI = 0; cI < 4; ++cI)
+                    {
+                        if
+                        (
+                            !vc2->vertex(cI)->constrained()
+                         && vc2->neighbor(cI) != vc1
+                         && !is_infinite(vc2->neighbor(cI))
+                         &&
+                            (
+                                vc2->neighbor(cI)->featurePointExternalCell()
+                             || vc2->neighbor(cI)->featurePointInternalCell()
+                            )
+                         && vc2->neighbor(cI)->hasConstrainedPoint()
+                        )
+                        {
+                            drawDelaunayCell
+                            (
+                                cellStr,
+                                vc2->neighbor(cI),
+                                vcount++
+                            );
+
+                            Info<< "        neighbour " << cI << " "
+                                << vc2->neighbor(cI)->dual() << endl;
+                            for (label I = 0; I < 4; ++I)
+                            {
+                                Info<< "            "
+                                    << vc2->neighbor(cI)->vertex(I)->info();
+                            }
+                        }
+                    }
+
+                    for (label cI = 0; cI < 4; ++cI)
+                    {
+                        if
+                        (
+                            !vc2->vertex(cI)->constrained()
+                         && vc2->neighbor(cI) != vc1
+                         && !is_infinite(vc2->neighbor(cI))
+                         &&
+                            (
+                                vc2->neighbor(cI)->featurePointExternalCell()
+                             || vc2->neighbor(cI)->featurePointInternalCell()
+                            )
+                         && vc2->neighbor(cI)->hasConstrainedPoint()
+                        )
+                        {
+                            // check if shared edge is internal/internal
+                            if (boundaryDualFace(vc2, vc2->neighbor(cI)))
+                            {
+                                nextCell = vc2->neighbor(cI);
+                                break;
+                            }
+                        }
+                    }
+
+                    vc1 = vc2;
+                    vc2 = nextCell;
+
+                    iter++;
+                } while (vc1 != startCell && iter < 100);
+            }
+        }
+    }
+
     for
     (
         Delaunay::Finite_edges_iterator eit = finite_edges_begin();
@@ -2039,7 +2425,35 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
         Vertex_handle vA = c->vertex(eit->second);
         Vertex_handle vB = c->vertex(eit->third);
 
+        if (vA->constrained() && vB->constrained())
+        {
+            continue;
+        }
+
         if
+        (
+            (vA->constrained() && vB->internalOrBoundaryPoint())
+         || (vB->constrained() && vA->internalOrBoundaryPoint())
+        )
+        {
+            face newDualFace = buildDualFace(eit);
+
+            label own = -1;
+            label nei = -1;
+
+            if (ownerAndNeighbour(vA, vB, own, nei))
+            {
+                reverse(newDualFace);
+            }
+
+            // internal face
+            faces[dualFaceI] = newDualFace;
+            owner[dualFaceI] = own;
+            neighbour[dualFaceI] = nei;
+
+            dualFaceI++;
+        }
+        else if
         (
             (vA->internalOrBoundaryPoint() && !vA->referred())
          || (vB->internalOrBoundaryPoint() && !vB->referred())
@@ -2057,15 +2471,15 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                     reverse(newDualFace);
                 }
 
+                label patchIndex = -1;
+
+                pointFromPoint ptA = topoint(vA->point());
+                pointFromPoint ptB = topoint(vB->point());
+
                 if (nei == -1)
                 {
                     // boundary face
 
-                    pointFromPoint ptA = topoint(vA->point());
-                    pointFromPoint ptB = topoint(vB->point());
-
-                    label patchIndex = -1;
-
                     if (isProcBoundaryEdge(eit))
                     {
                         // One (and only one) of the points is an internal
@@ -2487,6 +2901,7 @@ void Foam::conformalVoronoiMesh::addPatches
 
     faces.setSize(nInternalFaces + nBoundaryFaces);
     owner.setSize(nInternalFaces + nBoundaryFaces);
+    boundaryFacesToRemove.setSize(nInternalFaces + nBoundaryFaces);
 
     label faceI = nInternalFaces;
 

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

@@ -40,45 +40,439 @@ void Foam::conformalVoronoiMesh::createEdgePointGroup
     DynamicList<Vb>& pts
 )
 {
-    label edgeI = edHit.index();
+    if (foamyHexMeshControls().circulateEdges())
-
-    extendedFeatureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
-
-    switch (edStatus)
     {
-        case extendedFeatureEdgeMesh::EXTERNAL:
+        createEdgePointGroupByCirculating(feMesh, edHit, pts);
+    }
+    else
+    {
+        label edgeI = edHit.index();
+
+        extendedFeatureEdgeMesh::edgeStatus edStatus =
+            feMesh.getEdgeStatus(edgeI);
+
+        switch (edStatus)
         {
-            createExternalEdgePointGroup(feMesh, edHit, pts);
+            case extendedFeatureEdgeMesh::EXTERNAL:
-            break;
+            {
-        }
+                createExternalEdgePointGroup(feMesh, edHit, pts);
-        case extendedFeatureEdgeMesh::INTERNAL:
+                break;
-        {
+            }
-            createInternalEdgePointGroup(feMesh, edHit, pts);
+            case extendedFeatureEdgeMesh::INTERNAL:
-            break;
+            {
-        }
+                createInternalEdgePointGroup(feMesh, edHit, pts);
-        case extendedFeatureEdgeMesh::FLAT:
+                break;
-        {
+            }
-            createFlatEdgePointGroup(feMesh, edHit, pts);
+            case extendedFeatureEdgeMesh::FLAT:
-            break;
+            {
-        }
+                createFlatEdgePointGroup(feMesh, edHit, pts);
-        case extendedFeatureEdgeMesh::OPEN:
+                break;
-        {
+            }
-            createOpenEdgePointGroup(feMesh, edHit, pts);
+            case extendedFeatureEdgeMesh::OPEN:
-            break;
+            {
-        }
+                createOpenEdgePointGroup(feMesh, edHit, pts);
-        case extendedFeatureEdgeMesh::MULTIPLE:
+                break;
-        {
+            }
-            createMultipleEdgePointGroup(feMesh, edHit, pts);
+            case extendedFeatureEdgeMesh::MULTIPLE:
-            break;
+            {
-        }
+                createMultipleEdgePointGroup(feMesh, edHit, pts);
-        case extendedFeatureEdgeMesh::NONE:
+                break;
-        {
+            }
-            break;
+            case extendedFeatureEdgeMesh::NONE:
+            {
+                break;
+            }
         }
     }
 }
 
 
+bool Foam::conformalVoronoiMesh::meshableRegion
+(
+    const plane::side side,
+    const extendedFeatureEdgeMesh::sideVolumeType volType
+) const
+{
+    switch (volType)
+    {
+        case extendedFeatureEdgeMesh::INSIDE:
+        {
+            return (side == plane::FLIP) ? true : false;
+        }
+        case extendedFeatureEdgeMesh::OUTSIDE:
+        {
+            return (side == plane::NORMAL) ? true : false;
+        }
+        case extendedFeatureEdgeMesh::BOTH:
+        {
+            return true;
+        }
+        case extendedFeatureEdgeMesh::NEITHER:
+        {
+            return false;
+        }
+    }
+}
+
+
+bool Foam::conformalVoronoiMesh::regionIsInside
+(
+    const extendedFeatureEdgeMesh::sideVolumeType volTypeA,
+    const vector& normalA,
+    const extendedFeatureEdgeMesh::sideVolumeType volTypeB,
+    const vector& normalB,
+    const vector& masterPtVec
+) const
+{
+    plane::side sideA
+    (
+        ((masterPtVec & normalA) <= 0) ? plane::FLIP : plane::NORMAL
+    );
+
+    plane::side sideB
+    (
+        ((masterPtVec & normalB) <= 0) ? plane::FLIP : plane::NORMAL
+    );
+
+    const bool meshableRegionA = meshableRegion(sideA, volTypeA);
+    const bool meshableRegionB = meshableRegion(sideB, volTypeB);
+
+    if (meshableRegionA == meshableRegionB)
+    {
+        return meshableRegionA;
+    }
+    else
+    {
+        WarningIn
+        (
+            "Foam::conformalVoronoiMesh::regionIsInside"
+            "(volTypeA, normalA, volTypeB, normalB, masterPtVec)"
+        )   << ""
+            << endl;
+
+        return false;
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Vb>& pts
+)
+{
+    typedef Foam::indexedVertexEnum::vertexType         vertexType;
+    typedef extendedFeatureEdgeMesh::sideVolumeType     sideVolumeType;
+
+    const Foam::point& edgePt = edHit.hitPoint();
+    const label edgeI = edHit.index();
+
+    scalar ppDist = pointPairDistance(edgePt);
+
+    const vectorField& feNormals = feMesh.normals();
+    const vector& edDir = feMesh.edgeDirections()[edgeI];
+    const labelList& edNormalIs = feMesh.edgeNormals()[edgeI];
+    const labelList& feNormalDirections = feMesh.normalDirections()[edgeI];
+
+    const PackedList<2>& normalVolumeTypes = feMesh.normalVolumeTypes();
+
+    const_circulator<labelList> circ(edNormalIs);
+    const_circulator<labelList> circNormalDirs(feNormalDirections);
+
+    Map<Foam::point> masterPoints;
+    Map<vertexType> masterPointsTypes;
+    Map<plane> masterPointReflectionsPrev;
+    Map<plane> masterPointReflectionsNext;
+
+//    Info<< "Edge = " << edHit << ", edDir = " << edDir << endl;
+//    Info<< "    edNorms = " << edNormalIs.size() << endl;
+
+    bool addedMasterPreviously = false;
+    label initialRegion = -1;
+
+    if (circ.size()) do
+    {
+        const sideVolumeType volType =
+            sideVolumeType(normalVolumeTypes[circ()]);
+        const sideVolumeType nextVolType =
+            sideVolumeType(normalVolumeTypes[circ.next()]);
+
+        const vector& normal = feNormals[circ()];
+        const vector& nextNormal = feNormals[circ.next()];
+
+        vector normalDir = (normal ^ edDir);
+        normalDir *= circNormalDirs()/mag(normalDir);
+
+        vector nextNormalDir = (nextNormal ^ edDir);
+        nextNormalDir *= circNormalDirs.next()/mag(nextNormalDir);
+
+//        Info<< "    " << circ() << " " << circ.next() << nl
+//            << "    " << circNormalDirs() << " " << circNormalDirs.next()
+//            << nl
+//            << "    normals = " << normal << " " << nextNormal << nl
+//            << " normalDirs = " << normalDir << " " << nextNormalDir << nl
+//            << "      cross = " << (normalDir ^ nextNormalDir) << nl
+//            << "    "
+//            << extendedFeatureEdgeMesh::sideVolumeTypeNames_[volType] << " "
+//            << extendedFeatureEdgeMesh::sideVolumeTypeNames_[nextVolType]
+//            << endl;
+
+        // Calculate master point
+        vector masterPtVec(normalDir + nextNormalDir);
+        masterPtVec /= mag(masterPtVec) + SMALL;
+
+        if (((normalDir ^ nextNormalDir) & edDir) < SMALL)
+        {
+//            Info<< "    IGNORE REGION" << endl;
+            addedMasterPreviously = false;
+
+            if (circ.size() == 2 && mag((normal & nextNormal) - 1) < SMALL)
+            {
+                // Add an extra point
+
+                const vector n = 0.5*(normal + nextNormal);
+
+                vector s = ppDist*(edDir ^ n);
+
+                createPointPair(ppDist, edgePt + s, n, pts);
+                createPointPair(ppDist, edgePt - s, n, pts);
+
+                break;
+            }
+
+            continue;
+        }
+
+        if (mag(masterPtVec) < SMALL)
+        {
+            if (circ.size() == 2)
+            {
+                // Add an extra point
+
+                // Average normal to remove any bias to one side, although as
+                // it is a flat edge, the normals should be essentially the same
+                const vector n = 0.5*(normal + nextNormal);
+
+                // Direction along the surface to the control point, sense of
+                // direction not important, as +s and -s can be used because it
+                // is a flat edge
+                vector s = ppDist*(edDir ^ n);
+
+                createPointPair(ppDist, edgePt + s, n, pts);
+                createPointPair(ppDist, edgePt - s, n, pts);
+
+                break;
+            }
+            else
+            {
+//                Info<< "    IGNORE REGION" << endl;
+                addedMasterPreviously = false;
+                continue;
+            }
+        }
+
+        const Foam::point masterPt = edgePt + ppDist*masterPtVec;
+
+        // Check that region is inside or outside
+        const bool inside =
+            regionIsInside
+            (
+                volType,
+                normal,
+                nextVolType,
+                nextNormal,
+                masterPtVec
+            );
+
+        // Specialise for size = 1 && baffle
+        if (mag((normalDir & nextNormalDir) - 1) < SMALL)
+        {
+            if (inside)
+            {
+//                Info<< "Specialise for size 1 and baffle" << endl;
+
+                vector s = ppDist*(edDir ^ normal);
+
+                plane facePlane(edgePt, normal);
+
+                Foam::point pt1 = edgePt + s + ppDist*normal;
+                Foam::point pt2 = edgePt - s + ppDist*normal;
+
+                Foam::point pt3 = reflectPointInPlane(pt1, facePlane);
+                Foam::point pt4 = reflectPointInPlane(pt2, facePlane);
+
+                pts.append(Vb(pt1, Vb::vtInternalFeatureEdge));
+                pts.append(Vb(pt2, Vb::vtInternalFeatureEdge));
+                pts.append(Vb(pt3, Vb::vtInternalFeatureEdge));
+                pts.append(Vb(pt4, Vb::vtInternalFeatureEdge));
+
+                break;
+            }
+            else
+            {
+                WarningIn
+                (
+                    "Foam::conformalVoronoiMesh::"
+                    "createEdgePointGroupByCirculating"
+                    "("
+                    "    const extendedFeatureEdgeMesh&,"
+                    "    const pointIndexHit&,"
+                    "    DynamicList<Vb>&"
+                    ")"
+                )   << "Faces are parallel but master point is not inside"
+                    << endl;
+            }
+        }
+
+        if (!addedMasterPreviously)
+        {
+            if (initialRegion == -1)
+            {
+                initialRegion = circ.nRotations();
+            }
+
+            addedMasterPreviously = true;
+
+            masterPoints.insert(circ.nRotations(), masterPt);
+            masterPointsTypes.insert
+            (
+                circ.nRotations(),
+                inside
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            );
+
+            masterPointReflectionsPrev.insert
+            (
+                circ.nRotations(),
+                plane(edgePt, normal)
+            );
+
+            masterPointReflectionsNext.insert
+            (
+                circ.nRotations(),
+                plane(edgePt, nextNormal)
+            );
+        }
+        else if (addedMasterPreviously)
+        {
+            addedMasterPreviously = true;
+
+            masterPointReflectionsNext.erase(circ.nRotations() - 1);
+
+            // Shift the master point to be normal to the plane between it and
+            // the previous master point
+            // Should be the intersection of the normal and the plane with the
+            // new master point in it.
+
+            plane p(masterPoints[circ.nRotations() - 1], normalDir);
+            plane::ray r(edgePt, masterPt - edgePt);
+
+            scalar cutPoint = p.normalIntersect(r);
+
+            masterPoints.insert
+            (
+                circ.nRotations(),
+                edgePt + cutPoint*(masterPt - edgePt)
+            );
+
+            masterPointsTypes.insert
+            (
+                circ.nRotations(),
+                inside
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            );
+
+            masterPointReflectionsNext.insert
+            (
+                circ.nRotations(),
+                plane(edgePt, nextNormal)
+            );
+        }
+
+        if
+        (
+            masterPoints.size() > 1
+         && inside
+         && circ.nRotations() == circ.size() - 1
+        )
+        {
+            if (initialRegion == 0)
+            {
+                plane p(masterPoints[initialRegion], nextNormalDir);
+                plane::ray r(edgePt, masterPt - edgePt);
+
+                scalar cutPoint = p.normalIntersect(r);
+
+                masterPoints[circ.nRotations()] =
+                    edgePt + cutPoint*(masterPt - edgePt);
+
+                // Remove the first reflection plane if we are no longer
+                // circulating
+                masterPointReflectionsPrev.erase(initialRegion);
+                masterPointReflectionsNext.erase(circ.nRotations());
+            }
+            else
+            {
+
+            }
+        }
+    }
+    while
+    (
+        circ.circulate(CirculatorBase::CLOCKWISE),
+        circNormalDirs.circulate(CirculatorBase::CLOCKWISE)
+    );
+
+
+    forAllConstIter(Map<Foam::point>, masterPoints, iter)
+    {
+        const Foam::point& pt = masterPoints[iter.key()];
+        const vertexType ptType = masterPointsTypes[iter.key()];
+
+//        Info<< "    Adding Master " << iter.key() << " " << pt << " "
+//            << indexedVertexEnum::vertexTypeNames_[ptType] << endl;
+
+        pts.append(Vb(pt, ptType));
+
+        const vertexType reflectedPtType =
+        (
+            ptType == Vb::vtInternalFeatureEdge
+          ? Vb::vtExternalFeatureEdge
+          : Vb::vtInternalFeatureEdge
+        );
+
+        if (masterPointReflectionsPrev.found(iter.key()))
+        {
+            const Foam::point reflectedPt =
+                reflectPointInPlane(pt, masterPointReflectionsPrev[iter.key()]);
+
+//            Info<< "        Adding Prev " << reflectedPt << " "
+//                << indexedVertexEnum::vertexTypeNames_[reflectedPtType]
+//                << endl;
+
+            pts.append(Vb(reflectedPt, reflectedPtType));
+        }
+
+        if (masterPointReflectionsNext.found(iter.key()))
+        {
+            const Foam::point reflectedPt =
+                reflectPointInPlane(pt, masterPointReflectionsNext[iter.key()]);
+
+//            Info<< "        Adding Next " << reflectedPt << " "
+//                << indexedVertexEnum::vertexTypeNames_[reflectedPtType]
+//                << endl;
+
+            pts.append(Vb(reflectedPt, reflectedPtType));
+        }
+    }
+
+//    pts.append(Vb(edgePt, Vb::vtExternalFeatureEdge));
+}
+
+
 void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
 (
     const extendedFeatureEdgeMesh& feMesh,
@@ -142,7 +536,13 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
 
     pts.append
     (
-        Vb(refPt, Vb::vtInternalFeatureEdge)
+        Vb
+        (
+            refPt,
+            vertexCount() + pts.size(),
+            Vb::vtInternalFeatureEdge,
+            Pstream::myProcNo()
+        )
     );
 
     // Insert the slave points by reflecting refPt in both faces.
@@ -151,13 +551,26 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
     Foam::point reflectedA = refPt + 2*ppDist*nA;
     pts.append
     (
-        Vb(reflectedA, Vb::vtExternalFeatureEdge)
+        Vb
+        (
+            reflectedA,
+            vertexCount() + pts.size(),
+            Vb::vtExternalFeatureEdge,
+            Pstream::myProcNo()
+        )
     );
 
     Foam::point reflectedB = refPt + 2*ppDist*nB;
     pts.append
     (
-        Vb(reflectedB, Vb::vtExternalFeatureEdge)
+        Vb
+        (
+            reflectedB,
+            vertexCount() + pts.size(),
+            Vb::vtExternalFeatureEdge,
+            Pstream::myProcNo()
+        )
+    );
     );
 }
 
@@ -184,7 +597,6 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
     {
         // The normals are nearly parallel, so this is too sharp a feature to
         // conform to.
-
         return;
     }
 
@@ -262,13 +674,38 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
     // Master A is inside.
     pts.append
     (
-        Vb(reflectedA, Vb::vtInternalFeatureEdge)
+        Vb
+        (
+            reflectedA,
+            vertexCount() + pts.size(),
+            Vb::vtInternalFeatureEdge,
+            Pstream::myProcNo()
+        )
     );
 
     // Master B is inside.
     pts.append
     (
-        Vb(reflectedB, Vb::vtInternalFeatureEdge)
+        Vb
+        (
+            reflectedB,
+            vertexCount() + pts.size(),
+            Vb::vtInternalFeatureEdge,
+            Pstream::myProcNo()
+        )
+    );
+
+    // Slave is outside.
+    pts.append
+    (
+        Vb
+        (
+            reflMasterPt,
+            vertexCount() + pts.size(),
+            Vb::vtExternalFeatureEdge,
+            Pstream::myProcNo()
+        )
+    );
     );
 
     if (nAddPoints == 1)
@@ -277,7 +714,13 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
         // i.e. the original reference point
         pts.append
         (
-            Vb(refPt, Vb::vtInternalFeatureEdge)
+            Vb
+            (
+                refPt,
+                vertexCount() + pts.size(),
+                Vb::vtInternalFeatureEdge,
+                Pstream::myProcNo()
+            )
         );
     }
     else if (nAddPoints == 2)
@@ -285,21 +728,27 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
         Foam::point reflectedAa = refPt + ppDist*nB;
         pts.append
         (
-            Vb(reflectedAa, Vb::vtInternalFeatureEdge)
+            Vb
+            (
+                reflectedAa,
+                vertexCount() + pts.size(),
+                Vb::vtInternalFeatureEdge,
+                Pstream::myProcNo()
+            )
         );
 
         Foam::point reflectedBb = refPt + ppDist*nA;
         pts.append
         (
-            Vb(reflectedBb, Vb::vtInternalFeatureEdge)
+            Vb
+            (
+                reflectedBb,
+                vertexCount() + pts.size(),
+                Vb::vtInternalFeatureEdge,
+                Pstream::myProcNo()
+            )
         );
     }
-
-    // Slave is outside.
-    pts.append
-    (
-        Vb(reflMasterPt, Vb::vtExternalFeatureEdge)
-    );
 }
 
 
@@ -342,31 +791,43 @@ void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
     DynamicList<Vb>& pts
 )
 {
-//    // Assume it is a baffle and insert flat edge point pairs
+    // Assume it is a baffle and insert flat edge point pairs
-//    const Foam::point& edgePt = edHit.hitPoint();
+    const Foam::point& edgePt = edHit.hitPoint();
-//
-//    const scalar ppDist = pointPairDistance(edgePt);
-//
-//    const vectorField& feNormals = feMesh.normals();
-//    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
-//
-//    // As this is a flat edge, there are two normals by definition
-//    const vector& nA = feNormals[edNormalIs[0]];
-//    const vector& nB = feNormals[edNormalIs[1]];
-//
-//    // Average normal to remove any bias to one side, although as this
-//    // is a flat edge, the normals should be essentially the same
-//    const vector n = 0.5*(nA + nB);
-//
-//    // Direction along the surface to the control point, sense of edge
-//    // direction not important, as +s and -s can be used because this
-//    // is a flat edge
-//    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
-//
-//    createBafflePointPair(ppDist, edgePt + s, n, pts);
-//    createBafflePointPair(ppDist, edgePt - s, n, pts);
 
-    Info<< "NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
+    const scalar ppDist = pointPairDistance(edgePt);
+
+    const vectorField& feNormals = feMesh.normals();
+    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+
+    if (edNormalIs.size() == 1)
+    {
+//        Info<< "Inserting open edge point group around " << edgePt << endl;
+//        Info<< "    ppDist    = " << ppDist << nl
+//            << "    edNormals = " << edNormalIs
+//            << endl;
+
+        const vector& n = feNormals[edNormalIs[0]];
+
+        vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+
+        plane facePlane(edgePt, n);
+
+        Foam::point pt1 = edgePt + s + ppDist*n;
+        Foam::point pt2 = edgePt - s + ppDist*n;
+
+        Foam::point pt3 = reflectPointInPlane(pt1, facePlane);
+        Foam::point pt4 = reflectPointInPlane(pt2, facePlane);
+
+        pts.append(Vb(pt1, Vb::vtInternalSurface));
+        pts.append(Vb(pt2, Vb::vtInternalSurface));
+        pts.append(Vb(pt3, Vb::vtInternalSurface));
+        pts.append(Vb(pt4, Vb::vtInternalSurface));
+    }
+    else
+    {
+        Info<< "NOT INSERTING OPEN EDGE POINT GROUP WITH MORE THAN 1 "
+            << "EDGE NORMAL, NOT IMPLEMENTED" << endl;
+    }
 }
 
 
@@ -377,7 +838,33 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
     DynamicList<Vb>& pts
 )
 {
-    Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
+//    Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
+
+    const Foam::point& edgePt = edHit.hitPoint();
+
+    const scalar ppDist = pointPairDistance(edgePt);
+
+    const vectorField& feNormals = feMesh.normals();
+    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+
+    Info<< edNormalIs.size() << endl;
+    Info<< feMesh.edgeBaffles(edHit.index()) << endl;
+
+    // As this is a flat edge, there are two normals by definition
+    const vector& nA = feNormals[edNormalIs[0]];
+    const vector& nB = feNormals[edNormalIs[1]];
+
+    // Average normal to remove any bias to one side, although as this
+    // is a flat edge, the normals should be essentially the same
+    const vector n = 0.5*(nA + nB);
+
+    // Direction along the surface to the control point, sense of edge
+    // direction not important, as +s and -s can be used because this
+    // is a flat edge
+    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+
+    createPointPair(ppDist, edgePt + s, n, pts);
+    createPointPair(ppDist, edgePt - s, n, pts);
 }
 
 
@@ -453,7 +940,7 @@ void Foam::conformalVoronoiMesh::createMixedFeaturePoints
                 );
             }
 
-            if (!specialisedSuccess)
+            if (!specialisedSuccess && foamyHexMeshControls().edgeAiming())
             {
                 // Specialisations available for some mixed feature points.  For
                 // non-specialised feature points, inserting mixed internal and
@@ -462,34 +949,34 @@ void Foam::conformalVoronoiMesh::createMixedFeaturePoints
                 // Skipping unsupported mixed feature point types
                 bool skipEdge = false;
 
-                forAll(pEds, e)
+//                forAll(pEds, e)
-                {
+//                {
-                    const label edgeI = pEds[e];
+//                    const label edgeI = pEds[e];
-
+//
-                    const extendedFeatureEdgeMesh::edgeStatus edStatus
+//                    const extendedFeatureEdgeMesh::edgeStatus edStatus
-                        = feMesh.getEdgeStatus(edgeI);
+//                        = feMesh.getEdgeStatus(edgeI);
-
+//
-                    if
+//                    if
-                    (
+//                    (
-                        edStatus == extendedFeatureEdgeMesh::OPEN
+//                        edStatus == extendedFeatureEdgeMesh::OPEN
-                     || edStatus == extendedFeatureEdgeMesh::MULTIPLE
+//                     || edStatus == extendedFeatureEdgeMesh::MULTIPLE
-                    )
+//                    )
-                    {
+//                    {
-                        Info<< "Edge type " << edStatus
+//                        Info<< "Edge type " << edStatus
-                            << " found for mixed feature point " << ptI
+//                            << " found for mixed feature point " << ptI
-                            << ". Not supported."
+//                            << ". Not supported."
-                            << endl;
+//                            << endl;
-
+//
-                        skipEdge = true;
+//                        skipEdge = true;
-                    }
+//                    }
-                }
+//                }
-
+//
-                if (skipEdge)
+//                if (skipEdge)
-                {
+//                {
-                    Info<< "Skipping point " << ptI << nl << endl;
+//                    Info<< "Skipping point " << ptI << nl << endl;
-
+//
-                    continue;
+//                    continue;
-                }
+//                }
 
 //                createFeaturePoints(feMesh, ptI, pts, types);
 
@@ -576,6 +1063,16 @@ 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();
@@ -591,7 +1088,7 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints()
     {
         Vb& pt = pts[pI];
 
-        if (pt.featureEdgePoint())
+        //if (pt.featureEdgePoint())
         {
             if (pt.internalBoundaryPoint())
             {
@@ -721,22 +1218,25 @@ void Foam::conformalVoronoiMesh::addMasterAndSlavePoints
     {
         // Append master to the list of points
 
-//        OFstream strMasters("fpm_" + name(ptI) + ".obj");
-//        OFstream strSlaves("fps_" + name(ptI) + ".obj");
-
         const Foam::point& masterPt = masterPoints[pI];
         const vertexType masterType = masterPointsTypes[pI];
 
+//        Info<< "    Master = " << masterPt << endl;
+
         pts.append
         (
             Vb
             (
                 masterPt,
-                masterType
+                vertexCount() + pts.size(),
+                masterType,
+                Pstream::myProcNo()
             )
         );
 
-//        meshTools::writeOBJ(strMasters, masterPt);
+        const label masterIndex = pts[pts.size() - 1].index();
+
+        //meshTools::writeOBJ(strMasters, masterPt);
 
         const planeDynList& masterPointPlanes = masterPointReflections[pI];
 
@@ -749,6 +1249,8 @@ void Foam::conformalVoronoiMesh::addMasterAndSlavePoints
             const Foam::point slavePt =
                 reflectPointInPlane(masterPt, reflPlane);
 
+//            Info<< "        Slave " << planeI << " = " << slavePt << endl;
+
             const vertexType slaveType =
             (
                 masterType == Vb::vtInternalFeaturePoint
@@ -761,11 +1263,14 @@ void Foam::conformalVoronoiMesh::addMasterAndSlavePoints
                 Vb
                 (
                     slavePt,
-                    slaveType
+                    vertexCount() + pts.size(),
+                    slaveType,
+                    Pstream::myProcNo()
                 )
             );
 
-//            meshTools::writeOBJ(strSlaves, slavePt);
+
+            //meshTools::writeOBJ(strSlaves, slavePt);
         }
     }
 }
@@ -802,7 +1307,7 @@ void Foam::conformalVoronoiMesh::createMasterAndSlavePoints
 
     const Foam::point& featPt = feMesh.points()[ptI];
 
-    if (!positionOnThisProc(featPt))
+    if (!positionOnThisProc(featPt) || geometryToConformTo_.outside(featPt))
     {
         return;
     }
@@ -816,18 +1321,23 @@ void Foam::conformalVoronoiMesh::createMasterAndSlavePoints
 
     const labelList& featPtEdges = feMesh.featurePointEdges()[ptI];
 
-    const_circulator<labelList> circ(featPtEdges);
+    pointFeatureEdgesTypes pointEdgeTypes(ptI);
 
-//    Info<< "Point = " << ptI << endl;
+    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 eStatusPrev = feMesh.getEdgeStatus(circ.prev());
         const edgeStatus eStatusCurr = feMesh.getEdgeStatus(circ());
 //        const edgeStatus eStatusNext = feMesh.getEdgeStatus(circ.next());
 
-//        Info<< "Prev = "
+//        Info<< "    Prev = "
 //            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusPrev]
 //            << " Curr = "
 //            << extendedFeatureEdgeMesh::edgeStatusNames_[eStatusCurr]
@@ -843,10 +1353,13 @@ void Foam::conformalVoronoiMesh::createMasterAndSlavePoints
 
         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;
@@ -861,7 +1374,7 @@ void Foam::conformalVoronoiMesh::createMasterAndSlavePoints
               : Vb::vtInternalFeaturePoint // false
             );
 
-            Info<< "    " << " " << firstPlane << endl;
+            //Info<< "    " << " " << firstPlane << endl;
 
 //            const Foam::point reflectedPoint = reflectPointInPlane
 //            (
@@ -952,11 +1465,33 @@ void Foam::conformalVoronoiMesh::createFeaturePoints(DynamicList<Vb>& pts)
         (
             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
+                    )
+                );
+            }
+        }
     }
 }
 
@@ -1151,10 +1686,20 @@ Foam::vector Foam::conformalVoronoiMesh::sharedFaceNormal
     const vector& B1 = feMesh.normals()[nextEdgeInormals[0]];
     const vector& B2 = feMesh.normals()[nextEdgeInormals[1]];
 
-    const scalar A1B1 = mag(A1 ^ B1);
+//    Info<< "        A1 = " << A1 << endl;
-    const scalar A1B2 = mag(A1 ^ B2);
+//    Info<< "        A2 = " << A2 << endl;
-    const scalar A2B1 = mag(A2 ^ B1);
+//    Info<< "        B1 = " << B1 << endl;
-    const scalar A2B2 = mag(A2 ^ B2);
+//    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)
     {

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

@@ -62,10 +62,6 @@ Foam::cvControls::cvControls
         surfDict.lookup("featureEdgeExclusionDistanceCoeff")
     );
 
-    specialiseFeaturePoints_ = Switch
-    (
-        surfDict.lookup("specialiseFeaturePoints")
-    );
 
     surfaceSearchDistanceCoeff_ = readScalar
     (
@@ -85,6 +81,40 @@ Foam::cvControls::cvControls
         readLabel(surfDict.lookup("surfaceConformationRebuildFrequency"))
     );
 
+
+    const dictionary& featurePointControlsDict
+    (
+        surfDict.subDict("featurePointControls")
+    );
+
+    specialiseFeaturePoints_ = Switch
+    (
+        featurePointControlsDict.lookup("specialiseFeaturePoints")
+    );
+
+    guardFeaturePoints_ = Switch
+    (
+        featurePointControlsDict.lookup("guardFeaturePoints")
+    );
+
+    edgeAiming_ = Switch
+    (
+        featurePointControlsDict.lookup("edgeAiming")
+    );
+
+    if (!guardFeaturePoints_)
+    {
+        snapFeaturePoints_ = Switch
+        (
+            featurePointControlsDict.lookup("snapFeaturePoints")
+        );
+    }
+
+    circulateEdges_ = Switch
+    (
+        featurePointControlsDict.lookup("circulateEdges")
+    );
+
     // Controls for coarse surface conformation
 
     const dictionary& conformationControlsDict

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

@@ -78,8 +78,6 @@ class cvControls
             //  fraction of the local target cell size
             scalar featureEdgeExclusionDistanceCoeff_;
 
-            //- Switch for using specialised feature points
-            Switch specialiseFeaturePoints_;
 
             //- Surface search distance coefficient - fraction of the local
             //  target cell size
@@ -97,6 +95,23 @@ class cvControls
             //- Now often to rebuild the surface conformation
             label surfaceConformationRebuildFrequency_;
 
+            // Controls for feature point conformation
+
+                //-
+                Switch guardFeaturePoints_;
+
+                //-
+                Switch edgeAiming_;
+
+                //-
+                Switch snapFeaturePoints_;
+
+                //-
+                Switch circulateEdges_;
+
+                //- Switch for using specialised feature points
+                Switch specialiseFeaturePoints_;
+
 
             // Controls for coarse surface conformation
 
@@ -231,6 +246,18 @@ public:
             //- Return the featureEdgeExclusionDistanceCoeff
             inline scalar featureEdgeExclusionDistanceCoeff() const;
 
+            //-
+            inline Switch guardFeaturePoints() const;
+
+            //-
+            inline Switch edgeAiming() const;
+
+            //-
+            inline Switch snapFeaturePoints() const;
+
+            //-
+            inline Switch circulateEdges() const;
+
             //- Return the surfacePtExclusionDistanceCoeff
             inline scalar surfacePtExclusionDistanceCoeff() const;
 

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

@@ -54,6 +54,28 @@ inline Foam::scalar Foam::cvControls::featureEdgeExclusionDistanceCoeff() const
     return featureEdgeExclusionDistanceCoeff_;
 }
 
+
+inline Foam::Switch Foam::cvControls::guardFeaturePoints() const
+{
+    return guardFeaturePoints_;
+}
+
+inline Foam::Switch Foam::cvControls::edgeAiming() const
+{
+    return edgeAiming_;
+}
+
+inline Foam::Switch Foam::cvControls::snapFeaturePoints() const
+{
+    return snapFeaturePoints_;
+}
+
+inline Foam::Switch Foam::cvControls::circulateEdges() const
+{
+    return circulateEdges_;
+}
+
+
 inline Foam::scalar Foam::cvControls::surfacePtExclusionDistanceCoeff() const
 {
     return surfacePtExclusionDistanceCoeff_;