From e07309eea3ee8ee86a3a8986aaf5957c25cd461a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 14 Sep 2012 13:30:59 +0100
Subject: [PATCH] BUG: layering: disable mesh movement if extrusion disabled

---
 .../autoHexMeshDriver/autoLayerDriver.C       |  89 ++++++++++-
 .../autoHexMeshDriver/autoLayerDriver.H       |  14 +-
 .../autoHexMeshDriver/autoLayerDriverShrink.C | 139 +++++++++++++++---
 .../autoHexMeshDriver/autoSnapDriver.C        |  23 ++-
 .../autoHexMeshDriver/autoSnapDriverFeature.C |  70 +++++++++
 .../layerParameters/layerParameters.C         |   2 +-
 .../layerParameters/layerParameters.H         |   4 +
 .../meshRefinementProblemCells.C              |  13 +-
 8 files changed, 316 insertions(+), 38 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index abe3bc487c..26b4fb0187 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -48,6 +48,10 @@ Description
 #include "globalIndex.H"
 #include "DynamicField.H"
 #include "PatchTools.H"
+#include "slipPointPatchFields.H"
+#include "fixedValuePointPatchFields.H"
+#include "calculatedPointPatchFields.H"
+#include "cyclicSlipPointPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -71,7 +75,7 @@ void Foam::autoLayerDriver::dumpDisplacement
 )
 {
     OFstream dispStr(prefix + "_disp.obj");
-    Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
+    Info<< "Writing all displacements to " << dispStr.name() << endl;
 
     label vertI = 0;
 
@@ -87,7 +91,7 @@ void Foam::autoLayerDriver::dumpDisplacement
 
 
     OFstream illStr(prefix + "_illegal.obj");
-    Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
+    Info<< "Writing invalid displacements to " << illStr.name() << endl;
 
     vertI = 0;
 
@@ -801,6 +805,81 @@ void Foam::autoLayerDriver::setNumLayers
 }
 
 
+// Construct pointVectorField with correct boundary conditions for adding
+// layers
+Foam::tmp<Foam::pointVectorField>
+Foam::autoLayerDriver::makeLayerDisplacementField
+(
+    const pointMesh& pMesh,
+    const labelList& numLayers
+)
+{
+    // Construct displacement field.
+    const pointBoundaryMesh& pointPatches = pMesh.boundary();
+
+    wordList patchFieldTypes
+    (
+        pointPatches.size(),
+        slipPointPatchVectorField::typeName
+    );
+
+    forAll(numLayers, patchI)
+    {
+        //  0 layers: do not allow lslip so fixedValue 0
+        // >0 layers: fixedValue which gets adapted
+        if (numLayers[patchI] >= 0)
+        {
+            patchFieldTypes[patchI] = fixedValuePointPatchVectorField::typeName;
+        }
+    }
+
+    forAll(pointPatches, patchI)
+    {
+        if (isA<processorPointPatch>(pointPatches[patchI]))
+        {
+            patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
+        }
+        else if (isA<cyclicPointPatch>(pointPatches[patchI]))
+        {
+            patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName;
+        }
+    }
+
+
+//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl;
+//forAll(patchFieldTypes, patchI)
+//{
+//    Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name()
+//        << " type:" << patchFieldTypes[patchI]
+//        << " nLayers:" << numLayers[patchI]
+//        << endl;
+//}
+
+    const polyMesh& mesh = pMesh();
+
+    // Note: time().timeName() instead of meshRefinement::timeName() since
+    // postprocessable field.
+    tmp<pointVectorField> tfld
+    (
+        new pointVectorField
+        (
+            IOobject
+            (
+                "pointDisplacement",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            pMesh,
+            dimensionedVector("displacement", dimLength, vector::zero),
+            patchFieldTypes
+        )
+    );
+    return tfld;
+}
+
+
 void Foam::autoLayerDriver::growNoExtrusion
 (
     const indirectPrimitivePatch& pp,
@@ -2392,10 +2471,10 @@ void Foam::autoLayerDriver::addLayers
             mesh,
             pp(),
             patchIDs,
-            meshRefinement::makeDisplacementField
+            makeLayerDisplacementField
             (
                 pointMesh::New(mesh),
-                patchIDs
+                layerParams.numLayers()
             ),
             motionDict
         )
@@ -3186,7 +3265,7 @@ void Foam::autoLayerDriver::doLayers
     // Merge coplanar boundary faces
     mergePatchFacesUndo(layerParams, motionDict);
 
-    // Per patch the number of layers (0 if no layer)
+    // Per patch the number of layers (-1 or 0 if no layer)
     const labelList& numLayers = layerParams.numLayers();
 
     // Patches that need to get a layer
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index cb720dc64f..d36d15a4cb 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -195,6 +195,19 @@ class autoLayerDriver
                     List<extrudeMode>& extrudeStatus
                 ) const;
 
