ENH: foamyHexMesh: add face/cell zoning

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

@@ -21,6 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
     -IvectorTools
 
 EXE_LIBS = \

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

@@ -279,11 +279,9 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
     const Triangulation& t
 )
 {
-    tmp<pointField> tpts(new pointField(t.number_of_vertices(), point::max));
+    tmp<pointField> tpts(new pointField(t.vertexCount(), point::max));
     pointField& pts = tpts();
 
-    label nVert = 0;
-
     for
     (
         typename Triangulation::Finite_vertices_iterator vit =
@@ -292,9 +290,9 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
         ++vit
     )
     {
-        if (vit->internalOrBoundaryPoint())
+        if (vit->internalOrBoundaryPoint() && !vit->referred())
         {
-            pts[nVert++] = topoint(vit->point());
+            pts[vit->index()] = topoint(vit->point());
         }
     }
 

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

@@ -22,6 +22,7 @@ EXE_INC = \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
     -IPrintTable \
     -I../vectorTools
 
@@ -32,4 +33,5 @@ LIB_LIBS = \
     -ltriSurface \
     -ldynamicMesh \
     -lsurfMesh \
-    -lsampling
+    -lsampling \
+    -lautoMesh

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

@@ -600,6 +600,29 @@ private:
             PackedBoolList& boundaryFacesToRemove
         );
 
+        //- From meshRefinementBaffles.C. Use insidePoint for a surface to
+        //  determine the cell zone.
+        void findCellZoneInsideWalk
+        (
+            const polyMesh& mesh,
+            const labelList& locationSurfaces,
+            const labelList& faceToSurface,
+            labelList& cellToSurface
+        ) const;
+
+        //- Calculate the cell zones from cellCentres using all closed surfaces
+        labelList calcCellZones(const pointField& cellCentres) const;
+
+        //- Calculate the face zones
+        void calcFaceZones
+        (
+            const polyMesh& mesh,
+            const pointField& cellCentres,
+            const labelList& cellToSurface,
+            labelList& faceToSurface,
+            boolList& flipMap
+        ) const;
+
         //- Tet mesh calculation
         void calcTetMesh
         (
@@ -712,9 +735,6 @@ private:
             bool includeEmptyPatches = false
         ) const;
 
-        //- Create the cell centres to use for the mesh
-        void createCellCentres(pointField& cellCentres) const;
-
         //- Sort the faces, owner and neighbour lists into
         //  upper-triangular order.  For internal faces only, use
         //  before adding patch faces

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

