ENH: code cleanup

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

@@ -65,6 +65,7 @@ public:
 
         typedef typename Triangulation::Cell_handle Cell_handle;
         typedef typename Triangulation::Vertex_handle Vertex_handle;
+        typedef typename Triangulation::Edge Edge;
         typedef typename Triangulation::Point Point;
         typedef typename Triangulation::Facet Facet;
 

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

@@ -885,19 +885,9 @@ Foam::DistributedDelaunayMesh<Triangulation>::rangeInsertReferredWithInfo
             {
                 Cell_handle conflictingCell = V[i];
 
-                bool hasNullVertex = false;
-                for (label vertexI = 0; vertexI < 4; ++vertexI)
-                {
-                    if (conflictingCell->vertex(vertexI) == NULL)
-                    {
-                        hasNullVertex = true;
-                    }
-                }
-
                 if
                 (
-
+                    Triangulation::dimension() < 3 // 2D triangulation
-                    hasNullVertex
                  ||
                     (
                         !Triangulation::is_infinite(conflictingCell)

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

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

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

@@ -44,9 +44,9 @@ Foam::cellAspectRatioControl::cellAspectRatioControl
         )
     )
 {
-    Info<< nl << "    Cell Aspect Ratio Control" << nl
+    Info<< nl << "Cell Aspect Ratio Control" << nl
-        << "        Ratio     : " << aspectRatio_ << nl
+        << "    Ratio     : " << aspectRatio_ << nl
-        << "        Direction : " << aspectRatioDirection_
+        << "    Direction : " << aspectRatioDirection_
         << endl;
 }
 

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

@@ -124,7 +124,17 @@ Foam::scalar Foam::cellShapeControl::cellSize(const point& pt) const
 //            }
 //        }
 
-        // Find nearest surface
+//        cellShapeControlMesh::Vertex_handle nearV =
+//            shapeControlMesh_.nearest_vertex_in_cell
+//            (
+//                toPoint<cellShapeControlMesh::Point>(pt),
+//                ch
+//            );
+//
+//        size = nearV->targetCellSize();
+
+        // Find nearest surface. This can be quite slow if there are a lot of
+        // surfaces
         size = sizeAndAlignment_.cellSize(pt);
     }
     else

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

@@ -174,7 +174,7 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
     searchableSurface_(geometryToConformTo.geometry()[surfaceName_]),
     geometryToConformTo_(geometryToConformTo),
     cellSizeFunctions_(1),
-    regionToCellSizeFunctions_(geometryToConformTo_.patchNames().size(), 0),
+    regionToCellSizeFunctions_(geometryToConformTo_.patchNames().size(), -1),
     maxPriority_(-1)
 {
     Info<< indent << "Master settings:" << endl;
@@ -287,6 +287,18 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
 
         cellSizeFunctions_.transfer(regionCellSizeFunctions);
     }
+    else
+    {
+        const wordList& regionNames = geometryToConformTo_.patchNames();
+
+        forAll(regionNames, regionI)
+        {
+            if (regionToCellSizeFunctions_[regionI] == -1)
+            {
+                regionToCellSizeFunctions_[regionI] = 0;
+            }
+        }
+    }
 
 
     forAll(cellSizeFunctions_, funcI)

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

@@ -99,6 +99,16 @@ Foam::cellSizeFunction::cellSizeFunction
 
         sideMode_ = rmBothsides;
     }
+
+    if (debug)
+    {
+        Info<< nl
+            << "Cell size function for surface " << surface.name()
+            << ", " << mode
+            << ", priority = " << priority_
+            << ", regions = " << regionIndices_
+            << endl;
+    }
 }
 
 

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

@@ -514,7 +514,7 @@ void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
         << (cellSizeMesh.is_valid() ? "valid" : "not valid!" )
         << endl;
 
-    if (!Pstream::parRun())
+    if (!Pstream::parRun() && foamyHexMeshControls().objOutput())
     {
         //cellSizeMesh.writeTriangulation();
         cellSizeMesh.write();
@@ -709,6 +709,44 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
             );
         }
     }