+                //- Helper function to make a pointVectorField with correct
+                //  bcs for layer addition:
+                //  - numLayers > 0         : fixedValue
+                //  - numLayers == 0        : fixedValue (always zero)
+                //  - processor             : calculated (so free to move)
+                //  - cyclic/wedge/symmetry : slip
+                //  - other                 : slip
+                static tmp<pointVectorField> makeLayerDisplacementField
+                (
+                    const pointMesh& pMesh,
+                    const labelList& numLayers
+                );
+
                 //- Grow no-extrusion layer.
                 void growNoExtrusion
                 (
@@ -444,7 +457,6 @@ class autoLayerDriver
                     const PackedBoolList& isMasterEdge,
                     const labelList& meshEdges,
                     const scalar minCosLayerTermination,
-                    scalarField& field,
                     List<extrudeMode>& extrudeStatus,
                     pointField& patchDisp,
                     labelList& patchNLayers
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index 35133f062d..58725470e1 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -99,6 +99,76 @@ void Foam::autoLayerDriver::sumWeights
 
 
 // Smooth field on moving patch
+//void Foam::autoLayerDriver::smoothField
+//(
+//    const motionSmoother& meshMover,
+//    const PackedBoolList& isMasterEdge,
+//    const labelList& meshEdges,
+//    const scalarField& fieldMin,
+//    const label nSmoothDisp,
+//    scalarField& field
+//) const
+//{
+//    const indirectPrimitivePatch& pp = meshMover.patch();
+//    const edgeList& edges = pp.edges();
+//    const labelList& meshPoints = pp.meshPoints();
+//
+//    scalarField invSumWeight(pp.nPoints());
+//    sumWeights
+//    (
+//        isMasterEdge,
+//        meshEdges,
+//        meshPoints,
+//        edges,
+//        invSumWeight
+//    );
+//
+//    // Get smoothly varying patch field.
+//    Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
+//
+//    for (label iter = 0; iter < nSmoothDisp; iter++)
+//    {
+//        scalarField average(pp.nPoints());
+//        averageNeighbours
+//        (
+//            meshMover.mesh(),
+//            isMasterEdge,
+//            meshEdges,
+//            meshPoints,
+//            pp.edges(),
+//            invSumWeight,
+//            field,
+//            average
+//        );
+//
+//        // Transfer to field
+//        forAll(field, pointI)
+//        {
+//            //full smoothing neighbours + point value
+//            average[pointI] = 0.5*(field[pointI]+average[pointI]);
+//
+//            // perform monotonic smoothing
+//            if
+//            (
+//                average[pointI] < field[pointI]
+//             && average[pointI] >= fieldMin[pointI]
+//            )
+//            {
+//                field[pointI] = average[pointI];
+//            }
+//        }
+//
+//        // Do residual calculation every so often.
+//        if ((iter % 10) == 0)
+//        {
+//            Info<< "    Iteration " << iter << "   residual "
+//                <<  gSum(mag(field-average))
+//                   /returnReduce(average.size(), sumOp<label>())
+//                << endl;
+//        }
+//    }
+//}
+//XXXXXXXXX
 void Foam::autoLayerDriver::smoothField
 (
     const motionSmoother& meshMover,
@@ -126,9 +196,15 @@ void Foam::autoLayerDriver::smoothField
     // Get smoothly varying patch field.
     Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
 
-    for (label iter = 0; iter < nSmoothDisp; iter++)
+
+    const scalar lambda = 0.33;
+    const scalar mu = -0.34;
+
+    for (label iter = 0; iter < 90; iter++)
     {
         scalarField average(pp.nPoints());
+
+        // Calculate average of field
         averageNeighbours
         (
             meshMover.mesh(),
@@ -141,23 +217,37 @@ void Foam::autoLayerDriver::smoothField
             average
         );
 
-        // Transfer to field
-        forAll(field, pointI)
+        forAll(field, i)
         {
-            //full smoothing neighbours + point value
-            average[pointI] = 0.5*(field[pointI]+average[pointI]);
-
-            // perform monotonic smoothing
-            if
-            (
-                average[pointI] < field[pointI]
-             && average[pointI] >= fieldMin[pointI]
-            )
+            if (field[i] >= fieldMin[i])
             {
-                field[pointI] = average[pointI];
+                field[i] = (1-lambda)*field[i]+lambda*average[i];
             }
         }
 
+
+        // Calculate average of field
+        averageNeighbours
+        (
+            meshMover.mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            pp.edges(),
+            invSumWeight,
+            field,
+            average
+        );
+
+        forAll(field, i)
+        {
+            if (field[i] >= fieldMin[i])
+            {
+                field[i] = (1-mu)*field[i]+mu*average[i];
+            }
+        }
+
+
         // Do residual calculation every so often.
         if ((iter % 10) == 0)
         {
@@ -168,7 +258,7 @@ void Foam::autoLayerDriver::smoothField
         }
     }
 }
-
+//XXXXXXXXX
 
 // Smooth normals on moving patch.
 void Foam::autoLayerDriver::smoothPatchNormals
@@ -480,7 +570,6 @@ void Foam::autoLayerDriver::findIsolatedRegions
     const PackedBoolList& isMasterEdge,
     const labelList& meshEdges,
     const scalar minCosLayerTermination,
-    scalarField& field,
     List<extrudeMode>& extrudeStatus,
     pointField& patchDisp,
     labelList& patchNLayers
@@ -661,7 +750,6 @@ void Foam::autoLayerDriver::findIsolatedRegions
                         )
                     )
                     {
-                        field[f[fp]] = 0.0;
                         nPointCounter++;
                     }
                 }