@@ -2852,32 +2852,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 }
 
 
-void Foam::conformalVoronoiMesh::createCellCentres
-(
-    pointField& cellCentres
-) const
-{
-    cellCentres.setSize(number_of_vertices(), point::max);
-
-    label vertI = 0;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->internalOrBoundaryPoint())
-        {
-            cellCentres[vertI++] = topoint(vit->point());
-        }
-    }
-
-    cellCentres.setSize(vertI);
-}
-
-
 void Foam::conformalVoronoiMesh::sortFaces
 (
     faceList& faces,
@@ -3095,7 +3069,7 @@ Foam::labelList Foam::conformalVoronoiMesh::removeUnusedCells
 {
     Info<< nl << "Removing unused cells" << endl;
 
-    PackedBoolList cellUsed(number_of_vertices(), false);
+    PackedBoolList cellUsed(vertexCount(), false);
 
     // Scan all faces to find all of the cells that are used
 

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

@@ -35,6 +35,11 @@ License
 #include "pointMesh.H"
 #include "indexedVertexOps.H"
 #include "DelaunayMeshTools.H"
+#include "surfaceZonesInfo.H"
+#include "polyModifyCell.H"
+#include "polyModifyFace.H"
+#include "syncTools.H"
+#include "regionSplit.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -402,6 +407,376 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 }
 
 
+void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
+(
+    const polyMesh& mesh,
+    const labelList& locationSurfaces,  // indices of surfaces with inside point
+    const labelList& faceToSurface, // per face index of named surface
+    labelList& cellToSurface
+) const
+{
+    // Analyse regions. Reuse regionsplit
+    boolList blockedFace(mesh.nFaces());
+    //selectSeparatedCoupledFaces(blockedFace);
+
+    forAll(faceToSurface, faceI)
+    {
+        if (faceToSurface[faceI] == -1)
+        {
+            blockedFace[faceI] = false;
+        }
+        else
+        {
+            blockedFace[faceI] = true;
+        }
+    }
+    // No need to sync since namedSurfaceIndex already is synced
+
+    // Set region per cell based on walking
+    regionSplit cellRegion(mesh, blockedFace);
+    blockedFace.clear();
+
+
+    // Force calculation of face decomposition (used in findCell)
+    (void)mesh.tetBasePtIs();
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    // For all locationSurface find the cell
+    forAll(locationSurfaces, i)
+    {
+        label surfI = locationSurfaces[i];
+
+        const Foam::point& insidePoint = surfZones[surfI].zoneInsidePoint();
+
+        const word& surfName = geometryToConformTo().geometry().names()[surfI];
+
+        Info<< "    For surface " << surfName
+            << " finding inside point " << insidePoint
+            << endl;
+
+        // Find the region containing the insidePoint
+        label keepRegionI = -1;
+
+        label cellI = mesh.findCell(insidePoint);
+
+        if (cellI != -1)
+        {
+            keepRegionI = cellRegion[cellI];
+        }
+        reduce(keepRegionI, maxOp<label>());
+
+        Info<< "    For surface " << surfName
+            << " found point " << insidePoint << " in cell " << cellI
+            << " in global region " << keepRegionI
+            << " out of " << cellRegion.nRegions() << " regions." << endl;
+
+        if (keepRegionI == -1)
+        {
+            FatalErrorIn
+            (
+                "conformalVoronoiMesh::findCellZoneInsideWalk"
+                "(const polyMesh&, const labelList&"
+                ", const labelList&, labelList&)"
+            )   << "Point " << insidePoint
+                << " is not inside the mesh." << nl
+                << "Bounding box of the mesh:" << mesh.bounds()
+                << exit(FatalError);
+        }
+
+        // Set all cells with this region
+        forAll(cellRegion, cellI)
+        {
+            if (cellRegion[cellI] == keepRegionI)
+            {
+                if (cellToSurface[cellI] == -2)
+                {
+                    cellToSurface[cellI] = surfI;
+                }
+                else if (cellToSurface[cellI] != surfI)
+                {
+                    WarningIn
+                    (
+                        "conformalVoronoiMesh::findCellZoneInsideWalk"
+                        "(const labelList&, const labelList&"
+                        ", const labelList&, const labelList&)"
+                    )   << "Cell " << cellI
+                        << " at " << mesh.cellCentres()[cellI]
+                        << " is inside surface " << surfName
+                        << " but already marked as being in zone "
+                        << cellToSurface[cellI] << endl
+                        << "This can happen if your surfaces are not"
+                        << " (sufficiently) closed."
+                        << endl;
+                }
+            }
+        }
+    }
+}
+
+
+Foam::labelList Foam::conformalVoronoiMesh::calcCellZones
+(
+    const pointField& cellCentres
+) const
+{
+    labelList cellToSurface(cellCentres.size(), -1);
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    // Get list of closed surfaces
+    labelList closedNamedSurfaces
+    (
+        surfaceZonesInfo::getAllClosedNamedSurfaces
+        (
+            surfZones,
+            geometryToConformTo().geometry(),
+            geometryToConformTo().surfaces()
+        )
+    );
+
+    forAll(closedNamedSurfaces, i)
+    {
+        label surfI = closedNamedSurfaces[i];
+
+        const searchableSurface& surface =
+            allGeometry()[geometryToConformTo().surfaces()[surfI]];
+
+        const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
+            surfZones[surfI].zoneInside();
+
+        if
+        (
+            selectionMethod != surfaceZonesInfo::INSIDE
+         && selectionMethod != surfaceZonesInfo::OUTSIDE
+         && selectionMethod != surfaceZonesInfo::INSIDEPOINT
+        )
+        {
+            FatalErrorIn("conformalVoronoiMesh::calcCellZones(..)")
+                << "Trying to use surface "
+                << surface.name()
+                << " which has non-geometric inside selection method "
+                << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
+                << exit(FatalError);
+        }
+
+        if (surface.hasVolumeType())
+        {
+            List<volumeType> volType;
+            surface.getVolumeType(cellCentres, volType);
+
+            bool selectInside = true;
+            if (selectionMethod == surfaceZonesInfo::INSIDEPOINT)
+            {
+                List<volumeType> volTypeInsidePoint;
+                surface.getVolumeType
+                (
+                    pointField(1, surfZones[surfI].zoneInsidePoint()),
+                    volTypeInsidePoint
+                );
+
+                if (volTypeInsidePoint[0] == volumeType::OUTSIDE)
+                {
+                    selectInside = false;
+                }
+            }
+            else if (selectionMethod == surfaceZonesInfo::OUTSIDE)
+            {
+                selectInside = false;
+            }
+
+            forAll(volType, pointI)
+            {
+                if (cellToSurface[pointI] == -1)
+                {
+                    if
+                    (
+                        (
+                            volType[pointI] == volumeType::INSIDE
+                         && selectInside
+                        )
+                     || (
+                            volType[pointI] == volumeType::OUTSIDE
+                         && !selectInside
+                        )
+                    )
+                    {
+                        cellToSurface[pointI] = surfI;
+                    }
+                }
+            }
+        }
+    }
+
+    return cellToSurface;
+}
+
+
+void Foam::conformalVoronoiMesh::calcFaceZones
+(
+    const polyMesh& mesh,
+    const pointField& cellCentres,
+    const labelList& cellToSurface,
+    labelList& faceToSurface,
+    boolList& flipMap
+) const
+{
+    faceToSurface.setSize(mesh.nFaces(), -1);
+    flipMap.setSize(mesh.nFaces(), false);
+
+    const faceList& faces = mesh.faces();
+    const labelList& faceOwner = mesh.faceOwner();
+    const labelList& faceNeighbour = mesh.faceNeighbour();
+
+    forAll(faces, faceI)
+    {
+        const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
+
+        if (mesh.isInternalFace(faceI))
+        {
+            const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
+
+            flipMap[faceI] =
+                (
+                    ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
+                  ? false
+                  : true
+                );
+
+            if
+            (
+                (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
+             && ownerSurfaceI != neiSurfaceI
+            )
+            {
+                if (ownerSurfaceI > neiSurfaceI)
+                {
+                    faceToSurface[faceI] = ownerSurfaceI;
+                }
+                else
+                {
+                    faceToSurface[faceI] = neiSurfaceI;
+                }
+            }
+        }
+        else
+        {
+            if (ownerSurfaceI >= 0)
+            {
+                faceToSurface[faceI] = ownerSurfaceI;
+            }
+        }
+    }
+
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    labelList insidePointNamedSurfaces
+    (
+        surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+    );
+
+    // Use intersection of cellCentre connections
+    forAll(faces, faceI)
+    {
+        if
+        (
+            mesh.isInternalFace(faceI)
+         && faceToSurface[faceI] < 0
+        )
+        {
+            const label own = faceOwner[faceI];
+            const label nei = faceNeighbour[faceI];
+
+            pointIndexHit surfHit;
+            label hitSurface;
+
+            geometryToConformTo().findSurfaceAnyIntersection
+            (
+                cellCentres[own],
+                cellCentres[nei],
+                surfHit,
+                hitSurface
+            );
+
+            if (surfHit.hit())
+            {
+                if (findIndex(insidePointNamedSurfaces, hitSurface) != -1)
+                {
+                    faceToSurface[faceI] = hitSurface;
+
+                    vectorField norm;
+                    geometryToConformTo().getNormal
+                    (
+                        hitSurface,
+                        List<pointIndexHit>(1, surfHit),
+                        norm
+                    );
+
+                    const vector fN = faces[faceI].normal(mesh.points());
+
+                    if ((norm[0] & fN/(mag(fN) + SMALL)) < 0)
+                    {
+                        flipMap[faceI] = true;
+                    }
+                    else
+                    {
+                        flipMap[faceI] = false;
+                    }
+                }
+            }
+        }
+    }
+
+
+    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                label ownSurface = cellToSurface[faceOwner[faceI]];
+                neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                label ownSurface = cellToSurface[faceOwner[faceI]];
+                label neiSurface = neiCellSurface[faceI-mesh.nInternalFaces()];
+
+                if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
+                {
+                    // Give face the max cell zone
+                    faceToSurface[faceI] =  max(ownSurface, neiSurface);
+                }
+            }
+        }
+    }
+
+    // Sync
+    syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
+}
+
+
 Foam::autoPtr<Foam::fvMesh> Foam::conformalVoronoiMesh::createDummyMesh
 (
     const IOobject& io,
@@ -883,7 +1258,11 @@ void Foam::conformalVoronoiMesh::writeMesh
             // Check that the patch is not empty on every processor
             reduce(totalPatchSize, sumOp<label>());
 
-            if (totalPatchSize > 0)
+            if
+            (
+                totalPatchSize > 0
+//             && !geometryToConformTo().surfZones().set(p)
+            )
             {
                 patches[nValidPatches] = polyPatch::New
                 (
@@ -898,6 +1277,145 @@ void Foam::conformalVoronoiMesh::writeMesh
         }
     }
 
+    patches.setSize(nValidPatches);
+
+    mesh.addFvPatches(patches);
+
+
+    // Add zones to mesh
+    {
+        Info<< "    Adding zones to mesh" << endl;
+
+        const PtrList<surfaceZonesInfo>& surfZones =
+            geometryToConformTo().surfZones();
+
+        labelList cellToSurface(calcCellZones(cellCentres));
+
+        labelList faceToSurface;
+        boolList flipMap;
+
+        calcFaceZones
+        (
+            mesh,
+            cellCentres,
+            cellToSurface,
+            faceToSurface,
+            flipMap