ENH: snappyHexMesh: added debug switch

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -418,9 +418,9 @@ addLayersControls
         // default is 0.
         //nSmoothDisplacement 90;
 
-        // Optional: do not extrude any point where all surrounding faces
+        // (wip)Optional: do not extrude a point if none of the surrounding points is
-        // have at least one point not extruded. Default is true.
+        // not extruded. Default is false.
-        //detectExtrusionIsland false;
+        //detectExtrusionIsland true;
 
 
     // Mesh shrinking
@@ -488,7 +488,9 @@ meshQualityControls
 //(
 //    mesh            // write intermediate meshes
 //    intersections   // write current mesh intersections as .obj files
-//    featureSeeds    // write information about explicit feature edge refinement
+//    featureSeeds    // write information about explicit feature edge
+//                    // refinement
+//    attraction      // write attraction as .obj files
 //    layerInfo       // write information about layers
 //);
 //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -41,7 +41,7 @@ Description
 #include "mapPolyMesh.H"
 #include "addPatchCellLayer.H"
 #include "mapDistributePolyMesh.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "layerParameters.H"
 #include "combineFaces.H"
 #include "IOmanip.H"
@@ -80,37 +80,25 @@ void Foam::autoLayerDriver::dumpDisplacement
     const List<extrudeMode>& extrudeStatus
 )
 {
-    OFstream dispStr(prefix + "_disp.obj");
+    OBJstream dispStr(prefix + "_disp.obj");
     Info<< "Writing all displacements to " << dispStr.name() << endl;
 
-    label vertI = 0;
-
     forAll(patchDisp, patchPointI)
     {
         const point& pt = pp.localPoints()[patchPointI];
-
+        dispStr.write(linePointRef(pt, pt + patchDisp[patchPointI]));
-        meshTools::writeOBJ(dispStr, pt); vertI++;
-        meshTools::writeOBJ(dispStr, pt + patchDisp[patchPointI]); vertI++;
-
-        dispStr << "l " << vertI-1 << ' ' << vertI << nl;
     }
 
 
-    OFstream illStr(prefix + "_illegal.obj");
+    OBJstream illStr(prefix + "_illegal.obj");
     Info<< "Writing invalid displacements to " << illStr.name() << endl;
 
-    vertI = 0;
-
     forAll(patchDisp, patchPointI)
     {
         if (extrudeStatus[patchPointI] != EXTRUDE)
         {
             const point& pt = pp.localPoints()[patchPointI];
-
+            illStr.write(linePointRef(pt, pt + patchDisp[patchPointI]));
-            meshTools::writeOBJ(illStr, pt); vertI++;
-            meshTools::writeOBJ(illStr, pt + patchDisp[patchPointI]); vertI++;
-
-            illStr << "l " << vertI-1 << ' ' << vertI << nl;
         }
     }
 }
@@ -441,13 +429,12 @@ void Foam::autoLayerDriver::handleFeatureAngle
             point::max          // null value
         );
 