@@ -670,7 +758,8 @@ void Foam::autoLayerDriver::findIsolatedRegions
     }
 
     reduce(nPointCounter, sumOp<label>());
-    Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
+    Info<< "Number isolated points extrusion stopped : "<< nPointCounter
+        << endl;
 }
 
 
@@ -875,11 +964,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
+            //const pointPatchVectorField& pvf =
+            //    meshMover.displacement().boundaryField()[patchI];
 
             if
             (
                 !pp.coupled()
              && !isA<emptyPolyPatch>(pp)
+             //&&  pvf.constraintType() != word::null //exclude fixedValue
              && !adaptPatches.found(patchI)
             )
             {
@@ -1168,12 +1260,22 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         isMasterEdge,
         meshEdges,
         minCosLayerTermination,
-        thickness,
+
         extrudeStatus,
         patchDisp,
         patchNLayers
     );
 
+    // Update thickess for changed extrusion
+    forAll(thickness, patchPointI)
+    {
+        if (extrudeStatus[patchPointI] == NOEXTRUDE)
+        {
+            thickness[patchPointI] = 0.0;
+        }
+    }
+
+
     // smooth layer thickness on moving patch
     smoothField
     (
@@ -1182,6 +1284,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         meshEdges,
         minThickness,
         nSmoothThickness,
+
         thickness
     );
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 7f612246f4..edbd6a0474 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -506,7 +506,7 @@ void Foam::autoSnapDriver::dumpMove
 )
 {
     // Dump direction of growth into file
-    Pout<< nl << "Dumping move direction to " << fName << endl;
+    Info<< "Dumping move direction to " << fName << endl;
 
     OFstream nearestStream(fName);
 
@@ -721,14 +721,14 @@ void Foam::autoSnapDriver::preSmoothPatch
     if (debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing patch smoothed mesh to time "
+        Info<< "Writing patch smoothed mesh to time "
             << meshRefiner_.timeName() << '.' << endl;
         meshRefiner_.write
         (
             debug,
             mesh.time().path()/meshRefiner_.timeName()
         );
-        Pout<< "Dumped mesh in = "
+        Info<< "Dumped mesh in = "
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
     }
 
@@ -998,7 +998,7 @@ void Foam::autoSnapDriver::smoothDisplacement
     if (debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
+        Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
             << endl;
 
         // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
@@ -1010,13 +1010,12 @@ void Foam::autoSnapDriver::smoothDisplacement
             debug,
             mesh.time().path()/meshRefiner_.timeName()
         );
-
-        Pout<< "Writing displacement field ..." << endl;
+        Info<< "Writing displacement field ..." << endl;
         disp.write();
         tmp<pointScalarField> magDisp(mag(disp));
         magDisp().write();
 
-        Pout<< "Writing actual patch displacement ..." << endl;
+        Info<< "Writing actual patch displacement ..." << endl;
         vectorField actualPatchDisp(disp, pp.meshPoints());
         dumpMove
         (
@@ -1068,11 +1067,11 @@ bool Foam::autoSnapDriver::scaleMesh
         if (debug)
         {
             const_cast<Time&>(mesh.time())++;
-            Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName()
+            Info<< "Writing scaled mesh to time " << meshRefiner_.timeName()
                 << endl;
             mesh.write();
 
-            Pout<< "Writing displacement field ..." << endl;
+            Info<< "Writing displacement field ..." << endl;
             meshMover.displacement().write();
             tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
             magDisp().write();
@@ -1521,14 +1520,14 @@ void Foam::autoSnapDriver::doSnap
             if (debug)
             {
                 const_cast<Time&>(mesh.time())++;
-                Pout<< "Writing scaled mesh to time "
+                Info<< "Writing scaled mesh to time "
                     << meshRefiner_.timeName() << endl;
                 meshRefiner_.write
                 (
                     debug,
                     mesh.time().path()/meshRefiner_.timeName()
                 );
-                Pout<< "Writing displacement field ..." << endl;
+                Info<< "Writing displacement field ..." << endl;
                 meshMover.displacement().write();
                 tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
                 magDisp().write();
@@ -1564,7 +1563,7 @@ void Foam::autoSnapDriver::doSnap
     if (nChanged > 0 && debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing patchFace merged mesh to time "
+        Info<< "Writing patchFace merged mesh to time "
             << meshRefiner_.timeName() << endl;
         meshRefiner_.write
         (
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index e526f8faf9..88ec55b1c8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -2565,6 +2565,76 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         }
     }
 
+
+    // Collect additionally 'normal' boundary faces for boundaryPoints of pp
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // points on the boundary of pp should pick up non-pp normals
+    // as well for the feature-reconstruction to behave correctly.
+    // (the movement is already constrained outside correctly so it
+    //  is only that the unconstrained attraction vector is calculated
+    //  correctly)
+    {
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        labelList patchID(pbm.patchID());
+
+        // Unmark all non-coupled boundary faces
+        forAll(pbm, patchI)
+        {
+            const polyPatch& pp = pbm[patchI];
+
+            if (pp.coupled() || isA<emptyPolyPatch>(pp))
+            {
+                forAll(pp, i)
+                {
+                    label meshFaceI = pp.start()+i;
+                    patchID[meshFaceI-mesh.nInternalFaces()] = -1;
+                }
+            }
+        }
+
+        // Remove any meshed faces
+        PackedBoolList ppFaces(mesh.nFaces());
+        forAll(pp.addressing(), i)
+        {
+            label meshFaceI = pp.addressing()[i];
+            patchID[meshFaceI-mesh.nInternalFaces()] = -1;
+        }
+
+        // See if pp point uses any non-meshed boundary faces
+
+        const labelList& boundaryPoints = pp.boundaryPoints();
+        forAll(boundaryPoints, i)
+        {
+            label pointI = boundaryPoints[i];
+            label meshPointI = pp.meshPoints()[pointI];
+            const point& pt = mesh.points()[meshPointI];
+            const labelList& pFaces = mesh.pointFaces()[meshPointI];
+
+            List<point>& pNormals = pointFaceSurfNormals[pointI];
+            List<point>& pDisp = pointFaceDisp[pointI];
+            List<point>& pFc = pointFaceCentres[pointI];
+            labelList& pFid = pointFacePatchID[pointI];
+
+            forAll(pFaces, i)
+            {
+                label meshFaceI = pFaces[i];
+                if (!mesh.isInternalFace(meshFaceI))
+                {
+                    label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
+
+                    if (patchI != -1)
+                    {
+                        vector fn = mesh.faceAreas()[meshFaceI];
+                        pNormals.append(fn/mag(fn));
+                        pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
+                        pFc.append(mesh.faceCentres()[meshFaceI]);
+                        pFid.append(patchI);
+                    }
+                }
+            }
+        }
+    }
+
     syncTools::syncPointList
     (
         mesh,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index 95d9bf0dd4..cc466df4be 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -43,7 +43,7 @@ Foam::layerParameters::layerParameters
     const polyBoundaryMesh& boundaryMesh
 )
 :
-    numLayers_(boundaryMesh.size(), 0),
+    numLayers_(boundaryMesh.size(), -1),
     expansionRatio_
     (
         boundaryMesh.size(),
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
index 1a7913a8ef..4a48784ee9 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
@@ -131,6 +131,10 @@ public:
             // Per patch information
 
                 //- How many layers to add.
+                //  -1 : no specification. Assume 0 layers but allow sliding
+                //       to make layers
+                //   0 : specified to have 0 layers. No sliding allowed.
+                //  >0 : number of layers
                 const labelList& numLayers() const
                 {
                     return numLayers_;
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index f2ca0857dc..5561b21829 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -571,7 +571,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             hitInfo
         );
 
-        // Start of from current points
+        // Start off from current points
         newPoints = mesh_.points();
 
         forAll(hitInfo, i)
@@ -581,6 +581,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                 newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
             }
         }
+
+        if (debug)
+        {
+            const_cast<Time&>(mesh_.time())++;
+            pointField oldPoints(mesh_.points());
+            mesh_.movePoints(newPoints);
+            Pout<< "Writing newPoints mesh to time " << timeName()
+                << endl;
+            write(debug, mesh_.time().path()/"newPoints");
+            mesh_.movePoints(oldPoints);
+        }
     }