+
+//    OFstream str(runTime_.path()/"alignments_internal.obj");
+//
+//    for
+//    (
+//        Finite_vertices_iterator vit = finite_vertices_begin();
+//        vit != finite_vertices_end();
+//        ++vit
+//    )
+//    {
+//        if (!vit->farPoint())
+//        {
+//            // Write alignments
+//            const tensor& alignment = vit->alignment();
+//            pointFromPoint pt = topoint(vit->point());
+//
+//            if
+//            (
+//                alignment.x() == triad::unset[0]
+//             || alignment.y() == triad::unset[0]
+//             || alignment.z() == triad::unset[0]
+//            )
+//            {
+//                Info<< "Bad alignment = " << vit->info();
+//
+//                vit->alignment() = tensor::I;
+//
+//                Info<< "New alignment = " << vit->info();
+//
+//                continue;
+//            }
+//
+//            meshTools::writeOBJ(str, pt, alignment.x() + pt);
+//            meshTools::writeOBJ(str, pt, alignment.y() + pt);
+//            meshTools::writeOBJ(str, pt, alignment.z() + pt);
+//        }
+//    }
+
 }
 
 
@@ -754,7 +792,6 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
         }
 
         label cc1I = cc1->cellIndex();
-
         label cc2I = cc2->cellIndex();
 
         if (cc1I != cc2I)
@@ -768,7 +805,6 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
         }
 
         cc1++;
-
         cc2++;
 
     } while (cc1 != ccStart);
@@ -854,14 +890,20 @@ bool Foam::conformalVoronoiMesh::ownerAndNeighbour
 
     if (!vA->internalOrBoundaryPoint() || vA->referred())
     {
-        dualCellIndexA = -1;
+        if (!vA->constrained())
+        {
+            dualCellIndexA = -1;
+        }
     }
 
     label dualCellIndexB = vB->index();
 
     if (!vB->internalOrBoundaryPoint() || vB->referred())
     {
-        dualCellIndexB = -1;
+        if (!vB->constrained())
+        {
+            dualCellIndexB = -1;
+        }
     }
 
     if (dualCellIndexA == -1 && dualCellIndexB == -1)
@@ -1065,6 +1107,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
             Foam::indexedVertexEnum::vtExternalFeaturePoint
         );
     }
+
+    //writeFixedPoints("fixedPointsStart.obj");
 }
 
 
@@ -1532,7 +1576,7 @@ void Foam::conformalVoronoiMesh::move()
     pointsToInsert.shrink();
 
     // Save displacements to file.
-    if (foamyHexMeshControls().objOutput() && runTime_.outputTime())
+    if (foamyHexMeshControls().objOutput() && time().outputTime())
     {
         Pout<< "Writing point displacement vectors to file." << endl;
         OFstream str("displacements_" + runTime_.timeName() + ".obj");
@@ -1567,38 +1611,62 @@ void Foam::conformalVoronoiMesh::move()
 
     Info<< nl << "Inserting displaced tessellation" << endl;
 
-    insertInternalPoints(pointsToInsert, true);
-
     reinsertFeaturePoints(true);
 
-    // Remove internal points that have been inserted outside the surface.
+    insertInternalPoints(pointsToInsert, true);
-//    label internalPtIsOutside = 0;
+
-//
+    {
-//    for
+        // Remove internal points that have been inserted outside the surface.
-//    (
+        label internalPtIsOutside = 0;
-//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+
-//        vit != finite_vertices_end();
+        autoPtr<OBJstream> str;
-//        ++vit
+
-//    )
+        if (foamyHexMeshControls().objOutput() && time().outputTime())
-//    {
+        {
-//        if (vit->internalPoint() && !vit->referred())
+            str.set
-//        {
+            (
-//            bool inside = geometryToConformTo_.inside
+                new OBJstream(time().path()/"internalPointsOutsideDomain.obj")
-//            (
+            );
-//                topoint(vit->point())
+        }
-//            );
+
-//
+        for
-//            if (!inside)
+        (
-//            {
+            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-//                remove(vit);
+            vit != finite_vertices_end();
-//                internalPtIsOutside++;
+            ++vit
-//            }
+        )
-//        }
+        {
-//    }
+            if
-//
+            (
-//    Info<< "    " << internalPtIsOutside
+                (vit->internalPoint() || vit->internalBoundaryPoint())
-//        << " internal points were inserted outside the domain. "
+             && !vit->referred()
-//        << "They have been removed." << endl;
+            )
+            {
+                bool inside = geometryToConformTo_.inside
+                (
+                    topoint(vit->point())
+                );
+
+                if (!inside)
+                {
+                    if
+                    (
+                        foamyHexMeshControls().objOutput()
+                     && time().outputTime()
+                    )
+                    {
+                        str().write(topoint(vit->point()));
+                    }
+                    remove(vit);
+                    internalPtIsOutside++;
+                }
+            }
+        }
+
+        Info<< "    " << internalPtIsOutside
+            << " internal points were inserted outside the domain. "
+            << "They have been removed." << endl;
+    }
 
     // Fix points that have not been significantly displaced
 //    for
@@ -1649,14 +1717,14 @@ void Foam::conformalVoronoiMesh::move()
     {
         writePoints
         (
-            "internalPoints_" + runTime_.timeName() + ".obj",
+            "internalPoints_" + time().timeName() + ".obj",
             Foam::indexedVertexEnum::vtInternal
         );
     }
 
-    if (foamyHexMeshControls().objOutput() && runTime_.outputTime())
+    if (foamyHexMeshControls().objOutput() && time().outputTime())
     {
-        writeBoundaryPoints("boundaryPoints_" + runTime_.timeName() + ".obj");
+        writeBoundaryPoints("boundaryPoints_" + time().timeName() + ".obj");
     }
 
     timeCheck("After conformToSurface");
@@ -1667,9 +1735,35 @@ void Foam::conformalVoronoiMesh::move()
     }
 
     // Write the intermediate mesh, do not filter the dual faces.
