Merge branch 'feature/cvMesh' of /home/dm4/OpenFOAM/OpenFOAM-dev into feature/cvMesh

applications/utilities/mesh/generation/Allwmake

@@ -12,7 +12,7 @@ wmake snappyHexMesh
 if [ -d "$CGAL_ARCH_PATH" ]
 then
     foamyHexMesh/Allwmake
-    cv2DMesh/Allwmake
+    foamyHex2DMesh/Allwmake
 fi
 
 # ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/foamyHex2DMesh/Make/files

@@ -7,6 +7,6 @@ insertSurfaceNearPointPairs.C
 insertBoundaryConformPointPairs.C
 CV2DIO.C
 shortEdgeFilter2D.C
-cv2DMesh.C
+foamyHex2DMesh.C
 
-EXE = $(FOAM_APPBIN)/cv2DMesh
+EXE = $(FOAM_APPBIN)/foamyHex2DMesh

applications/utilities/mesh/generation/foamyHex2DMesh/Make/options

@@ -12,7 +12,7 @@ EXE_INC = \
     -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
     -I../cvMesh/vectorTools \
     -IconformalVoronoi2DMesh/lnInclude \
-    -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
+    -I$(FOAM_APP)/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/utilities/mesh/generation/foamyHex2DMesh/foamyHex2DMesh.C

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cv2DMesh
+    foamyHex2DMesh
 
 Description
     Conformal-Voronoi 2D extruding automatic mesher with grid or read

applications/utilities/mesh/generation/foamyHex2DMesh/foamyHex2DMeshDict

@@ -5,21 +5,14 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-
 FoamFile
 {
     version         2.0;
     format          ascii;
-
-    root            "";
-    case            "";
-    instance        "";
-    local           "";
-
     class           dictionary;
-    object          cv2DMeshDict;
+    location        "system";
+    object          foamyHex2DMeshDict;
 }
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
@@ -47,42 +40,41 @@ geometry
 
 surfaceConformation
 {
-    locationInMesh (-2.8 0.7 0.5);
+    locationInMesh              (-2.8 0.7 0.5);
 
-    pointPairDistanceCoeff    0.005;
+    pointPairDistanceCoeff      0.005;
 
-    minEdgeLenCoeff           0.005;
+    minEdgeLenCoeff             0.005;
 
-    maxNotchLenCoeff          0.003;
+    maxNotchLenCoeff            0.003;
 
-    minNearPointDistCoeff     0.0025;
+    minNearPointDistCoeff       0.0025;
 
-    maxQuadAngle 125;
+    maxQuadAngle                125;
 
     // Insert near-boundary point mirror or point-pairs
     insertSurfaceNearestPointPairs yes;
 
     // Mirror near-boundary points rather than insert point-pairs
-    mirrorPoints   no;
+    mirrorPoints                no;
 
     // Insert point-pairs vor dual-cell vertices very near the surface
     insertSurfaceNearPointPairs yes;
 
     // Maximum number of iterations used in boundaryConform.
-    maxBoundaryConformingIter 5;
+    maxBoundaryConformingIter   5;
 
     geometryToConformTo
     {
         laurence_clean_preciser
         {
-            featureMethod extendedFeatureEdgeMesh;
+            featureMethod           extendedFeatureEdgeMesh;
             extendedFeatureEdgeMesh "laurence_clean_preciser.extendedFeatureEdgeMesh";
         }
     }
 
     additionalFeatures
-    {
-    }
+    {}
 
     // Choose if to randomise the initial grid created by insertGrid.
     randomiseInitialGrid yes;
@@ -148,27 +140,26 @@ motionControl
 
     adaptiveLinearCoeffs
     {
-        relaxationStart 0.5;
+        relaxationStart 1.0;
         relaxationEnd   0.0;
     }
 
-    objOutput                   no;
+    objOutput           no;
 
     // Near-wall region where cells are aligned with the wall specified as a number
     // of cell layers
     nearWallAlignedDist 3;
-
 }
 
 shortEdgeFilter
 {
     // Factor to multiply the average of a face's edge lengths by.
     // If an edge of that face is smaller than that value then delete it.
-    shortEdgeFilterFactor     0.2;
+    shortEdgeFilterFactor           0.2;
 
     // Weighting for the lengths of edges that are attached to the boundaries.
     // Used when calculating the length of an edge. Default 2.0.
-    edgeAttachedToBoundaryFactor 2.0;
+    edgeAttachedToBoundaryFactor    2.0;
 }
 
 extrusion
@@ -206,3 +197,6 @@ extrusion
 
     thickness       0.1;
 }
+
+
+// ************************************************************************* //

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

@@ -64,10 +64,40 @@ void Foam::DelaunayMesh<Triangulation>::reset()
 {
     Info<< "Clearing triangulation" << endl;
 
+    DynamicList<Vb> vertices;
+
+    for
+    (
+        Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
+        vit != Triangulation::finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->fixed())
+        {
+            vertices.append
+            (
+                Vb
+                (
+                    vit->point(),
+                    vit->index(),
+                    vit->type(),
+                    vit->procIndex()
+                )
+            );
+
+            vertices.last().fixed() = vit->fixed();
+        }
+    }
+
     this->clear();
 
     resetVertexCount();
     resetCellCount();
+
+    insertPoints(vertices);
+
+    Info<< "Inserted " << vertexCount() << " fixed points" << endl;
 }
 
 

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

@@ -354,8 +354,8 @@ Foam::DelaunayMesh<Triangulation>::createMesh
     List<DynamicList<face> > patchFaces(1, DynamicList<face>());
     List<DynamicList<label> > patchOwners(1, DynamicList<label>());
 
-    vertexMap.setSize(vertexCount());
-    cellMap.setSize(Triangulation::number_of_finite_cells());
+    vertexMap.setSize(vertexCount(), -1);
+    cellMap.setSize(Triangulation::number_of_finite_cells(), -1);
 
     // Calculate pts and a map of point index to location in pts.
     label vertI = 0;

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

@@ -19,6 +19,7 @@ EXE_INC = \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -IPrintTable \
@@ -30,4 +31,5 @@ LIB_LIBS = \
     -lfileFormats \
     -ltriSurface \
     -ldynamicMesh \
+    -lsurfMesh \
     -lsampling

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

@@ -62,11 +62,6 @@ bool Foam::cellSizeAndAlignmentControls::evalCellSizeFunctions
                 const searchableSurfaceControl& sSC =
                     refCast<const searchableSurfaceControl>(cSF);
 
-                if (debug)
-                {
-                    Info<< "size function " << sSC.name() << endl;
-                }
-
                 anyFunctionFound = sSC.cellSize(pt, minSize, previousPriority);
 
                 if (previousPriority > maxPriority)
