ENH: autoLayerDriver: do not disable extrusions around edges

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

@@ -439,6 +439,86 @@ void Foam::autoLayerDriver::handleFeatureAngle
 }
 
 
+//Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::undistortedEdgeLength
+//(
+//    const indirectPrimitivePatch& pp,
+//    const bool relativeSizes,
+//    const bool faceSize
+//)
+//{
+//    const fvMesh& mesh = meshRefiner_.mesh();
+//
+//    tmp<scalarField> tfld(new scalarField());
+//    scalarField& fld = tfld();
+//
+//    if (faceSize)
+//    {
+//        fld.setSize(pp.size());
+//    }
+//    else
+//    {
+//        fld.setSize(pp.nPoints());
+//    }
+//
+//
+//    if (relativeSizes)
+//    {
+//        const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
+//        const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
+//
+//        if (faceSize)
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.addressing()[i];
+//                label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
+//                fld[i] = edge0Len/(1<<ownLevel);
+//            }
+//        }
+//        else
+//        {
+//            // Determine per point the max cell level of connected cells
+//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+//            labelList maxPointLevel(pp.nPoints(), labelMin);
+//
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.addressing()[i];
+//                label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
+//
+//                const face& f = pp.localFaces()[i];
+//                forAll(f, fp)
+//                {
+//                    maxPointLevel[f[fp]] =
+//                        max(maxPointLevel[f[fp]], ownLevel);
+//                }
+//            }
+//
+//            syncTools::syncPointList
+//            (
+//                mesh,
+//                pp.meshPoints(),
+//                maxPointLevel,
+//                maxEqOp<label>(),
+//                labelMin            // null value
+//            );
+//
+//
+//            forAll(maxPointLevel, pointI)
+//            {
+//                // Find undistorted edge size for this level.
+//                fld[i] = edge0Len/(1<<maxPointLevel[pointI]);
+//            }
+//        }
+//    }
+//    else
+//    {
+//        // Use actual cell size
+//    }
+//}
+
+
 // No extrusion on cells with warped faces. Calculates the thickness of the
 // layer and compares it to the space the warped face takes up. Disables
 // extrusion if layer thickness is more than faceRatio of the thickness of
@@ -2383,22 +2463,27 @@ void Foam::autoLayerDriver::addLayers
 
         // Disable extrusion on warped faces
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // It is hard to calculate some length scale if not in relative
+        // mode so disable this check.
+        if (layerParams.relativeSizes())
+        {
+            // Undistorted edge length
+            const scalar edge0Len =
+                meshRefiner_.meshCutter().level0EdgeLength();
+            const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
 
-        // Undistorted edge length
+            handleWarpedFaces
-        const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
+            (
-        const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
+                pp,
+                layerParams.maxFaceThicknessRatio(),
+                edge0Len,
+                cellLevel,
 
-        handleWarpedFaces
+                patchDisp,
-        (
+                patchNLayers,
-            pp,
+                extrudeStatus
-            layerParams.maxFaceThicknessRatio(),
+            );
-            edge0Len,
+        }
-            cellLevel,
-
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
 
         //// Disable extrusion on cells with multiple patch faces
         //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2585,6 +2670,21 @@ void Foam::autoLayerDriver::addLayers
         dimensionedScalar("medialRatio", dimless, 0.0)
     );
 
+    pointVectorField medialVec
+    (
+        IOobject
+        (
+            "medialVec",
+            meshRefiner_.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        meshMover().pMesh(),
+        dimensionedVector("medialVec", dimLength, vector::zero)
+    );
+
     // Setup information for medial axis smoothing. Calculates medial axis
     // and a smoothed displacement direction.
     // - pointMedialDist : distance to medial axis
@@ -2600,7 +2700,8 @@ void Foam::autoLayerDriver::addLayers
 
         dispVec,
         medialRatio,
-        pointMedialDist
+        pointMedialDist,
+        medialVec
     );
 
 
@@ -2692,6 +2793,7 @@ void Foam::autoLayerDriver::addLayers
                 dispVec,
                 medialRatio,
                 pointMedialDist,
+                medialVec,
 
                 extrudeStatus,
                 patchDisp,
@@ -2830,7 +2932,7 @@ void Foam::autoLayerDriver::addLayers
                 mesh.name(),
                 static_cast<polyMesh&>(mesh).instance(),
                 mesh.time(),  // register with runTime
-                static_cast<polyMesh&>(mesh).readOpt(),
+                IOobject::NO_READ,
                 static_cast<polyMesh&>(mesh).writeOpt()
             ),              // io params from original mesh but new name
             mesh,           // original mesh

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -460,7 +460,8 @@ class autoLayerDriver
 
                     pointVectorField& dispVec,
                     pointScalarField& medialRatio,