-    if (runTime_.outputTime())
+    if (time().outputTime())
     {
-        writeMesh(runTime_.timeName());
+        writeMesh(time().timeName());
+
+//        label cellI = 0;
+//        for
+//        (
+//            Finite_cells_iterator cit = finite_cells_begin();
+//            cit != finite_cells_end();
+//            ++cit
+//        )
+//        {
+//            if
+//            (
+//                !cit->hasFarPoint()
+//             && !is_infinite(cit)
+//            )
+//            {
+//                cit->cellIndex() = cellI++;
+//            }
+//        }
+//
+//        labelList vertexMap;
+//        labelList cellMap;
+//        autoPtr<fvMesh> tetMesh =
+//            createMesh("tetMesh", runTime_, vertexMap, cellMap);
+//
+//        tetMesh().write();
+        //writeFixedPoints("fixedPointsStart_" + runTime_.timeName() + ".obj");
     }
 
     updateSizesAndAlignments(pointsToInsert);

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

@@ -28,6 +28,7 @@ License
 #include "vectorTools.H"
 #include "indexedCellChecks.H"
 #include "IOmanip.H"
+#include "OBJstream.H"
 
 using namespace Foam::vectorTools;
 
@@ -519,7 +520,11 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                     }
                 }
             }
-            else if (vit->externalBoundaryPoint())
+            else if
+            (
+                vit->externalBoundaryPoint()
+             || (vit->externalBoundaryPoint() && vit->referred())
+            )
             {
                 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
 
@@ -618,6 +623,63 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             }
         }
 
+//        for
+//        (
+//            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+//            cit != finite_cells_end();
+//            ++cit
+//        )
+//        {
+//            if (cit->boundaryDualVertex() && !cit->parallelDualVertex())
+//            {
+//                const Foam::point& pt = cit->dual();
+//
+//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
+//                pointIndexHit surfHit;
+//                label hitSurface;
+//
+//                geometryToConformTo_.findSurfaceNearest
+//                (
+//                    pt,
+//                    sqr(targetCellSize(pt)),//surfaceSearchDistanceSqr(pt),
+//                    surfHit,
+//                    hitSurface
+//                );
+//
+//                if (!surfHit.hit())
+//                {
+//                    geometryToConformTo_.findSurfaceNearest
+//                    (
+//                        pt,
+//                        GREAT,
+//                        surfHit,
+//                        hitSurface
+//                    );
+//
+//                    if (surfHit.hit())
+//                    {
+//                        surfaceIntersections.append
+//                        (
+//                            pointIndexHitAndFeature(surfHit, hitSurface)
+//                        );
+//
+//                        addSurfaceAndEdgeHits
+//                        (
+//                            pt,
+//                            surfaceIntersections,
+//                            surfacePtReplaceDistCoeffSqr,
+//                            edgeSearchDistCoeffSqr,
+//                            surfaceHits,
+//                            featureEdgeHits,
+//                            surfaceToTreeShape,
+//                            edgeToTreeShape,
+//                            false
+//                        );
+//                    }
+//                }
+//            }
+//        }
+
 //        for
 //        (
 //            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -2281,6 +2343,7 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
             !vit->referred()
          && vit->boundaryPoint()
          && !vit->featurePoint()