@@ -127,6 +122,21 @@ Foam::cellSizeAndAlignmentControls::cellSizeAndAlignmentControls
 
         functionI++;
     }
+
+    // Sort controlFunctions_ by maxPriority
+    SortableList<label> functionPriorities(functionI);
+
+    forAll(controlFunctions_, funcI)
+    {
+        functionPriorities[funcI] = controlFunctions_[funcI].maxPriority();
+    }
+
+    functionPriorities.reverseSort();
+
+    labelList invertedFunctionPriorities =
+        invert(functionPriorities.size(), functionPriorities.indices());
+
+    controlFunctions_.reorder(invertedFunctionPriorities);
 }
 
 

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

@@ -288,6 +288,7 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
         cellSizeFunctions_.transfer(regionCellSizeFunctions);
     }
 
+
     forAll(cellSizeFunctions_, funcI)
     {
         const label funcPriority = cellSizeFunctions_[funcI].priority();
@@ -298,6 +299,21 @@ Foam::searchableSurfaceControl::searchableSurfaceControl
         }
     }
 
+    // Sort controlFunctions_ by maxPriority
+    SortableList<label> functionPriorities(cellSizeFunctions_.size());
+
+    forAll(cellSizeFunctions_, funcI)
+    {
+        functionPriorities[funcI] = cellSizeFunctions_[funcI].priority();
+    }
+
+    functionPriorities.reverseSort();
+
+    labelList invertedFunctionPriorities =
+        invert(functionPriorities.size(), functionPriorities.indices());
+
+    cellSizeFunctions_.reorder(invertedFunctionPriorities);
+
     Info<< nl << "There are " << cellSizeFunctions_.size()
         << " region control functions" << endl;
 }
@@ -340,6 +356,8 @@ void Foam::searchableSurfaceControl::initialVertices
             infoFeature
         );
 
+        scalar limitedCellSize = GREAT;
+
         autoPtr<triad> pointAlignment;
 
         if (info.hit())
@@ -398,6 +416,26 @@ void Foam::searchableSurfaceControl::initialVertices
                 searchableSurface_.getNormal(infoList, normals);
 
                 pointAlignment.set(new triad(normals[0]));
+
+                // Limit cell size
+                const vector vN =
+                    infoList[0].hitPoint()
+                  - 2.0*normals[0]*defaultCellSize_;
+
+                List<pointIndexHit> intersectionList;
+                searchableSurface_.findLineAny
+                (
+                    ptField,
+                    pointField(1, vN),
+                    intersectionList
+                );
+
+                if (intersectionList[0].hit())
+                {
+                    scalar dist = mag(intersectionList[0].hitPoint() - pts[pI]);
+
+                    //limitedCellSize = dist/2.0;
+                }
             }
         }
 
@@ -412,6 +450,8 @@ void Foam::searchableSurfaceControl::initialVertices
                 << abort(FatalError);
         }
 
+        sizes[pI] = min(limitedCellSize, sizes[pI]);
+
         alignments[pI] = pointAlignment();
     }
 }

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

@@ -277,9 +277,14 @@ 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
@@ -387,7 +392,7 @@ void Foam::controlMeshRefinement::initialMeshPopulation
 
             if (debug)
             {
-                Info<< "    size difference = " << sizeDiff
+                Info<< "    size difference = " << sizeDiff << nl
                     << ", alignment difference = " << alignmentDiff << endl;
             }
 
@@ -568,30 +573,35 @@ void Foam::controlMeshRefinement::initialMeshPopulation
             {
                 // Check the priority
 
-                cellShapeControlMesh::Cell_handle ch =
-                    mesh_.locate(toPoint<cellShapeControlMesh::Point>(pt));
+//                cellShapeControlMesh::Cell_handle ch =
+//                    mesh_.locate(toPoint<cellShapeControlMesh::Point>(pt));
 
-                const label newPtPriority = controlFunction.maxPriority();
+//                if (mesh_.is_infinite(ch))
+//                {
+//                    continue;
+//                }
 
-                label highestPriority = -1;
-                for (label cI = 0; cI < 4; ++cI)
-                {
-                    if (mesh_.is_infinite(ch->vertex(cI)))
-                    {
-                        continue;
-                    }
+//                const label newPtPriority = controlFunction.maxPriority();
 
-                    const label vertPriority =
-                        priorityMap[ch->vertex(cI)->index()];
+//                label highestPriority = -1;
+//                for (label cI = 0; cI < 4; ++cI)
+//                {
+//                    if (mesh_.is_infinite(ch->vertex(cI)))
+//                    {
+//                        continue;
+//                    }
 
-                    if (vertPriority > highestPriority)
-                    {
-                        highestPriority = vertPriority;
-                    }
-                }
+//                    const label vertPriority =
+//                        priorityMap[ch->vertex(cI)->index()];
 
-                if (newPtPriority >= highestPriority)
-                {
+//                    if (vertPriority > highestPriority)
+//                    {
+//                        highestPriority = vertPriority;
+//                    }
+//                }
+
+//                if (newPtPriority >= highestPriority)
+//                {
                     const label oldSize = mesh_.vertexCount();
 
                     cellShapeControlMesh::Vertex_handle insertedVert =
@@ -603,15 +613,15 @@ void Foam::controlMeshRefinement::initialMeshPopulation
                             Vb::vtInternal
                         );
 
-                    if (oldSize == mesh_.vertexCount() - 1)
-                    {
-                        priorityMap.insert
-                        (
-                            insertedVert->index(),
-                            newPtPriority
-                        );
-                    }
-                }
+//                    if (oldSize == mesh_.vertexCount() - 1)
+//                    {
+//                        priorityMap.insert
+//                        (
+//                            insertedVert->index(),
+//                            newPtPriority
+//                        );
+//                    }
+//                }
             }
         }
 

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

@@ -27,6 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "triSurfaceMesh.H"
 #include "vtkSurfaceWriter.H"
+#include "primitivePatchInterpolation.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -68,9 +69,9 @@ void Foam::automatic::smoothField(triSurfaceScalarField& surf)
                 const scalar faceValue = surf[faceLabel];
                 const scalar distance = mag(faceCentre - fC);
 
-                newValue += faceValue/distance;
+                newValue += faceValue/(distance + SMALL);
 
-                totalDist += 1.0/distance;
+                totalDist += 1.0/(distance + SMALL);
 
                 if (value < faceValue)
                 {
@@ -81,7 +82,7 @@ void Foam::automatic::smoothField(triSurfaceScalarField& surf)
             // Do not smooth out the peak values
             if (nFaces == faceFaces.size())
             {
-                //continue;
+                continue;
             }
 
             surf[sI] = newValue/totalDist;
@@ -107,7 +108,6 @@ Foam::automatic::automatic
         defaultCellSize
     ),
     coeffsDict_(cellSizeCalcTypeDict.subDict(typeName + "Coeffs")),