-                    pointScalarField& medialDist
+                    pointScalarField& medialDist,
+                    pointVectorField& medialVec
                 ) const;
 
                 //- Main routine to shrink mesh
@@ -481,6 +482,7 @@ class autoLayerDriver
                     const pointVectorField& dispVec,
                     const pointScalarField& medialRatio,
                     const pointScalarField& medialDist,
+                    const pointVectorField& medialVec,
 
                     List<extrudeMode>& extrudeStatus,
                     pointField& patchDisp,

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

@@ -689,7 +689,8 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
 
     pointVectorField& dispVec,
     pointScalarField& medialRatio,
-    pointScalarField& medialDist
+    pointScalarField& medialDist,
+    pointVectorField& medialVec
 ) const
 {
 
@@ -929,7 +930,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
         forAll(pointMedialDist, pointI)
         {
             medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
-            //medialVec[pointI] = pointMedialDist[pointI].origin();
+            medialVec[pointI] = pointMedialDist[pointI].origin();
         }
     }
 
@@ -966,14 +967,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
             << " : normalised direction of nearest displacement" << nl
             << "    " << medialDist.name()
             << " : distance to medial axis" << nl
-            //<< "    " << medialVec.name()
+            << "    " << medialVec.name()
-            //<< " : nearest point on medial axis" << nl
+            << " : nearest point on medial axis" << nl
             << "    " << medialRatio.name()
             << " : ratio of medial distance to wall distance" << nl
             << endl;
         dispVec.write();
         medialDist.write();
-        //medialVec.write();
+        medialVec.write();
         medialRatio.write();
     }
 }
@@ -996,7 +997,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
     const pointVectorField& dispVec,
     const pointScalarField& medialRatio,
     const pointScalarField& medialDist,
-    //const pointVectorField& medialVec,
+    const pointVectorField& medialVec,
 
     List<extrudeMode>& extrudeStatus,
     pointField& patchDisp,
@@ -1066,6 +1067,24 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
             << str().name() << endl;
     }
 
+    autoPtr<OFstream> medialVecStr;
+    label medialVertI = 0;
+    if (debug)
+    {
+        medialVecStr.reset
+        (
+            new OFstream
+            (
+                mesh.time().path()
+              / "thicknessRatioExcludeMedialVec_"
+              + meshRefiner_.timeName()
+              + ".obj"
+            )
+        );
+        Info<< "Writing points with too large a extrusion distance to "
+            << medialVecStr().name() << endl;
+    }
+
     forAll(meshPoints, patchPointI)
     {
         if (extrudeStatus[patchPointI] != NOEXTRUDE)
@@ -1082,12 +1101,9 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
             vector n =
                 patchDisp[patchPointI]
               / (mag(patchDisp[patchPointI]) + VSMALL);
-            //vector mVec = mesh.points()[pointI]-medialVec[pointI];
+            vector mVec = mesh.points()[pointI]-medialVec[pointI];
-            //scalar mDist = mag(mVec);
+            mVec /= mag(mVec)+VSMALL;
-            //scalar thicknessRatio =
+            thicknessRatio *= (n&mVec);
-            //    (n&mVec)
-            //   *thickness[patchPointI]
-            //   /(mDist+VSMALL);
 
             if (thicknessRatio > maxThicknessToMedialRatio)
             {
@@ -1103,8 +1119,9 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
                                 minThickness[patchPointI]
                                +thickness[patchPointI]
                             )
-                        //<< " since near medial at:" << medialVec[pointI]
+                        << " medial direction:" << mVec
-                        //<< " with thicknessRatio:" << thicknessRatio
+                        << " extrusion direction:" << n
+                        << " with thicknessRatio:" << thicknessRatio
                         << endl;
                 }
 
@@ -1124,6 +1141,16 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
                     vertI++;
                     str()<< "l " << vertI-1 << ' ' << vertI << nl;
                 }
+                if (medialVecStr.valid())
+                {
+                    const point& pt = mesh.points()[pointI];
+                    meshTools::writeOBJ(medialVecStr(), pt);
+                    medialVertI++;
+                    meshTools::writeOBJ(medialVecStr(), medialVec[pointI]);
+                    medialVertI++;
+                    medialVecStr()<< "l " << medialVertI-1
+                        << ' ' << medialVertI << nl;
+                }
             }
         }
     }
@@ -1228,6 +1255,15 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
                 )
             )()
         );
+
+        // Above move will have changed the instance only on the points (which
+        // is correct).
+        // However the previous mesh written will be the mesh with layers
+        // (see autoLayerDriver.C) so we now have to force writing all files
+        // so we can easily step through time steps. Note that if you
+        // don't write the mesh with layers this is not necessary.
+        meshRefiner_.mesh().setInstance(meshRefiner_.timeName());
+
         meshRefiner_.write
         (
             debug,