+         && !vit->constrained()
         )
         {
             tempSurfaceVertices.append

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

@@ -334,23 +334,48 @@ inline void Foam::conformalVoronoiMesh::createPointPair
 {
     vector ppDistn = ppDist*n;
 
-    pts.append
+    const Foam::point internalPt = surfPt - ppDistn;
-    (
+    const Foam::point externalPt = surfPt + ppDistn;
-        Vb
-        (
-            surfPt - ppDistn,
-            Vb::vtInternalSurface
-        )
-    );
 
-    pts.append
+    if
     (
-        Vb
+        geometryToConformTo_.inside(internalPt)
+     && geometryToConformTo_.outside(externalPt)
+    )
+    {
+        pts.append
         (
-            surfPt + ppDistn,
+            Vb
-            Vb::vtExternalSurface
+            (
-        )
+                surfPt - ppDistn,
-    );
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                surfPt + ppDistn,
+                vertexCount() + pts.size(),
+                Vb::vtExternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+    }
+    else
+    {
+        Info<< "Warning: point pair not inside/outside" << nl
+            << "    surfPt = " << surfPt << nl
+            << "    internal = "
+            << internalPt << " " << geometryToConformTo_.inside(internalPt)
+            << nl
+            << "    external = "
+            << externalPt << " " << geometryToConformTo_.outside(externalPt)
+            << endl;
+    }
 }
 
 
@@ -386,12 +411,25 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
 
     pts.append
     (
-        Vb(surfPt - ppDistn, Vb::vtInternalSurface)
+        Vb
+        (
+            surfPt - ppDistn,
+            vertexCount() + 1 + pts.size(),
+            Vb::vtInternalSurface,
+            Pstream::myProcNo()
+        )
     );
 
     pts.append
     (
-        Vb(surfPt + ppDistn, Vb::vtInternalSurface)
+        Vb
+        (
+            surfPt + ppDistn,
+            vertexCount() + 1 + pts.size(),
+            Vb::vtInternalSurface,
+            Pstream::myProcNo()
+        )
+    );
     );
 }
 

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

@@ -1177,91 +1177,93 @@ void Foam::conformalVoronoiMesh::writeMesh
     }
 
 
-    volTensorField alignments
+//    volTensorField alignments
-    (
+//    (
-        IOobject
+//        IOobject
-        (
+//        (
-            "alignmentsField",
+//            "alignmentsField",
-            runTime_.timeName(),
+//            runTime_.timeName(),
-            runTime_,
+//            runTime_,
-            IOobject::NO_READ,
+//            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
+//            IOobject::AUTO_WRITE
-        ),
+//        ),
-        mesh,
+//        mesh,
-        tensor::zero
+//        tensor::zero
-    );
+//    );
+//
+//    forAll(mesh.cellCentres(), pI)
+//    {
+//        Vertex_handle nearV =
+//            nearest_vertex
+//            (
+//                toPoint<Point>(mesh.cellCentres()[pI])
+//            );
+//        alignments[pI] = nearV->alignment();
+//    }
+//    alignments.write();
+//
+//    {
+//        volVectorField alignmentx
+//        (
+//            IOobject
+//            (
+//                "alignmentsx",
+//                runTime_.timeName(),
+//                runTime_,
+//                IOobject::NO_READ,
+//                IOobject::AUTO_WRITE
+//            ),
+//            mesh,
+//            vector::zero
+//        );
+//        forAll(alignmentx, aI)
+//        {
+//            alignmentx[aI] = alignments[aI].x();
+//        }
+//        alignmentx.write();
+//    }
+//    {
+//        volVectorField alignmenty
+//        (
+//            IOobject
+//            (
+//                "alignmentsy",
+//                runTime_.timeName(),
+//                runTime_,
+//                IOobject::NO_READ,
+//                IOobject::AUTO_WRITE
+//            ),
+//            mesh,
+//            vector::zero
+//        );
+//        forAll(alignmenty, aI)
+//        {
+//            alignmenty[aI] = alignments[aI].y();
+//        }
+//        alignmenty.write();
+//    }
+//    {
+//        volVectorField alignmentz
+//        (
+//            IOobject
+//            (
+//                "alignmentsz",
+//                runTime_.timeName(),
+//                runTime_,
+//                IOobject::NO_READ,
+//                IOobject::AUTO_WRITE
+//            ),
+//            mesh,
+//            vector::zero
+//        );
+//        forAll(alignmentz, aI)
+//        {
+//            alignmentz[aI] = alignments[aI].z();
+//        }
+//        alignmentz.write();
+//    }
 