-    surface_(surface),
     surfaceName_(surface.searchableSurface::name()),
     readCurvature_(Switch(coeffsDict_.lookup("curvature"))),
     curvatureFile_(coeffsDict_.lookup("curvatureFile")),
@@ -128,31 +128,36 @@ Foam::automatic::automatic
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::triSurfaceScalarField Foam::automatic::load()
+Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
 {
-    Info<< indent << "Calculating cell size on surface: "
-        << surfaceName_ << endl;
+    Info<< indent
+        << "Calculating cell size on surface: " << surfaceName_ << endl;
 
-    triSurfaceScalarField surfaceCellSize
+    tmp<triSurfacePointScalarField> tPointCellSize
     (
-        IOobject
+        new triSurfacePointScalarField
         (
-            surfaceName_ + ".cellSize",
-            surface_.searchableSurface::time().constant(),
-            "triSurface",
-            surface_.searchableSurface::time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        surface_,
-        dimLength,
-        scalarField(surface_.size(), maximumCellSize_)
+            IOobject
+            (
+                surfaceName_ + ".cellSize",
+                surface_.searchableSurface::time().constant(),
+                "triSurface",
+                surface_.searchableSurface::time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surface_,
+            dimLength,
+            scalarField(surface_.nPoints(), maximumCellSize_)
+        )
     );
 
+    triSurfacePointScalarField& pointCellSize = tPointCellSize();
+
     if (readCurvature_)
     {
-        Info<< indent << "Reading curvature         : "
-            << curvatureFile_ << endl;
+        Info<< indent
+            << "Reading curvature         : " << curvatureFile_ << endl;
 
         triSurfacePointScalarField curvature
         (
@@ -170,42 +175,33 @@ Foam::triSurfaceScalarField Foam::automatic::load()
             true
         );
 
-        const List<labelledTri>& localFaces = surface_.localFaces();
-        const labelList& meshPoints = surface_.meshPoints();
-
-        forAll(surfaceCellSize, fI)
+        forAll(pointCellSize, pI)
         {
-            const labelList& facePoints = localFaces[fI].triFaceFace();
-
-            scalar interpolatedCurvatureToFace = 0.0;
-
-            forAll(facePoints, fpI)
-            {
-                interpolatedCurvatureToFace
-                    += curvature[meshPoints[facePoints[fpI]]];
-            }
-
-            interpolatedCurvatureToFace /= facePoints.size();
-
-            surfaceCellSize[fI] =
+            pointCellSize[pI] =
                 min
                 (
                     1.0
                    /max
                     (
-                        (1.0/curvatureCellSizeCoeff_)
-                       *interpolatedCurvatureToFace,
+                        (1.0/curvatureCellSizeCoeff_)*mag(curvature[pI]),
                         1.0/maximumCellSize_
                     ),
-                    surfaceCellSize[fI]
+                    pointCellSize[pI]
                 );
         }
     }
 
+    PrimitivePatchInterpolation
+    <
+        PrimitivePatch<labelledTri, ::Foam::List, pointField, point>
+    > patchInterpolate(surface_);
+
+    const Map<label>& meshPointMap = surface_.meshPointMap();
+
     if (readInternalCloseness_)
     {
-        Info<< indent << "Reading internal closeness: "
-            << internalClosenessFile_ << endl;
+        Info<< indent
+            << "Reading internal closeness: " << internalClosenessFile_ << endl;
 
         triSurfaceScalarField internalCloseness
         (
@@ -223,21 +219,24 @@ Foam::triSurfaceScalarField Foam::automatic::load()
             true
         );
 
-        forAll(surfaceCellSize, fI)
+        scalarField internalClosenessPointField =
+            patchInterpolate.faceToPointInterpolate(internalCloseness);
+
+        forAll(pointCellSize, pI)
         {
-            surfaceCellSize[fI] =
+            pointCellSize[pI] =
                 min
                 (
-                    internalCloseness[fI],
-                    surfaceCellSize[fI]
+                    internalClosenessPointField[meshPointMap[pI]],
+                    pointCellSize[pI]
                 );
         }
     }
 
     if (readFeatureProximity_)
     {
-        Info<< indent << "Reading feature proximity : "
-            << featureProximityFile_ << endl;
+        Info<< indent
+            << "Reading feature proximity : " << featureProximityFile_ << endl;
 
         triSurfaceScalarField featureProximity
         (
@@ -255,20 +254,23 @@ Foam::triSurfaceScalarField Foam::automatic::load()
             true
         );
 
-        forAll(surfaceCellSize, fI)
+        scalarField featureProximityPointField =
+            patchInterpolate.faceToPointInterpolate(featureProximity);
+
+        forAll(pointCellSize, pI)
         {
-            surfaceCellSize[fI] =
+            pointCellSize[pI] =
                 min
                 (
-                    featureProximity[fI],
-                    surfaceCellSize[fI]
+                    featureProximityPointField[meshPointMap[pI]],
+                    pointCellSize[pI]
                 );
         }
     }
 
-    smoothField(surfaceCellSize);
+    //smoothField(surfaceCellSize);
 
-    surfaceCellSize.write();
+    pointCellSize.write();
 
     if (debug)
     {
@@ -282,17 +284,17 @@ Foam::triSurfaceScalarField Foam::automatic::load()
         vtkSurfaceWriter().write
         (
             surface_.searchableSurface::time().constant()/"triSurface",
-            surfaceName_,
+            surfaceName_.lessExt().name(),
             surface_.points(),
             faces,
             "cellSize",
-            surfaceCellSize,
-            false,
+            pointCellSize,
+            true,
             true
         );
     }
 
-    return surfaceCellSize;
+    return tPointCellSize;
 }
 
 

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

@@ -62,11 +62,8 @@ private:
         //- Dictionary of coefficients for automatic cell sizing
         const dictionary& coeffsDict_;
 
-        //- Reference to the triSurfaceMesh
-        const triSurfaceMesh& surface_;
-
         //- Name of the surface. Used to write the cell size field
-        const word surfaceName_;
+        const fileName surfaceName_;
 
         const Switch readCurvature_;
         const word curvatureFile_;
@@ -115,7 +112,7 @@ public:
     // Member Functions
 
         //- Load the cell size field
-        virtual triSurfaceScalarField load();
+        virtual tmp<triSurfacePointScalarField> load();
 };
 
 

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

@@ -46,6 +46,7 @@ Foam::cellSizeCalculationType::cellSizeCalculationType
 )
 :
     cellSizeCalculationTypeDict_(cellSizeCalculationTypeDict),
+    surface_(surface),
     defaultCellSize_(defaultCellSize)
 {}
 

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