-        label vertI = 0;
+        autoPtr<OBJstream> str;
-        autoPtr<OFstream> str;
         if (debug&meshRefinement::MESH)
         {
             str.reset
             (
-                new OFstream
+                new OBJstream
                 (
                     mesh.time().path()
                   / "featureEdges_"
@@ -497,11 +484,9 @@ void Foam::autoLayerDriver::handleFeatureAngle
 
                     if (str.valid())
                     {
-                        meshTools::writeOBJ(str(), pp.localPoints()[e[0]]);
+                        const point& p0 = pp.localPoints()[e[0]];
-                        vertI++;
+                        const point& p1 = pp.localPoints()[e[1]];
-                        meshTools::writeOBJ(str(), pp.localPoints()[e[1]]);
+                        str().write(linePointRef(p0, p1));
-                        vertI++;
-                        str()<< "l " << vertI-1 << ' ' << vertI << nl;
                     }
                 }
             }
@@ -1542,7 +1527,7 @@ void Foam::autoLayerDriver::getPatchDisplacement
 
             if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
             {
-                if (debug&meshRefinement::OBJINTERSECTIONS)
+                if (debug&meshRefinement::ATTRACTION)
                 {
                     Pout<< "No valid normal for point " << meshPointI
                         << ' ' << pp.points()[meshPointI]
@@ -1582,7 +1567,7 @@ void Foam::autoLayerDriver::getPatchDisplacement
 
             if (nPoints > 0)
             {
-                if (debug&meshRefinement::OBJINTERSECTIONS)
+                if (debug&meshRefinement::ATTRACTION)
                 {
                     Pout<< "Displacement at illegal point "
                         << localPoints[patchPointI]

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -94,7 +94,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
             (
                 meshRefiner_.refineCandidates
                 (
-                    refineParams.keepPoints()[0],    // For now only use one.
+                    refineParams.keepPoints(),
                     refineParams.curvature(),
                     refineParams.planarAngle(),
 
@@ -207,7 +207,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
         (
             meshRefiner_.refineCandidates
             (
-                refineParams.keepPoints()[0],
+                refineParams.keepPoints(),
                 refineParams.curvature(),
                 refineParams.planarAngle(),
 
@@ -341,7 +341,7 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
         (
             meshRefiner_.refineCandidates
             (
-                refineParams.keepPoints()[0],
+                refineParams.keepPoints(),
                 refineParams.curvature(),
                 refineParams.planarAngle(),
 
@@ -753,7 +753,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
         (
             meshRefiner_.refineCandidates
             (
-                refineParams.keepPoints()[0],
+                refineParams.keepPoints(),
                 refineParams.curvature(),
                 refineParams.planarAngle(),
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -1279,7 +1279,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
 
     autoPtr<OBJstream> gapStr;
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         gapStr.reset
         (
@@ -3089,7 +3089,7 @@ void Foam::autoSnapDriver::doSnap
             meshMover.setDisplacement(disp);
 
 
-            if (debug&meshRefinement::OBJINTERSECTIONS)
+            if (debug&meshRefinement::ATTRACTION)
             {
                 dumpMove
                 (
@@ -3120,13 +3120,6 @@ void Foam::autoSnapDriver::doSnap
                     << "    surfaces but the"
                     << " resulting mesh will not satisfy your quality"
                     << " constraints" << nl << endl;
-                //Info<< "Did not succesfully snap mesh. Giving up."
-                //    << nl << endl;
-                //
-                //// Use current mesh as base mesh
-                //meshMover.correct();
-                //
-                //break;
             }
 
             if (debug&meshRefinement::MESH)

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C

@@ -422,7 +422,7 @@ void Foam::autoSnapDriver::calcNearestFace
           ^ faceSurfaceNormal[faceI];
     }
 
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         dumpMove
         (
@@ -1037,7 +1037,7 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 {
     autoPtr<OBJstream> feStr;
     autoPtr<OBJstream> fpStr;
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         feStr.reset
         (
@@ -1331,7 +1331,7 @@ void Foam::autoSnapDriver::releasePointsNextToMultiPatch
 ) const
 {
     autoPtr<OBJstream> multiPatchStr;
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         multiPatchStr.reset
         (
@@ -1809,7 +1809,7 @@ void Foam::autoSnapDriver::determineFeatures
     autoPtr<OBJstream> missedEdgeStr;
     autoPtr<OBJstream> featurePointStr;
 
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         featureEdgeStr.reset
         (
@@ -2263,7 +2263,7 @@ void Foam::autoSnapDriver::determineBaffleFeatures
 
 
     autoPtr<OBJstream> baffleEdgeStr;
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         baffleEdgeStr.reset
         (
@@ -2865,7 +2865,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
 
     // Dump
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         OBJstream featureEdgeStr
         (
@@ -2962,7 +2962,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
         );
     }
 
-    if (debug&meshRefinement::OBJINTERSECTIONS)
+    if (debug&meshRefinement::ATTRACTION)
     {
         dumpMove
         (
@@ -3380,7 +3380,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         // Re-add normal component
         tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
 
-        if (debug&meshRefinement::OBJINTERSECTIONS)
+        if (debug&meshRefinement::ATTRACTION)
         {
             dumpMove
             (

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -67,13 +67,14 @@ namespace Foam
     const char* Foam::NamedEnum
     <
         Foam::meshRefinement::IOdebugType,
-        4
+        5
     >::names[] =
     {
         "mesh",
         //"scalarLevels",
         "intersections",
         "featureSeeds",
+        "attraction",
         "layerInfo"
     };
 
@@ -102,7 +103,7 @@ namespace Foam
 
 }
 
-const Foam::NamedEnum<Foam::meshRefinement::IOdebugType, 4>
+const Foam::NamedEnum<Foam::meshRefinement::IOdebugType, 5>
 Foam::meshRefinement::IOdebugTypeNames;
 
 const Foam::NamedEnum<Foam::meshRefinement::IOoutputType, 1>

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -92,15 +92,17 @@ public:
             //IOSCALARLEVELS,
             IOOBJINTERSECTIONS,
             IOFEATURESEEDS,
+            IOATTRACTION,
             IOLAYERINFO
         };
-        static const NamedEnum<IOdebugType, 4> IOdebugTypeNames;
+        static const NamedEnum<IOdebugType, 5> IOdebugTypeNames;
         enum debugType
         {
             MESH = 1<<IOMESH,
             //SCALARLEVELS = 1<<IOSCALARLEVELS,
             OBJINTERSECTIONS = 1<<IOOBJINTERSECTIONS,
             FEATURESEEDS = 1<<IOFEATURESEEDS,
+            ATTRACTION = 1<< IOATTRACTION,
             LAYERINFO = 1<<IOLAYERINFO
         };
 
@@ -256,7 +258,7 @@ private:
             //  (or -1 if not passed through). Uses tracking.
             void markFeatureCellLevel
             (
-                const point& keepPoint,
+                const pointField& keepPoints,
                 labelList& maxFeatureLevel
             ) const;
 
@@ -264,7 +266,7 @@ private:
             //  features.
             label markFeatureRefinement
             (
-                const point& keepPoint,
+                const pointField& keepPoints,
                 const label nAllowRefine,
 
                 labelList& refineCell,
@@ -779,7 +781,7 @@ public:
             //- Calculate list of cells to refine.
             labelList refineCandidates
             (
-                const point& keepPoint,
+                const pointField& keepPoints,
                 const scalar curvature,
                 const scalar planarAngle,
 
@@ -914,7 +916,7 @@ public:
             //- Select coupled faces that are not collocated
             void selectSeparatedCoupledFaces(boolList&) const;
 
-            //- Find region cell is in. Uses optional perturbation to re-test.
+            //- Find region point is in. Uses optional perturbation to re-test.
             static label findRegion
             (
                 const polyMesh&,
@@ -1018,6 +1020,8 @@ public:
                     const labelHashSet& set
                 ) const;
 
+                //- Merge edges, maintain mesh quality. Return global number
+                //  of edges merged
                 label mergeEdgesUndo
                 (
                     const scalar minCos,

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -313,7 +313,7 @@ bool Foam::meshRefinement::markForRefine
 
 void Foam::meshRefinement::markFeatureCellLevel
 (
-    const point& keepPoint,
+    const pointField& keepPoints,
     labelList& maxFeatureLevel
 ) const
 {
@@ -345,6 +345,10 @@ void Foam::meshRefinement::markFeatureCellLevel
     // what to seed. Do this on only the processor that
     // holds the keepPoint.
 
+    forAll(keepPoints, i)
+    {
+        const point& keepPoint = keepPoints[i];
+
         label cellI = -1;
         label tetFaceI = -1;
         label tetPtI = -1;
@@ -449,6 +453,7 @@ void Foam::meshRefinement::markFeatureCellLevel
                 }
             }
         }
+    }
 
 
     // Largest refinement level of any feature passed through
@@ -628,7 +633,7 @@ void Foam::meshRefinement::markFeatureCellLevel
 // Calculates list of cells to refine based on intersection with feature edge.
 Foam::label Foam::meshRefinement::markFeatureRefinement
 (
-    const point& keepPoint,
+    const pointField& keepPoints,
     const label nAllowRefine,
 
     labelList& refineCell,
@@ -637,7 +642,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 {
     // Largest refinement level of any feature passed through
     labelList maxFeatureLevel;
-    markFeatureCellLevel(keepPoint, maxFeatureLevel);
+    markFeatureCellLevel(keepPoints, maxFeatureLevel);
 
     // See which cells to refine. maxFeatureLevel will hold highest level
     // of any feature edge that passed through.
@@ -2009,7 +2014,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
 // hitting overall limit maxGlobalCells.
 Foam::labelList Foam::meshRefinement::refineCandidates
 (
-    const point& keepPoint,
+    const pointField& keepPoints,
     const scalar curvature,
     const scalar planarAngle,
 
@@ -2077,7 +2082,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         {
             label nFeatures = markFeatureRefinement
             (
-                keepPoint,
+                keepPoints,
                 nAllowRefine,
 
                 refineCell,