-    forAll(mesh.cellCentres(), pI)
-    {
-        Vertex_handle nearV =
-            nearest_vertex
-            (
-                toPoint<Point>(mesh.cellCentres()[pI])
-            );
-        alignments[pI] = nearV->alignment();
-    }
-    alignments.write();
 
-    {
-        volVectorField alignmentx
-        (
-            IOobject
-            (
-                "alignmentsx",
-                runTime_.timeName(),
-                runTime_,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            vector::zero
-        );
-        forAll(alignmentx, aI)
-        {
-            alignmentx[aI] = alignments[aI].x();
-        }
-        alignmentx.write();
-    }
-    {
-        volVectorField alignmenty
-        (
-            IOobject
-            (
-                "alignmentsy",
-                runTime_.timeName(),
-                runTime_,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            vector::zero
-        );
-        forAll(alignmenty, aI)
-        {
-            alignmenty[aI] = alignments[aI].y();
-        }
-        alignmenty.write();
-    }
-    {
-        volVectorField alignmentz
-        (
-            IOobject
-            (
-                "alignmentsz",
-                runTime_.timeName(),
-                runTime_,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            vector::zero
-        );
-        forAll(alignmentz, aI)
-        {
-            alignmentz[aI] = alignments[aI].z();
-        }
-        alignmentz.write();
-    }
     labelIOList boundaryIOPts
     (
         IOobject
@@ -1275,9 +1277,6 @@ void Foam::conformalVoronoiMesh::writeMesh
         boundaryPts
     );
 
-
-
-
     // Dump list of boundary points
     forAll(mesh.boundaryMesh(), patchI)
     {

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

@@ -115,6 +115,7 @@ class indexedCell
 
 public:
 
+    typedef typename Cb::Triangulation_data_structure    Tds;
     typedef typename Cb::Vertex_handle                   Vertex_handle;
     typedef typename Cb::Cell_handle                     Cell_handle;
 
@@ -183,6 +184,8 @@ public:
 
         inline bool hasInternalPoint() const;
 
+        inline bool hasConstrainedPoint() const;
+
         //- Does the Dual vertex form part of a processor patch
         inline bool parallelDualVertex() const;
 
@@ -227,6 +230,9 @@ public:
 
         inline bool potentialCoplanarCell() const;
 
+        inline bool featurePointExternalCell() const;
+        inline bool featurePointInternalCell() const;
+
 
     // Info
 

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

@@ -229,6 +229,19 @@ inline bool CGAL::indexedCell<Gt, Cb>::hasInternalPoint() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasConstrainedPoint() const
+{
+    return
+    (
+        this->vertex(0)->constrained()
+     || this->vertex(1)->constrained()
+     || this->vertex(2)->constrained()
+     || this->vertex(3)->constrained()
+    );
+}
+
+
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
 {
@@ -548,4 +561,128 @@ inline bool CGAL::indexedCell<Gt, Cb>::potentialCoplanarCell() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::featurePointExternalCell() const
+{
+    int featureVertex = -1;
+    for (int i = 0; i < 4; ++i)
+    {
+        if (this->vertex(i)->constrained())
+        {
+            featureVertex = i;
+        }
+    }
+
+    // Pick cell with a face attached to an infinite cell
+
+    if (featureVertex != -1)
+    {
+        Vertex_handle v1 =
+            this->vertex(Tds::vertex_triple_index(featureVertex, 0));
+        Vertex_handle v2 =
+            this->vertex(Tds::vertex_triple_index(featureVertex, 1));
+        Vertex_handle v3 =
+            this->vertex(Tds::vertex_triple_index(featureVertex, 2));
+
+        if (v1->internalBoundaryPoint())
+        {
+            if
+            (
+                v2->externalBoundaryPoint()
+             && v3->externalBoundaryPoint()
+            )
+            {
+                return true;
+            }
+        }
+        else if (v2->internalBoundaryPoint())
+        {
+            if
+            (
+                v1->externalBoundaryPoint()
+             && v3->externalBoundaryPoint()
+            )
+            {
+                return true;
+            }
+        }
+        else if (v3->internalBoundaryPoint())
+        {
+            if
+            (
+                v1->externalBoundaryPoint()
+             && v2->externalBoundaryPoint()
+            )
+            {
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::featurePointInternalCell() const
+{
+    int featureVertex = -1;
+    for (int i = 0; i < 4; ++i)
+    {
+        if (this->vertex(i)->constrained())
+        {
+            featureVertex = i;
+        }
+    }
+
+    // Pick cell with a face attached to an infinite cell
+
+    if (featureVertex != -1)
+    {
+        Vertex_handle v1 =
+            this->vertex(Tds::vertex_triple_index(featureVertex, 0));
+        Vertex_handle v2 =
+            this->vertex(Tds::vertex_triple_index(featureVertex, 1));
+        Vertex_handle v3 =
+            this->vertex(Tds::vertex_triple_index(featureVertex, 2));
+
+        if (v1->externalBoundaryPoint())
+        {
+            if
+            (
+                v2->internalBoundaryPoint()
+             && v3->internalBoundaryPoint()
+            )
+            {
+                return true;
+            }
+        }
+        else if (v2->externalBoundaryPoint())
+        {
+            if
+            (
+                v1->internalBoundaryPoint()
+             && v3->internalBoundaryPoint()
+            )
+            {
+                return true;
+            }
+        }
+        else if (v3->externalBoundaryPoint())
+        {
+            if
+            (
+                v1->internalBoundaryPoint()
+             && v2->internalBoundaryPoint()
+            )
+            {
+                return true;
+            }
+        }
+    }
+
+    return false;
+}
+
+
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //

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

@@ -248,6 +248,8 @@ public:
 
         inline bool externalBoundaryPoint() const;
 
+        inline bool constrained() const;
+
         //- Is the vertex fixed or movable
         inline bool fixed() const;
 

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

@@ -30,7 +30,7 @@ License
 
 template<>
 const char*
-Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 10>::names[] =
+Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>::names[] =
 {
     "Unassigned",
     "Internal",
@@ -41,10 +41,11 @@ Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 10>::names[] =
     "ExternalSurface",
     "ExternalFeatureEdge",
     "ExternalFeaturePoint",
-    "Far"
+    "Far",
+    "Constrained"
 };
 
-const Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 10>
+const Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>
 Foam::indexedVertexEnum::vertexTypeNames_;
 
 

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

@@ -58,7 +58,8 @@ public:
         vtExternalSurface       = 6,
         vtExternalFeatureEdge   = 7,
         vtExternalFeaturePoint  = 8,
-        vtFar                   = 9
+        vtFar                   = 9,
+        vtConstrained           = 10
     };
 
     enum vertexMotion
@@ -67,7 +68,7 @@ public:
         movable = 1
     };
 
-    static const Foam::NamedEnum<vertexType, 10> vertexTypeNames_;
+    static const Foam::NamedEnum<vertexType, 11> vertexTypeNames_;
 
     static const Foam::NamedEnum<vertexMotion, 2> vertexMotionNames_;
 

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

@@ -322,6 +322,13 @@ inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
 }
 
 
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::constrained() const
+{
+    return type_ == vtConstrained;
+}
+
+
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::featurePoint() const
 {

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

@@ -510,6 +510,21 @@ Foam::conformationSurfaces::conformationSurfaces
         << "Testing for locationInMesh " << locationInMesh_ << endl;
 
     hasBoundedVolume(referenceVolumeTypes_);
+
+    if (debug)
+    {
+        Info<< "Names = " << allGeometry_.names() << endl;
+        Info<< "Surfaces = " << surfaces_ << endl;
+        Info<< "AllGeom to Surfaces = " << allGeometryToSurfaces_ << endl;
+        Info<< "Baffle Surfaces = " << baffleSurfaces_ << endl;
+        Info<< "Patch names = " << patchNames_ << endl;
+        Info<< "Region Offset = " << regionOffset_ << endl;
+
+        forAll(features_, fI)
+        {
+            Info<< features_[fI].name() << endl;
+        }
+    }
 }
 
 
@@ -597,7 +612,34 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 
         if (!baffleSurfaces_[regionI])
         {
-            surface.getVolumeType(samplePts, surfaceVolumeTests[s]);
+//            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]);
+            }
         }
     }
 
@@ -664,8 +706,8 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 //            Info<< surface.name() << " = "
 //                << volumeType::names[surfaceVolumeTests[s][i]] << endl;
 
-            if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
             //if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
+            if (surfaceVolumeTests[s][i] != volumeType::INSIDE)
             {
                 insidePoint[i] = false;