@@ -35,7 +35,6 @@ SourceFiles
 #ifndef cellSizeCalculationType_H
 #define cellSizeCalculationType_H
 
-#include "searchableSurface.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 #include "triSurfaceFields.H"
@@ -58,6 +57,9 @@ protected:
 
         const dictionary& cellSizeCalculationTypeDict_;
 
+        //- Reference to the triSurfaceMesh
+        const triSurfaceMesh& surface_;
+
         const scalar& defaultCellSize_;
 
 
@@ -124,7 +126,7 @@ public:
     // Member Functions
 
         //- Load the cell size
-        virtual triSurfaceScalarField load() = 0;
+        virtual tmp<triSurfacePointScalarField> load() = 0;
 };
 
 

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

@@ -58,7 +58,6 @@ Foam::fieldFromFile::fieldFromFile
         surface,
         defaultCellSize
     ),
-    surface_(surface),
     fileName_
     (
         cellSizeCalcTypeDict.subDict(typeName + "Coeffs").lookup("fieldFile")
@@ -68,27 +67,30 @@ Foam::fieldFromFile::fieldFromFile
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::triSurfaceScalarField Foam::fieldFromFile::load()
+Foam::tmp<Foam::triSurfacePointScalarField> Foam::fieldFromFile::load()
 {
     Info<< indent << "Loading: " << fileName_ << endl;
 
-    triSurfaceScalarField surfaceCellSize
+    tmp<triSurfacePointScalarField> pointCellSize
     (
-        IOobject
+        new triSurfacePointScalarField
         (
-            fileName_,
-            surface_.searchableSurface::time().constant(),
-            "triSurface",
-            surface_.searchableSurface::time(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        surface_,
-        dimLength,
-        true
+            IOobject
+            (
+                fileName_,
+                surface_.searchableSurface::time().constant(),
+                "triSurface",
+                surface_.searchableSurface::time(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            surface_,
+            dimLength,
+            true
+        )
     );
 
-    return surfaceCellSize;
+    return pointCellSize;
 }
 
 

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

@@ -58,9 +58,6 @@ private:
 
     // Private data
 
-        //- Reference to the triSurfaceMesh
-        const triSurfaceMesh& surface_;
-
         //- Name of the triSurfaceScalarField file to load in. Must be in
         //  constant/triSurface
         const word fileName_;
@@ -91,7 +88,7 @@ public:
     // Member Functions
 
         //- Load the cell size field
-        virtual triSurfaceScalarField load();
+        virtual tmp<triSurfacePointScalarField> load();
 };
 
 

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

@@ -60,21 +60,6 @@ Foam::nonUniformField::nonUniformField
         defaultCellSize
     ),
     surfaceTriMesh_(refCast<const triSurfaceMesh>(surface)),
-    surfaceCellSize_
-    (
-        IOobject
-        (
-            "surfaceCellSize.cellSize",
-            surface.time().constant(),
-            "triSurface",
-            surface.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        surfaceTriMesh_,
-        dimLength,
-        false
-    ),
     cellSizeCalculationType_
     (
         cellSizeCalculationType::New
@@ -84,17 +69,30 @@ Foam::nonUniformField::nonUniformField
             defaultCellSize
         )
     ),
-    pointCellSize_(),
-    patch_()
+    pointCellSize_
+    (
+        IOobject
+        (
+            "pointCellSize.cellSize",
+            surfaceTriMesh_.searchableSurface::time().constant(),
+            "triSurface",
+            surfaceTriMesh_.searchableSurface::time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surfaceTriMesh_,
+        dimLength,
+        false
+    )
 {
     Info<< incrIndent;
 
-    surfaceCellSize_ = cellSizeCalculationType_().load();
+    pointCellSize_ = cellSizeCalculationType_().load();
 
     Info<< indent << "Cell size field statistics:" << nl
-        << indent << "    Minimum: " << min(surfaceCellSize_).value() << nl
-        << indent << "    Average: " << average(surfaceCellSize_).value() << nl
-        << indent << "    Maximum: " << max(surfaceCellSize_).value() << endl;
+        << indent << "    Minimum: " << min(pointCellSize_).value() << nl
+        << indent << "    Average: " << average(pointCellSize_).value() << nl
+        << indent << "    Maximum: " << max(pointCellSize_).value() << endl;
 
     Info<< decrIndent;
 }
@@ -102,31 +100,12 @@ Foam::nonUniformField::nonUniformField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::scalar& Foam::nonUniformField::surfaceSize(const label index) const
-{
-    return surfaceCellSize_[index];
-}
-
-
-const Foam::scalar& Foam::nonUniformField::refineSurfaceSize(const label index)
-{
-    surfaceCellSize_[index] *= refinementFactor_;
-
-    return surfaceCellSize_[index];
-}
-
-
 Foam::scalar Foam::nonUniformField::interpolate
 (
     const point& pt,
     const label index
 ) const
 {
-    if (patch_.empty())
-    {
-        recalculateInterpolation();
-    }
-
     const face& faceHitByPt = surfaceTriMesh_.triSurface::operator[](index);
 
     const pointField& pts = surfaceTriMesh_.points();
@@ -143,34 +122,9 @@ Foam::scalar Foam::nonUniformField::interpolate
 
     tri.barycentric(pt, bary);
 
-    return pointCellSize_()[pMap[faceHitByPt[0]]]*bary[0]
-         + pointCellSize_()[pMap[faceHitByPt[1]]]*bary[1]
-         + pointCellSize_()[pMap[faceHitByPt[2]]]*bary[2];
-}
-
-
-void Foam::nonUniformField::recalculateInterpolation() const
-{
-    patch_.reset(new primitivePatchInterpolation(surfaceTriMesh_));
-
-    pointCellSize_.reset
-    (
-        new triSurfacePointScalarField
-        (
-            IOobject
-            (
-                "pointCellSize.cellSize",
-                surfaceTriMesh_.searchableSurface::time().constant(),
-                "triSurface",
-                surfaceTriMesh_.searchableSurface::time(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            surfaceTriMesh_,
-            dimLength,
-            patch_().faceToPointInterpolate(surfaceCellSize_)
-        )
-    );
+    return pointCellSize_[pMap[faceHitByPt[0]]]*bary[0]
+         + pointCellSize_[pMap[faceHitByPt[1]]]*bary[1]
+         + pointCellSize_[pMap[faceHitByPt[2]]]*bary[2];
 }
 
 

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

@@ -70,13 +70,9 @@ protected:
 
         const triSurfaceMesh& surfaceTriMesh_;
 
-        triSurfaceScalarField surfaceCellSize_;
-
         autoPtr<cellSizeCalculationType> cellSizeCalculationType_;
 
-        mutable autoPtr<triSurfacePointScalarField> pointCellSize_;
-
-        mutable autoPtr<primitivePatchInterpolation> patch_;
+        triSurfacePointScalarField pointCellSize_;
 
 
 public:
@@ -103,13 +99,15 @@ public:
 
     // Member Functions
 
-        virtual const scalar& surfaceSize(const label index) const;
+        // Query
 
-        virtual const scalar& refineSurfaceSize(const label index);
-
-        virtual scalar interpolate(const point& pt, const label index) const;
-
-        virtual void recalculateInterpolation() const;
+            //- Return the interpolated cell size for a point in the given
+            //  surface triangle
+            virtual scalar interpolate
+            (
+                const point& pt,
+                const label index
+            ) const;
 };
 
 

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

@@ -60,7 +60,7 @@ void Foam::conformalVoronoiMesh::cellSizeMeshOverlapsBackground() const
         ++vit
     )
     {
-        if (vit->internalOrBoundaryPoint())
+        if (vit->internalOrBoundaryPoint() && !vit->referred())
         {
             pts.append(topoint(vit->point()));
         }
@@ -141,7 +141,7 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
     Info<< "    " << nInserted << " points inserted"
         << ", failed to insert " << nPoints - nInserted
         << " ("
-        << 100.0*(nPoints - nInserted)/nInserted
+        << 100.0*(nPoints - nInserted)/(nInserted + SMALL)
         << " %)"<< endl;
 
     for
@@ -330,6 +330,28 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
     // Assume that the initial points method made the correct decision for
     // which processor each point should be on, so give distribute = false
     insertInternalPoints(initPts);
+
+    if (initialPointsMethod_->fixInitialPoints())
+    {
+        for
+        (
+            Finite_vertices_iterator vit = finite_vertices_begin();
+            vit != finite_vertices_end();
+            ++vit
+        )
+        {
+            vit->fixed() = true;
+        }
+    }
+
+    if (foamyHexMeshControls().objOutput())
+    {
+        writePoints
+        (
+            "initialPoints.obj",
+            Foam::indexedVertexEnum::vtInternal
+        );
+    }
 }
 
 
@@ -1187,8 +1209,8 @@ void Foam::conformalVoronoiMesh::move()
 
                 if
                 (
-                    vA->internalPoint() && !vA->referred()
-                 && vB->internalPoint() && !vB->referred()
+                    vA->internalPoint() && !vA->referred() && !vA->fixed()
+                 && vB->internalPoint() && !vB->referred() && !vB->fixed()
                 )
                 {
                     // Only insert a point at the midpoint of
@@ -1209,12 +1231,12 @@ void Foam::conformalVoronoiMesh::move()
                     }
                 }
 
-                if (vA->internalPoint() && !vA->referred())
+                if (vA->internalPoint() && !vA->referred() && !vA->fixed())
                 {
                     pointToBeRetained[vA->index()] = false;
                 }
 
-                if (vB->internalPoint() && !vB->referred())
+                if (vB->internalPoint() && !vB->referred() && !vB->fixed())
                 {
                     pointToBeRetained[vB->index()] = false;
                 }
@@ -1358,8 +1380,12 @@ void Foam::conformalVoronoiMesh::move()
                         // Point removal
                         if
                         (
-                            vA->internalPoint() && !vA->referred()
-                         && vB->internalPoint() && !vB->referred()
+                            vA->internalPoint()
+                         && !vA->referred()
+                         && !vA->fixed()
+                         && vB->internalPoint()
+                         && !vB->referred()
+                         && !vB->fixed()
                         )
                         {
                             // Only insert a point at the midpoint of
@@ -1379,26 +1405,60 @@ void Foam::conformalVoronoiMesh::move()
                             }
                         }
 
-                        if (vA->internalPoint() && !vA->referred())
+                        if
+                        (
+                            vA->internalPoint()
+                         && !vA->referred()
+                         && !vA->fixed()
+                        )
                         {
                             pointToBeRetained[vA->index()] = false;
                         }
 
-                        if (vB->internalPoint() && !vB->referred())
+                        if
+                        (
+                            vB->internalPoint()
+                         && !vB->referred()
+                         && !vB->fixed()
+                        )
                         {
                             pointToBeRetained[vB->index()] = false;
                         }
                     }
                     else
                     {
-                        if (vA->internalPoint() && !vA->referred())
+                        if
+                        (
+                            vA->internalPoint()
+                         && !vA->referred()
+                         && !vA->fixed()
+                        )
                         {
-                            displacementAccumulator[vA->index()] += delta;
+                            if (vB->fixed())
+                            {
+                                displacementAccumulator[vA->index()] += 2*delta;
+                            }
+                            else
+                            {
+                                displacementAccumulator[vA->index()] += delta;
+                            }
                         }
 
-                        if (vB->internalPoint() && !vB->referred())
+                        if
+                        (
+                            vB->internalPoint()
+                         && !vB->referred()
+                         && !vB->fixed()
+                        )
                         {
-                            displacementAccumulator[vB->index()] -= delta;
+                            if (vA->fixed())
+                            {
+                                displacementAccumulator[vB->index()] -= 2*delta;
+                            }
+                            else
+                            {
+                                displacementAccumulator[vB->index()] -= delta;
+                            }
                         }
                     }
                 }
@@ -1416,7 +1476,7 @@ void Foam::conformalVoronoiMesh::move()
         ++vit
     )
     {
-        if (vit->internalPoint() && !vit->referred())
+        if (vit->internalPoint() && !vit->referred() && !vit->fixed())
         {
             if (pointToBeRetained[vit->index()] == true)
             {
@@ -1443,7 +1503,7 @@ void Foam::conformalVoronoiMesh::move()
         ++vit
     )
     {
-        if (vit->internalPoint() && !vit->referred())
+        if (vit->internalPoint() && !vit->referred() && !vit->fixed())
         {
             if (pointToBeRetained[vit->index()] == true)
             {

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

@@ -181,6 +181,30 @@ void Foam::conformalVoronoiMesh::writePoints
 }
 
 
+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

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

@@ -118,7 +118,9 @@ Foam::Ostream& Foam::operator<<
         << static_cast<int>(p.type())
         << p.procIndex()
         << p.alignment()
-        << p.targetCellSize();
+        << p.targetCellSize()
+        << token::SPACE
+        << static_cast<int>(p.fixed());
 
     return os;
 }
@@ -136,13 +138,16 @@ Foam::Istream& Foam::operator>>
 
     int type;
     is  >> type;
-
     p.type() = static_cast<Foam::indexedVertexEnum::vertexType>(type);
 
     is  >> p.procIndex()
         >> p.alignment()
         >> p.targetCellSize();
 
+    int fixed;
+    is  >> fixed;
+    p.fixed() = static_cast<bool>(fixed);
+
     return is;
 }
 
@@ -163,6 +168,13 @@ Foam::Ostream& Foam::operator<<
         CGAL::to_double(iv.point().z())
     );
 
+    string fixed
+    (
+        iv.vertexFixed_
+      ? string(" fixed, ")
+      : string(" free, ")
+    );
+
     string referred
     (
         Pstream::myProcNo() == iv.processor_
@@ -175,6 +187,7 @@ Foam::Ostream& Foam::operator<<
         << " at:" << pt
         << " size:" << iv.targetCellSize_
         << " alignment:" << iv.alignment_
+        << fixed
         << referred.c_str()
         << endl;
 

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

@@ -124,7 +124,7 @@ class indexedVertex
         Foam::scalar targetCellSize_;
 
         //- Specify whether the vertex is fixed or movable.
-//        bool vertexFixed_;
+        bool vertexFixed_;
 
 
 public:
@@ -159,6 +159,14 @@ public:
             int processor
         );
 
+        inline indexedVertex
+        (
+            const Foam::point& p,
+            int index,
+            vertexType type,
+            int processor
+        );
+
         inline indexedVertex(const Point& p, Cell_handle f);
 
         inline indexedVertex(Cell_handle f);
@@ -240,11 +248,11 @@ public:
 
         inline bool externalBoundaryPoint() const;
 
-//        //- Is the vertex fixed or movable
-//        inline bool isVertexFixed() const;
-//
-//        //- Fix the vertex so that it can't be moved
-//        inline void setVertexFixed();
+        //- Is the vertex fixed or movable
+        inline bool fixed() const;
+
+        //- Fix the vertex so that it can't be moved
+        inline bool& fixed();
 
         inline indexedVertex& operator=(const indexedVertex& rhs)
         {
@@ -255,6 +263,7 @@ public:
             this->processor_ = rhs.procIndex();
             this->alignment_ = rhs.alignment();
             this->targetCellSize_ = rhs.targetCellSize();
+            this->vertexFixed_ = rhs.fixed();
 
             return *this;
         }
@@ -267,6 +276,7 @@ public:
                 this->type_ == rhs.type()
              && this->index_ == rhs.index()
              && this->processor_ == rhs.procIndex()
+             && this->vertexFixed_ == rhs.fixed()
             );
         }
 

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

@@ -35,7 +35,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex()
     index_(-1),
     processor_(Foam::Pstream::myProcNo()),
     alignment_(Foam::triad::unset),
-    targetCellSize_(0.0)
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -47,7 +48,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
     index_(-1),
     processor_(Foam::Pstream::myProcNo()),
     alignment_(Foam::triad::unset),
-    targetCellSize_(0.0)
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -63,7 +65,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
     index_(-1),
     processor_(Foam::Pstream::myProcNo()),
     alignment_(Foam::triad::unset),
-    targetCellSize_(0.0)
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -79,7 +82,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
     index_(-1),
     processor_(Foam::Pstream::myProcNo()),
     alignment_(Foam::triad::unset),
-    targetCellSize_(0.0)
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -96,8 +100,28 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
     type_(type),
     index_(index),
     processor_(processor),
-    alignment_(Foam::tensor::zero),
-    targetCellSize_(0.0)
+    alignment_(Foam::triad::unset),
+    targetCellSize_(0.0),
+    vertexFixed_(false)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
+(
+    const Foam::point& p,
+    int index,
+    vertexType type,
+    int processor
+)
+:
+    Vb(Point(p.x(), p.y(), p.z())),
+    type_(type),
+    index_(index),
+    processor_(processor),
+    alignment_(Foam::triad::unset),
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -109,7 +133,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p, Cell_handle f)
     index_(-1),
     processor_(Foam::Pstream::myProcNo()),
     alignment_(Foam::triad::unset),
-    targetCellSize_(0.0)
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -121,7 +146,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
     index_(-1),
     processor_(Foam::Pstream::myProcNo()),
     alignment_(Foam::triad::unset),
-    targetCellSize_(0.0)
+    targetCellSize_(0.0),
+    vertexFixed_(false)
 {}
 
 
@@ -317,18 +343,18 @@ inline bool CGAL::indexedVertex<Gt, Vb>::surfacePoint() const
 }
 
 
-//template<class Gt, class Vb>
-//inline bool CGAL::indexedVertex<Gt, Vb>::isVertexFixed() const
-//{
-//    return vertexFixed_;
-//}
-//
-//
-//template<class Gt, class Vb>
-//inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed()
-//{
-//    vertexFixed_ = true;
-//}
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::fixed() const
+{
+    return vertexFixed_;
+}
+
+
+template<class Gt, class Vb>
+inline bool& CGAL::indexedVertex<Gt, Vb>::fixed()
+{
+    return vertexFixed_;
+}
 
 
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //

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

@@ -100,8 +100,8 @@ void Foam::conformationSurfaces::hasBoundedVolume
     Info<< "    Sum of all the surface normals (if near zero, surface is"
         << " probably closed):" << nl
         << "    Note: Does not include baffle surfaces in calculation" << nl
-        << "        Sum = " << sum/totalTriangles << nl
-        << "        mag(Sum) = " << mag(sum)/totalTriangles
+        << "        Sum = " << sum/(totalTriangles + SMALL) << nl
+        << "        mag(Sum) = " << mag(sum)/(totalTriangles + SMALL)
         << endl;
 }
 
@@ -641,16 +641,16 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
         {
             const label regionI = regionOffset_[s];
 
-            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
+//            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
 
-            if
-            (
-                !surface.hasVolumeType()
-             && !surface.bounds().contains(samplePts[i])
-            )
-            {
-                continue;
-            }
+//            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
@@ -661,7 +661,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
                 continue;
             }
 
+//            Info<< surface.name() << " = "
+//                << volumeType::names[surfaceVolumeTests[s][i]] << endl;
+
             if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
+            //if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
             {
                 insidePoint[i] = false;
 
@@ -722,7 +726,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellOutside
     // distanceSquared
 
     // Assume that the point is wellInside until demonstrated otherwise.
-    Field<bool> outsidePoint(samplePts.size(), true);
+    Field<bool> outsidePoint(samplePts.size(), false);
 
     //Check if the points are inside the surface by the given distance squared
 
@@ -755,16 +759,16 @@ Foam::Field<bool> Foam::conformationSurfaces::wellOutside
 
         forAll(surfaces_, s)
         {
-            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
+//            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
 
-            if
-            (
-                !surface.hasVolumeType()
-             && !surface.bounds().contains(samplePts[i])
-            )
-            {
-                continue;
-            }
+//            if
+//            (
+//                !surface.hasVolumeType()
+//             //&& !surface.bounds().contains(samplePts[i])
+//            )
+//            {
+//                continue;
+//            }
 
             const label regionI = regionOffset_[s];
 
@@ -777,9 +781,9 @@ Foam::Field<bool> Foam::conformationSurfaces::wellOutside
                 continue;
             }
 
-            if (surfaceVolumeTests[s][i] == volumeType::INSIDE)
+            if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
             {
-                outsidePoint[i] = false;
+                outsidePoint[i] = true;
 
                 break;
             }

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

@@ -58,7 +58,8 @@ initialPointsMethod::initialPointsMethod
                 initialPointsDict.lookup("minimumSurfaceDistanceCoeff")
             )
         )
-    )
+    ),
+    fixInitialPoints_(Switch(initialPointsDict.lookup("fixInitialPoints")))
 {}
 
 

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

@@ -73,6 +73,9 @@ protected:
         //  the local target cell size.  Store square of value.
         scalar minimumSurfaceDistanceCoeffSqr_;
 
+        Switch fixInitialPoints_;
+
+
 private:
 
     // Private Member Functions
@@ -138,6 +141,11 @@ public:
             return detailsDict_;
         }
 
+        Switch fixInitialPoints() const
+        {
+            return fixInitialPoints_;
+        }
+
         //- Return the initial points for the conformalVoronoiMesh
         virtual List<Vb::Point> initialPoints() const = 0;
 };

applications/utilities/mesh/generation/foamyHexMesh/foamyHexMeshDict

@@ -10,6 +10,7 @@ FoamFile
     version         2.0;
     format          ascii;
     class           dictionary;
+    location        "system";
     object          foamyHexMeshDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/surface/surfaceHookUp/surfaceHookUp.C

@@ -201,6 +201,71 @@ void createBoundaryEdgeTrees
 }
 
 
+class findNearestOpSubset
+{
+    const indexedOctree<treeDataEdge>& tree_;
+
+    DynamicList<label>& shapeMask_;
+
+public:
+
+    findNearestOpSubset
+    (
+        const indexedOctree<treeDataEdge>& tree,
+        DynamicList<label>& shapeMask
+    )
+    :
+        tree_(tree),
+        shapeMask_(shapeMask)
+    {}
+
+    void operator()
+    (
+        const labelUList& indices,
+        const point& sample,
+
+        scalar& nearestDistSqr,
+        label& minIndex,
+        point& nearestPoint
+    ) const
+    {
+        const treeDataEdge& shape = tree_.shapes();
+
+        forAll(indices, i)
+        {
+            const label index = indices[i];
+            const label edgeIndex = shape.edgeLabels()[index];
+
+            if
+            (
+                !shapeMask_.empty()
+             && findIndex(shapeMask_, edgeIndex) != -1
+            )
+            {
+                continue;
+            }
+
+            const edge& e = shape.edges()[edgeIndex];
+
+            pointHit nearHit = e.line(shape.points()).nearestDist(sample);
+
+            // Only register hit if closest point is not an edge point
+            if (nearHit.hit())
+            {
+                scalar distSqr = sqr(nearHit.distance());
+
+                if (distSqr < nearestDistSqr)
+                {
+                    nearestDistSqr = distSqr;
+                    minIndex = index;
+                    nearestPoint = nearHit.rawPoint();
+                }
+            }
+        }
+    }
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -226,7 +291,8 @@ int main(int argc, char *argv[])
     const IOdictionary dict(dictIO);
 
     const scalar dist(args.argRead<scalar>(1));
-    const scalar matchTolerance(SMALL);
+    const scalar matchTolerance(max(1e-6*dist, SMALL));
+    const label maxIters = 100;
 
     Info<< "Hooking distance = " << dist << endl;
 
@@ -289,7 +355,7 @@ int main(int argc, char *argv[])
     }
 
     label nChanged = 0;
-    label nIters = 0;
+    label nIters = 1;
 
     do
     {
@@ -339,7 +405,7 @@ int main(int argc, char *argv[])
                         (
                             samplePt,
                             sqr(dist),
-                            treeDataEdge::findNearestOpSubset
+                            findNearestOpSubset
                             (
                                 bEdgeTree,
                                 shapeMask
@@ -439,14 +505,18 @@ int main(int argc, char *argv[])
 
                     nChanged++;
 
+                    label newPointI = -1;
+
                     // Keep the points in the same place and move the edge
-                    newPoints[hitSurfI].append(newPoints[surfI][pointI]);
-
-                    // Move the points to the edges
-                    //newPoints[pointI] = eHit.hitPoint();
-                    //newPoints.append(eHit.hitPoint());
-
-                    visitedFace[hitSurfI][faceI] = true;
+                    if (hitSurfI == surfI)
+                    {
+                        newPointI = pointI;
+                    }
+                    else
+                    {
+                        newPoints[hitSurfI].append(newPoints[surfI][pointI]);
+                        newPointI = newPoints[hitSurfI].size() - 1;
+                    }
 
                     // Split the other face.
                     greenRefine
@@ -454,22 +524,23 @@ int main(int argc, char *argv[])
                         hitSurf,
                         faceI,
                         eIndex,
-                        newPoints[hitSurfI].size() - 1,
+                        newPointI,
                         newFacesFromSplit
                     );
 
+                    visitedFace[hitSurfI][faceI] = true;
+
                     forAll(newFacesFromSplit, newFaceI)
                     {
+                        const labelledTri& fN = newFacesFromSplit[newFaceI];
+
                         if (newFaceI == 0)
                         {
-                            newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
+                            newFaces[hitSurfI][faceI] = fN;
                         }
                         else
                         {
-                            newFaces[hitSurfI].append
-                            (
-                                newFacesFromSplit[newFaceI]
-                            );
+                            newFaces[hitSurfI].append(fN);
                         }
                     }
                 }
@@ -502,7 +573,7 @@ int main(int argc, char *argv[])
             );
         }
 
-    } while (nChanged > 0);
+    } while (nChanged > 0 && nIters <= maxIters);
 
     Info<< endl;
 

src/meshTools/indexedOctree/treeDataEdge.C

@@ -105,17 +105,6 @@ Foam::treeDataEdge::findNearestOp::findNearestOp
 {}
 
 
-Foam::treeDataEdge::findNearestOpSubset::findNearestOpSubset
-(
-    const indexedOctree<treeDataEdge>& tree,
-    DynamicList<label>& shapeMask
-)
-:
-    tree_(tree),
-    shapeMask_(shapeMask)
-{}
-
-
 Foam::treeDataEdge::findIntersectOp::findIntersectOp
 (
     const indexedOctree<treeDataEdge>& tree
@@ -281,105 +270,6 @@ void Foam::treeDataEdge::findNearestOp::operator()
 }
 
 
-void Foam::treeDataEdge::findNearestOpSubset::operator()
-(
-    const labelUList& indices,
-    const point& sample,
-
-    scalar& nearestDistSqr,
-    label& minIndex,
-    point& nearestPoint
-) const
-{
-    const treeDataEdge& shape = tree_.shapes();
-
-    forAll(indices, i)
-    {
-        const label index = indices[i];
-        const label edgeIndex = shape.edgeLabels()[index];
-
-        if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
-        {
-            continue;
-        }
-
-        const edge& e = shape.edges()[edgeIndex];
-
-        pointHit nearHit = e.line(shape.points()).nearestDist(sample);
-
-        scalar distSqr = sqr(nearHit.distance());
-
-        if (distSqr < nearestDistSqr)
-        {
-            nearestDistSqr = distSqr;
-            minIndex = index;
-            nearestPoint = nearHit.rawPoint();
-        }
-    }
-}
-
-
-void Foam::treeDataEdge::findNearestOpSubset::operator()
-(
-    const labelUList& indices,
-    const linePointRef& ln,
-
-    treeBoundBox& tightest,
-    label& minIndex,
-    point& linePoint,
-    point& nearestPoint
-) const
-{
-    const treeDataEdge& shape = tree_.shapes();
-
-    // Best so far
-    scalar nearestDistSqr = magSqr(linePoint - nearestPoint);
-
-    forAll(indices, i)
-    {
-        const label index = indices[i];
-        const label edgeIndex = shape.edgeLabels()[index];
-
-        if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
-        {
-            continue;
-        }
-
-        const edge& e = shape.edges()[edgeIndex];
-
-        // Note: could do bb test ? Worthwhile?
-
-        // Nearest point on line
-        point ePoint, lnPt;
-        scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt);
-        scalar distSqr = sqr(dist);
-
-        if (distSqr < nearestDistSqr)
-        {
-            nearestDistSqr = distSqr;
-            minIndex = index;
-            linePoint = lnPt;
-            nearestPoint = ePoint;
-
-            {
-                point& minPt = tightest.min();
-                minPt = min(ln.start(), ln.end());
-                minPt.x() -= dist;
-                minPt.y() -= dist;
-                minPt.z() -= dist;
-            }
-            {
-                point& maxPt = tightest.max();
-                maxPt = max(ln.start(), ln.end());
-                maxPt.x() += dist;
-                maxPt.y() += dist;
-                maxPt.z() += dist;
-            }
-        }
-    }
-}
-
-
 bool Foam::treeDataEdge::findIntersectOp::operator()
 (
     const label index,

src/meshTools/indexedOctree/treeDataEdge.H

@@ -118,43 +118,6 @@ public:
         ) const;
     };
 
-    class findNearestOpSubset
-    {
-        const indexedOctree<treeDataEdge>& tree_;
-
-        DynamicList<label>& shapeMask_;
-
-    public:
-
-        findNearestOpSubset
-        (
-            const indexedOctree<treeDataEdge>& tree,
-            DynamicList<label>& shapeMask
-        );
-
-        void operator()
-        (
-            const labelUList& indices,
-            const point& sample,
-
-            scalar& nearestDistSqr,
-            label& minIndex,
-            point& nearestPoint
-        ) const;
-
-        //- Calculates nearest (to line) point in shape.
-        //  Returns point and distance (squared)
-        void operator()
-        (
-            const labelUList& indices,
-            const linePointRef& ln,
-
-            treeBoundBox& tightest,
-            label& minIndex,
-            point& linePoint,
-            point& nearestPoint
-        ) const;
-    };
 
     class findIntersectOp
     {

src/meshTools/searchableSurface/searchableSurface.C

@@ -39,7 +39,6 @@ defineRunTimeSelectionTable(searchableSurface, dict);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-
 Foam::autoPtr<Foam::searchableSurface> Foam::searchableSurface::New
 (
     const word& searchableSurfaceType,

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -738,9 +738,19 @@ void Foam::triSurfaceMesh::getVolumeType
     {
         const point& pt = points[pointI];
 
-        // - use cached volume type per each tree node
-        // - cheat conversion since same values
-        volType[pointI] = tree().getVolumeType(pt);
+        if (!tree().bb().contains(pt))
+        {
+            // Have to calculate directly as outside the octree
+            volType[pointI] = tree().shapes().getVolumeType(tree(), pt);
+        }
+        else
+        {
+            // - use cached volume type per each tree node
+            volType[pointI] = tree().getVolumeType(pt);
+        }
+
+//        Info<< "octree : " << pt << " = "
+//            << volumeType::names[volType[pointI]] << endl;
     }
 
     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;

tutorials/mesh/foamyHex2DMesh/OpenCFD/Allrun

@@ -7,7 +7,7 @@ cd ${0%/*} || exit 1    # run from this directory
 cp system/controlDict.mesher system/controlDict
 
 runApplication surfaceFeatureExtract
-runApplication cv2DMesh -overwrite
+runApplication foamyHex2DMesh -overwrite
 runApplication extrude2DMesh -overwrite polyMesh2D
 runApplication checkMesh -allGeometry -allTopology -constant -noZero
 

tutorials/mesh/foamyHex2DMesh/OpenCFD/system/foamyHex2DMeshDict

@@ -10,7 +10,8 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      cv2DMeshDict;
+    location    "system";
+    object      foamyHex2DMeshDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -98,7 +99,7 @@ motionControl
         {
             type                  searchableSurfaceControl;
             priority              1;
-            mode                  bothSides;
+            mode                  inside;
             cellSizeFunction      uniform;
 
             linearDistanceCoeffs

tutorials/mesh/foamyHex2DMesh/jaggedBoundary/Allrun

@@ -5,7 +5,7 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
 runApplication surfaceFeatureExtract
-runApplication cv2DMesh -overwrite
+runApplication foamyHex2DMesh -overwrite
 runApplication extrude2DMesh -overwrite MeshedSurface
 runApplication checkMesh -allGeometry -allTopology -constant -noZero
 

tutorials/mesh/foamyHex2DMesh/jaggedBoundary/system/foamyHex2DMeshDict

@@ -10,7 +10,8 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      cv2DMeshDict;
+    location    "system";
+    object      foamyHex2DMeshDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/mesh/foamyHex2DMesh/square/Allrun

@@ -5,7 +5,7 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
 runApplication surfaceFeatureExtract
-runApplication cv2DMesh -overwrite
+runApplication foamyHex2DMesh -overwrite
 runApplication extrude2DMesh -overwrite polyMesh2D
 runApplication checkMesh -allGeometry -allTopology -constant -noZero
 

tutorials/mesh/foamyHex2DMesh/square/system/foamyHex2DMeshDict

@@ -10,7 +10,8 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      DMeshDict;
+    location    "system";
+    object      foamyHex2DMeshDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //