From f4109dc68b3ac67356605ce0bfe775b752714e36 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 26 Oct 2012 16:18:45 +0100
Subject: [PATCH 01/10] ENH: Further updates to fan and jump BCs

---
 .../derived/fan/fanFvPatchField.C             |  9 ++++
 .../derived/fan/fanFvPatchField.H             |  7 +++
 .../derived/fan/fanFvPatchFields.C            | 53 +++++++++++--------
 .../derived/fixedJump/fixedJumpFvPatchField.C |  5 +-
 .../fixedJumpAMI/fixedJumpAMIFvPatchField.C   |  3 +-
 5 files changed, 53 insertions(+), 24 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
index 3a150c6f70..18d5b33242 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
@@ -84,4 +84,13 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::fanFvPatchField<Type>::jump() const
+{
+    return uniformJumpFvPatchField<Type>::jump();
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
index 463ebeb543..bf0762f478 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@@ -172,6 +172,9 @@ public:
 
         // Evaluation functions
 
+            //- Return the "jump" across the patch
+            tmp<Field<Type> > jump() const;
+
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 };
@@ -185,6 +188,10 @@ fanFvPatchField<scalar>::fanFvPatchField
     const DimensionedField<scalar, volMesh>&,
     const dictionary&
 );
+
+template<>
+tmp<Field<scalar> > fanFvPatchField<scalar>::jump() const;
+
 template<>
 void fanFvPatchField<scalar>::updateCoeffs();
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
index 668d43483d..7922b66fe3 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
@@ -30,18 +30,16 @@ License
 #include "Tuple2.H"
 #include "polynomial.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-makeTemplatePatchTypeField
-(
-    fvPatchScalarField,
-    fanFvPatchScalarField
-);
+    makeTemplatePatchTypeField
+    (
+        fvPatchScalarField,
+        fanFvPatchScalarField
+    );
+}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -115,15 +113,10 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//- Specialisation of the jump-condition for the pressure
 template<>
-void Foam::fanFvPatchField<Foam::scalar>::updateCoeffs()
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::fanFvPatchField<Foam::scalar>::jump() const
 {
-    if (this->updated())
-    {
-        return;
-    }
-
     if (this->cyclicPatch().owner())
     {
         const surfaceScalarField& phi =
@@ -139,15 +132,33 @@ void Foam::fanFvPatchField<Foam::scalar>::updateCoeffs()
             Un /= patch().lookupPatchField<volScalarField, scalar>("rho");
         }
 
-        this->jump_ = this->jumpTable_->value(Un);
+        return this->jumpTable_->value(Un);
+    }
+    else
+    {
+        return refCast<const fanFvPatchField<scalar> >
+        (
+            this->neighbourPatchField()
+        ).jump();
+    }
+}
+
+
+template<>
+void Foam::fanFvPatchField<Foam::scalar>::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    if (this->cyclicPatch().owner())
+    {
+        this->jump_ = jump();
     }
 
     uniformJumpFvPatchField<scalar>::updateCoeffs();
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
index e58270ec63..9ef4be21a4 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
@@ -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
@@ -138,7 +138,8 @@ template<class Type>
 void Foam::fixedJumpFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
-    os.writeKeyword("patchType") << "cyclic" << token::END_STATEMENT << nl;
+    os.writeKeyword("patchType") << this->interfaceFieldType()
+        << token::END_STATEMENT << nl;
     jump_.writeEntry("jump", os);
     this->writeEntry("value", os);
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
index d9cdc383a3..4c3290cf2d 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
@@ -141,7 +141,8 @@ template<class Type>
 void Foam::fixedJumpAMIFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
-    os.writeKeyword("patchType") << "cyclicAMI" << token::END_STATEMENT << nl;
+    os.writeKeyword("patchType") << this->interfaceFieldType()
+        << token::END_STATEMENT << nl;
     jump_.writeEntry("jump", os);
     this->writeEntry("value", os);
 }

From 725a1b3eb03759fefbd2fceae1a3b68edb19b4fd Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 26 Oct 2012 16:47:04 +0100
Subject: [PATCH 02/10] BUG: snappyHexMesh: multi-patch from surrounding
 regions, not raw patchIDs

---
 .../autoHexMeshDriver/autoSnapDriver.C        |  11 +-
 .../autoHexMeshDriver/autoSnapDriver.H        |  26 +-
 .../autoHexMeshDriver/autoSnapDriverFeature.C | 510 +++++++++---------
 .../snapParameters/snapParameters.C           |   6 +-
 .../snapParameters/snapParameters.H           |   7 +
 5 files changed, 283 insertions(+), 277 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index e358ba4ffa..54a78f477c 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1406,9 +1406,10 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
+        indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
+        const scalarField snapDist(calcSnapDistance(snapParams, pp));
 
 
         // Construct iterative mesh mover.
@@ -1420,7 +1421,7 @@ void Foam::autoSnapDriver::doSnap
         motionSmoother meshMover
         (
             mesh,
-            ppPtr(),
+            pp,
             adaptPatchIDs,
             meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
             motionDict
@@ -1475,7 +1476,7 @@ void Foam::autoSnapDriver::doSnap
             }
 
             // Check for displacement being outwards.
-            outwardsDisplacement(ppPtr(), disp);
+            outwardsDisplacement(pp, disp);
 
             // Set initial distribution of displacement field (on patches)
             // from patchDisp and make displacement consistent with b.c.
@@ -1489,8 +1490,8 @@ void Foam::autoSnapDriver::doSnap
                 (
                     mesh.time().path()
                   / "patchDisplacement_" + name(iter) + ".obj",
-                    ppPtr().localPoints(),
-                    ppPtr().localPoints() + disp
+                    pp.localPoints(),
+                    pp.localPoints() + disp
                 );
             }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index bd5dd86188..7cf1e79f61 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -141,14 +141,6 @@ class autoSnapDriver
                     vectorField& pointSurfaceNormal,
                     vectorField& pointRotation
                 ) const;
-                //void calcNearestFace
-                //(
-                //    const label iter,
-                //    const indirectPrimitivePatch& pp,
-                //    vectorField& faceDisp,
-                //    vectorField& faceSurfaceNormal,
-                //    vectorField& faceRotation
-                //) const;
                 void calcNearestFace
                 (
                     const label iter,
@@ -158,15 +150,19 @@ class autoSnapDriver
                     labelList& faceSurfaceRegion,
                     vectorField& faceRotation
                 ) const;
-                void interpolateFaceToPoint
+                void calcNearestFacePointProperties
                 (
                     const label iter,
                     const indirectPrimitivePatch& pp,
-                    const List<List<point> >& pointFaceDisp,
-                    const List<List<point> >& pointFaceRotation,
-                    const List<List<point> >& pointFaceCentres,
-                    vectorField& patchDisp
-                    //vectorField& patchRotationDisp
+
+                    const vectorField& faceDisp,
+                    const vectorField& faceSurfaceNormal,
+                    const labelList& faceSurfaceRegion,
+
+                    List<List<point> >& pointFaceSurfNormals,
+                    List<List<point> >& pointFaceDisp,
+                    List<List<point> >& pointFaceCentres,
+                    List<labelList>&    pointFacePatchID
                 ) const;
                 void correctAttraction
                 (
@@ -276,6 +272,7 @@ class autoSnapDriver
                 (
                     const label iter,
                     const scalar featureCos,
+                    const bool multiRegionFeatureSnap,
 
                     const indirectPrimitivePatch&,
                     const scalarField&,
@@ -339,6 +336,7 @@ class autoSnapDriver
                 (
                     const label iter,
                     const scalar featureCos,
+                    const bool multiRegionFeatureSnap,
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 8bcd9c6270..f32d9ce4bd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -321,7 +321,7 @@ void Foam::autoSnapDriver::calcNearestFace
     const indirectPrimitivePatch& pp,
     vectorField& faceDisp,
     vectorField& faceSurfaceNormal,
-    labelList& faceSurfaceRegion,
+    labelList& faceSurfaceGlobalRegion,
     vectorField& faceRotation
 ) const
 {
@@ -333,8 +333,8 @@ void Foam::autoSnapDriver::calcNearestFace
     faceDisp = vector::zero;
     faceSurfaceNormal.setSize(pp.size());
     faceSurfaceNormal = vector::zero;
-    faceSurfaceRegion.setSize(pp.size());
-    faceSurfaceRegion = -1;
+    faceSurfaceGlobalRegion.setSize(pp.size());
+    faceSurfaceGlobalRegion = -1;
 
     // Divide surfaces into zoned and unzoned
     labelList zonedSurfaces = surfaces.getNamedSurfaces();
@@ -419,7 +419,7 @@ void Foam::autoSnapDriver::calcNearestFace
                 label faceI = ppFaces[hitI];
                 faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
                 faceSurfaceNormal[faceI] = hitNormal[hitI];
-                faceSurfaceRegion[faceI] = surfaces.globalRegion
+                faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
                 (
                     hitSurface[hitI],
                     hitRegion[hitI]
@@ -477,7 +477,11 @@ void Foam::autoSnapDriver::calcNearestFace
             label faceI = ppFaces[hitI];
             faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
             faceSurfaceNormal[faceI] = hitNormal[hitI];
-            faceSurfaceRegion[faceI] = hitRegion[hitI];
+            faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
+            (
+                hitSurface[hitI],
+                hitRegion[hitI]
+            );
         }
     }
 
@@ -517,6 +521,169 @@ void Foam::autoSnapDriver::calcNearestFace
 }
 
 
+// Collect (possibly remote) per point data of all surrounding faces
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// - faceSurfaceNormal
+// - faceDisp
+// - faceCentres&faceNormal
+void Foam::autoSnapDriver::calcNearestFacePointProperties
+(
+    const label iter,
+    const indirectPrimitivePatch& pp,
+
+    const vectorField& faceDisp,
+    const vectorField& faceSurfaceNormal,
+    const labelList& faceSurfaceGlobalRegion,
+
+    List<List<point> >& pointFaceSurfNormals,
+    List<List<point> >& pointFaceDisp,
+    List<List<point> >& pointFaceCentres,
+    List<labelList>&    pointFacePatchID
+) const
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    // For now just get all surrounding face data. Expensive - should just
+    // store and sync data on coupled points only
+    // (see e.g PatchToolsNormals.C)
+
+    pointFaceSurfNormals.setSize(pp.nPoints());
+    pointFaceDisp.setSize(pp.nPoints());
+    pointFaceCentres.setSize(pp.nPoints());
+    pointFacePatchID.setSize(pp.nPoints());
+
+    // Fill local data
+    forAll(pp.pointFaces(), pointI)
+    {
+        const labelList& pFaces = pp.pointFaces()[pointI];
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
+        pNormals.setSize(pFaces.size());
+        List<point>& pDisp = pointFaceDisp[pointI];
+        pDisp.setSize(pFaces.size());
+        List<point>& pFc = pointFaceCentres[pointI];
+        pFc.setSize(pFaces.size());
+        labelList& pFid = pointFacePatchID[pointI];
+        pFid.setSize(pFaces.size());
+
+        forAll(pFaces, i)
+        {
+            label faceI = pFaces[i];
+            pNormals[i] = faceSurfaceNormal[faceI];
+            pDisp[i] = faceDisp[faceI];
+            pFc[i] = pp.faceCentres()[faceI];
+            //label meshFaceI = pp.addressing()[faceI];
+            //pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
+            pFid[i] = globalToPatch_[faceSurfaceGlobalRegion[faceI]];
+        }
+    }
+
+
+    // 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
+        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,
+        pp.meshPoints(),
+        pointFaceSurfNormals,
+        listPlusEqOp<point>(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceDisp,
+        listPlusEqOp<point>(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFaceCentres,
+        listPlusEqOp<point>(),
+        List<point>(),
+        listTransform()
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        pointFacePatchID,
+        listPlusEqOp<label>(),
+        List<label>()
+    );
+}
+
+
 // Gets passed in offset to nearest point on feature edge. Calculates
 // if the point has a different number of faces on either side of the feature
 // and if so attracts the point to that non-dominant plane.
@@ -580,56 +747,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findMultiPatchPoint
     }
     return pointIndexHit(false, vector::zero, labelMax);
 }
-////XXXXXXXX
-//void Foam::autoSnapDriver::attractMultiPatchPoint
-//(
-//    const label iter,
-//    const scalar featureCos,
-//
-//    const indirectPrimitivePatch& pp,
-//    const scalarField& snapDist,
-//    const label pointI,
-//
-//    const List<List<point> >& pointFaceSurfNormals,
-//    const labelListList& pointFaceSurfaceRegion,
-//    const List<List<point> >& pointFaceDisp,
-//    const List<List<point> >& pointFaceCentres,
-//    const labelListList& pointFacePatchID,
-//
-//    vector& patchAttraction,
-//    pointConstraint& patchConstraint
-//) const
-//{
-//    // Collect
-//
-//        );
-//
-//        if
-//        (
-//            (constraint.first() > patchConstraints[pointI].first())
-//         || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
-//        )
-//        {
-//            patchAttraction[pointI] = attraction;
-//            patchConstraints[pointI] = constraint;
-//
-//            // Check the number of directions
-//            if (patchConstraints[pointI].first() == 1)
-//            {
-//                // Flat surface. Check for different patchIDs
-//                pointIndexHit multiPatchPt
-//                (
-//                    findMultiPatchPoint
-//                    (
-//                        pt,
-//                        pointFacePatchID[pointI],
-//                        pointFaceCentres[pointI]
-//                    )
-//                );
-//                if (multiPatchPt.hit())
-//                {
-//                    // Behave like when having two surface normals so
-////XXXXXXXX
+
 
 void Foam::autoSnapDriver::binFeatureFace
 (
@@ -1361,6 +1479,7 @@ void Foam::autoSnapDriver::determineFeatures
 (
     const label iter,
     const scalar featureCos,
+    const bool multiRegionFeatureSnap,
 
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
@@ -1464,69 +1583,72 @@ void Foam::autoSnapDriver::determineFeatures
             if (patchConstraints[pointI].first() == 1)
             {
                 // Flat surface. Check for different patchIDs
-                pointIndexHit multiPatchPt
-                (
-                    findMultiPatchPoint
-                    (
-                        pt,
-                        pointFacePatchID[pointI],
-                        pointFaceCentres[pointI]
-                    )
-                );
-                if (multiPatchPt.hit())
+                if (multiRegionFeatureSnap)
                 {
-                    // Behave like when having two surface normals so
-                    // attract to nearest feature edge (with a guess for
-                    // the multipatch point as starting point)
-                    label featI = -1;
-                    pointIndexHit nearInfo = findNearFeatureEdge
+                    pointIndexHit multiPatchPt
                     (
-                        pp,
-                        snapDist,
-                        pointI,
-                        multiPatchPt.hitPoint(),        //estimatedPt
-
-                        featI,
-                        edgeAttractors,
-                        edgeConstraints,
-
-                        patchAttraction,
-                        patchConstraints
+                        findMultiPatchPoint
+                        (
+                            pt,
+                            pointFacePatchID[pointI],
+                            pointFaceCentres[pointI]
+                        )
                     );
+                    if (multiPatchPt.hit())
+                    {
+                        // Behave like when having two surface normals so
+                        // attract to nearest feature edge (with a guess for
+                        // the multipatch point as starting point)
+                        label featI = -1;
+                        pointIndexHit nearInfo = findNearFeatureEdge
+                        (
+                            pp,
+                            snapDist,
+                            pointI,
+                            multiPatchPt.hitPoint(),        //estimatedPt
 
-                    if (nearInfo.hit())
-                    {
-                        // Dump
-                        if (featureEdgeStr.valid())
+                            featI,
+                            edgeAttractors,
+                            edgeConstraints,
+
+                            patchAttraction,
+                            patchConstraints
+                        );
+
+                        if (nearInfo.hit())
                         {
-                            meshTools::writeOBJ(featureEdgeStr(), pt);
-                            featureEdgeVertI++;
-                            meshTools::writeOBJ
-                            (
-                                featureEdgeStr(),
-                                nearInfo.hitPoint()
-                            );
-                            featureEdgeVertI++;
-                            featureEdgeStr()
-                                << "l " << featureEdgeVertI-1 << ' '
-                                << featureEdgeVertI << nl;
+                            // Dump
+                            if (featureEdgeStr.valid())
+                            {
+                                meshTools::writeOBJ(featureEdgeStr(), pt);
+                                featureEdgeVertI++;
+                                meshTools::writeOBJ
+                                (
+                                    featureEdgeStr(),
+                                    nearInfo.hitPoint()
+                                );
+                                featureEdgeVertI++;
+                                featureEdgeStr()
+                                    << "l " << featureEdgeVertI-1 << ' '
+                                    << featureEdgeVertI << nl;
+                            }
                         }
-                    }
-                    else
-                    {
-                        if (missedEdgeStr.valid())
+                        else
                         {
-                            meshTools::writeOBJ(missedEdgeStr(), pt);
-                            missedVertI++;
-                            meshTools::writeOBJ
-                            (
-                                missedEdgeStr(),
-                                nearInfo.missPoint()
-                            );
-                            missedVertI++;
-                            missedEdgeStr()
-                                << "l " << missedVertI-1 << ' '
-                                << missedVertI << nl;
+                            if (missedEdgeStr.valid())
+                            {
+                                meshTools::writeOBJ(missedEdgeStr(), pt);
+                                missedVertI++;
+                                meshTools::writeOBJ
+                                (
+                                    missedEdgeStr(),
+                                    nearInfo.missPoint()
+                                );
+                                missedVertI++;
+                                missedEdgeStr()
+                                    << "l " << missedVertI-1 << ' '
+                                    << missedVertI << nl;
+                            }
                         }
                     }
                 }
@@ -1645,6 +1767,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 (
     const label iter,
     const scalar featureCos,
+    const bool multiRegionFeatureSnap,
 
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
@@ -1697,6 +1820,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     (
         iter,
         featureCos,
+        multiRegionFeatureSnap,
 
         pp,
         snapDist,
@@ -2131,6 +2255,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     //MEJ: any faces that have multi-patch points only keep the multi-patch
     //     points. The other points on the face will be dragged along
     //     (hopefully)
+    if (multiRegionFeatureSnap)
     {
         autoPtr<OFstream> multiPatchStr;
         if (debug&meshRefinement::OBJINTERSECTIONS)
@@ -2536,10 +2661,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 {
     const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
     const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
+    const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
 
     Info<< "Overriding displacement on features :" << nl
-        << "   implicit features : " << implicitFeatureAttraction << nl
-        << "   explicit features : " << explicitFeatureAttraction << nl
+        << "   implicit features    : " << implicitFeatureAttraction << nl
+        << "   explicit features    : " << explicitFeatureAttraction << nl
+        << "   multi-patch features : " << multiRegionFeatureSnap << nl
         << endl;
 
 
@@ -2547,6 +2674,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const pointField& localPoints = pp.localPoints();
     const fvMesh& mesh = meshRefiner_.mesh();
 
+
     // Displacement and orientation per pp face
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -2554,7 +2682,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     vectorField faceDisp(pp.size(), vector::zero);
     // normal of surface at point on surface
     vectorField faceSurfaceNormal(pp.size(), vector::zero);
-    labelList faceSurfaceRegion(pp.size(), -1);
+    labelList faceSurfaceGlobalRegion(pp.size(), -1);
     vectorField faceRotation(pp.size(), vector::zero);
 
     calcNearestFace
@@ -2563,7 +2691,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pp,
         faceDisp,
         faceSurfaceNormal,
-        faceSurfaceRegion,
+        faceSurfaceGlobalRegion,
         faceRotation
     );
 
@@ -2587,158 +2715,25 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
-    // - faceRotation
     // - faceCentres&faceNormal
-
-    // For now just get all surrounding face data. Expensive - should just
-    // store and sync data on coupled points only
-    // (see e.g PatchToolsNormals.C)
-
     List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
-    //List<List<point> > pointFaceRotation(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
     List<labelList>    pointFacePatchID(pp.nPoints());
 
-    // Fill local data
-    forAll(pp.pointFaces(), pointI)
-    {
-        const labelList& pFaces = pp.pointFaces()[pointI];
-        List<point>& pNormals = pointFaceSurfNormals[pointI];
-        pNormals.setSize(pFaces.size());
-        List<point>& pDisp = pointFaceDisp[pointI];
-        pDisp.setSize(pFaces.size());
-        //List<point>& pRot = pointFaceRotation[pointI];
-        //pRot.setSize(pFaces.size());
-        List<point>& pFc = pointFaceCentres[pointI];
-        pFc.setSize(pFaces.size());
-        labelList& pFid = pointFacePatchID[pointI];
-        pFid.setSize(pFaces.size());
-
-        forAll(pFaces, i)
-        {
-            label faceI = pFaces[i];
-            pNormals[i] = faceSurfaceNormal[faceI];
-            pDisp[i] = faceDisp[faceI];
-            //pRot[i] = faceRotation[faceI];
-            pFc[i] = pp.faceCentres()[faceI];
-            label meshFaceI = pp.addressing()[faceI];
-            pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
-        }
-    }
-
-
-    // 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
-        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
+    calcNearestFacePointProperties
     (
-        mesh,
-        pp.meshPoints(),
+        iter,
+        pp,
+
+        faceDisp,
+        faceSurfaceNormal,
+        faceSurfaceGlobalRegion,
+
         pointFaceSurfNormals,
-        listPlusEqOp<point>(),
-        List<point>(),
-        listTransform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
         pointFaceDisp,
-        listPlusEqOp<point>(),
-        List<point>(),
-        listTransform()
-    );
-    //syncTools::syncPointList
-    //(
-    //    mesh,
-    //    pp.meshPoints(),
-    //    pointFaceRotation,
-    //    listPlusEqOp<point>(),
-    //    List<point>(),
-    //    listTransform()
-    //);
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
         pointFaceCentres,
-        listPlusEqOp<point>(),
-        List<point>(),
-        listTransform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        pointFacePatchID,
-        listPlusEqOp<label>(),
-        List<label>()
+        pointFacePatchID
     );
 
 
@@ -2793,6 +2788,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         (
             iter,
             featureCos,
+            multiRegionFeatureSnap,
 
             pp,
             snapDist,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
index b9d6f3b0c6..77bff60187 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
@@ -36,7 +36,11 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
     nSnap_(readLabel(dict.lookup("nRelaxIter"))),
     nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
     explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
-    implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false))
+    implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false)),
+    multiRegionFeatureSnap_
+    (
+        dict.lookupOrDefault("multiRegionFeatureSnap", false)
+    )
 {}
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
index 8103cd4082..294d0de55d 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
@@ -68,6 +68,8 @@ class snapParameters
 
         const Switch implicitFeatureSnap_;
 
+        const Switch multiRegionFeatureSnap_;
+
 
     // Private Member Functions
 
@@ -134,6 +136,11 @@ public:
                 return implicitFeatureSnap_;
             }
 
+            Switch multiRegionFeatureSnap() const
+            {
+                return multiRegionFeatureSnap_;
+            }
+
 };
 
 

From ee3fd4c4800aa89d3a77db3835c26bc87be524a4 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Sat, 27 Oct 2012 21:34:09 +0100
Subject: [PATCH 03/10] Update $* to "$@"

---
 bin/foamNew | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/bin/foamNew b/bin/foamNew
index 4c5a996aad..5d61de6333 100755
--- a/bin/foamNew
+++ b/bin/foamNew
@@ -52,11 +52,11 @@ USAGE
 case "$1" in
 (-s | -source | source)
     shift
-    $WM_PROJECT_DIR/etc/codeTemplates/source/foamNewSource $*
+    $WM_PROJECT_DIR/etc/codeTemplates/source/foamNewSource "$@"
     ;;
 (-t | -template | template)
     shift
-    $WM_PROJECT_DIR/etc/codeTemplates/template/foamNewTemplate $*
+    $WM_PROJECT_DIR/etc/codeTemplates/template/foamNewTemplate "$@"
     ;;
 (*)
     usage "unknown type '$1'"

From 3a1113e51331bcf5cbd650d0deb93223f891f3fc Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Sat, 27 Oct 2012 21:34:40 +0100
Subject: [PATCH 04/10] Updated to avoid creating PINC and PLIBS environment
 variables

---
 etc/config/settings.csh | 12 +-----------
 etc/config/settings.sh  | 13 +------------
 2 files changed, 2 insertions(+), 23 deletions(-)

diff --git a/etc/config/settings.csh b/etc/config/settings.csh
index 069624f513..059ea30e33 100644
--- a/etc/config/settings.csh
+++ b/etc/config/settings.csh
@@ -391,22 +391,12 @@ case SYSTEMOPENMPI:
     # Use the system installed openmpi, get library directory via mpicc
     setenv FOAM_MPI openmpi-system
 
-    # Set compilation flags here instead of in wmake/rules/../mplibSYSTEMOPENMPI
-    setenv PINC "`mpicc --showme:compile`"
-    setenv PLIBS "`mpicc --showme:link`"
-    set libDir=`echo "$PLIBS" | sed -e 's/.*-L\([^ ]*\).*/\1/'`
+    set libDir=`mpicc --showme:link | sed -e 's/.*-L\([^ ]*\).*/\1/'`
 
     # Bit of a hack: strip off 'lib' and hope this is the path to openmpi
     # include files and libraries.
     setenv MPI_ARCH_PATH "${libDir:h}"
 
-    if ($?FOAM_VERBOSE && $?prompt) then
-        echo "Using system installed MPI:"
-        echo "    compile flags : $PINC"
-        echo "    link flags    : $PLIBS"
-        echo "    libmpi dir    : $libDir"
-    endif
-
     _foamAddLib     $libDir
     unset libDir
     breaksw
diff --git a/etc/config/settings.sh b/etc/config/settings.sh
index 1fc2136530..ed45429a65 100644
--- a/etc/config/settings.sh
+++ b/etc/config/settings.sh
@@ -418,23 +418,12 @@ SYSTEMOPENMPI)
     # Use the system installed openmpi, get library directory via mpicc
     export FOAM_MPI=openmpi-system
 
-    # Set compilation flags here instead of in wmake/rules/../mplibSYSTEMOPENMPI
-    export PINC="`mpicc --showme:compile`"
-    export PLIBS="`mpicc --showme:link`"
-    libDir=`echo "$PLIBS" | sed -e 's/.*-L\([^ ]*\).*/\1/'`
+    libDir=`mpicc --showme:link | sed -e 's/.*-L\([^ ]*\).*/\1/'`
 
     # Bit of a hack: strip off 'lib' and hope this is the path to openmpi
     # include files and libraries.
     export MPI_ARCH_PATH="${libDir%/*}"
 
-    if [ "$FOAM_VERBOSE" -a "$PS1" ]
-    then
-        echo "Using system installed MPI:" 1>&2
-        echo "    compile flags : $PINC" 1>&2
-        echo "    link flags    : $PLIBS" 1>&2
-        echo "    libmpi dir    : $libDir" 1>&2
-    fi
-
     _foamAddLib     $libDir
     unset libDir
     ;;

From 0ef5b1cbfbf2c9dabaafd4bbd06954561790e8b7 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 29 Oct 2012 12:22:04 +0000
Subject: [PATCH 05/10] ENH: Further updates to fan and jump BCs

---
 .../jumpCyclic/jumpCyclicFvPatchField.C       | 52 +++++-------------
 .../jumpCyclic/jumpCyclicFvPatchField.H       | 11 ++++
 .../jumpCyclic/jumpCyclicFvPatchFields.C      | 54 ++++++++++++++++++-
 .../jumpCyclicAMI/jumpCyclicAMIFvPatchField.C | 49 ++++++-----------
 .../jumpCyclicAMI/jumpCyclicAMIFvPatchField.H | 12 +++++
 .../jumpCyclicAMIFvPatchFields.C              | 44 +++++++++++++++
 .../derived/fan/fanFvPatchField.C             |  7 ---
 .../derived/fan/fanFvPatchField.H             | 11 +---
 .../derived/fan/fanFvPatchFields.C            |  5 --
 .../uniformJump/uniformJumpFvPatchField.C     |  2 +-
 10 files changed, 153 insertions(+), 94 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
index db6849d267..2287477bdf 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
@@ -142,43 +142,17 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes
 ) const
 {
-    scalarField pnf(this->size());
-
-    const labelUList& nbrFaceCells =
-        this->cyclicPatch().neighbFvPatch().faceCells();
-
-    // for AMG solve - only apply jump to finest level
-    if (psiInternal.size() == this->internalField().size())
-    {
-        Field<scalar> jf(this->jump()().component(cmpt));
-
-        if (!this->cyclicPatch().owner())
-        {
-            jf *= -1.0;
-        }
-
-        forAll(*this, facei)
-        {
-            pnf[facei] = psiInternal[nbrFaceCells[facei]] - jf[facei];
-        }
-    }
-    else
-    {
-        forAll(*this, facei)
-        {
-            pnf[facei] = psiInternal[nbrFaceCells[facei]];
-        }
-    }
-
-    // Transform according to the transformation tensors
-    this->transformCoupleField(pnf, cmpt);
-
-    // Multiply the field by coefficients and add into the result
-    const labelUList& faceCells = this->cyclicPatch().faceCells();
-    forAll(faceCells, elemI)
-    {
-        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
-    }
+    notImplemented
+    (
+        "void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix"
+        "("
+            "scalarField&, "
+            "const scalarField&, "
+            "const scalarField&, "
+            "const direction, "
+            "const Pstream::commsTypes"
+        ") const"
+    );
 }
 
 
@@ -196,8 +170,8 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
     const labelUList& nbrFaceCells =
         this->cyclicPatch().neighbFvPatch().faceCells();
 
-    // for AMG solve - only apply jump to finest level
-    if (psiInternal.size() == this->internalField().size())
+    // only apply jump to original field
+    if (&psiInternal == &this->internalField())
     {
         Field<Type> jf(this->jump());
 
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H
index a5fd2343bf..42d0c1805f 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H
@@ -144,6 +144,17 @@ public:
             ) const;
 };
 
+//- Update result field based on interface functionality
+template<>
+void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
+(
+    scalarField& result,
+    const scalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes commsType
+) const;
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
index 02e25e31da..adfc1326ff 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
@@ -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
@@ -36,6 +36,58 @@ namespace Foam
 
 makePatchFieldsTypeName(jumpCyclic);
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<>
+void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
+(
+    scalarField& result,
+    const scalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes
+) const
+{
+    scalarField pnf(this->size());
+
+    const labelUList& nbrFaceCells =
+        this->cyclicPatch().neighbFvPatch().faceCells();
+
+    // only apply jump to original field
+    if (&psiInternal == &this->internalField())
+    {
+        Field<scalar> jf(this->jump());
+
+        if (!this->cyclicPatch().owner())
+        {
+            jf *= -1.0;
+        }
+
+        forAll(*this, facei)
+        {
+            pnf[facei] = psiInternal[nbrFaceCells[facei]] - jf[facei];
+        }
+    }
+    else
+    {
+        forAll(*this, facei)
+        {
+            pnf[facei] = psiInternal[nbrFaceCells[facei]];
+        }
+    }
+
+    // Transform according to the transformation tensors
+    this->transformCoupleField(pnf, cmpt);
+
+    // Multiply the field by coefficients and add into the result
+    const labelUList& faceCells = this->cyclicPatch().faceCells();
+    forAll(faceCells, elemI)
+    {
+        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C
index 87d7c20efd..ff36ee9797 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C
@@ -129,33 +129,17 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes
 ) const
 {
-    const labelUList& nbrFaceCells =
-        this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
-
-    scalarField pnf(psiInternal, nbrFaceCells);
-
-    pnf = this->cyclicAMIPatch().interpolate(pnf);
-
-    // for AMG solve - only apply jump to finest level
-    if (psiInternal.size() == this->internalField().size())
-    {
-        tmp<Field<scalar> > tjf = jump()().component(cmpt);
-        if (!this->cyclicAMIPatch().owner())
-        {
-            tjf = -tjf;
-        }
-        pnf -= tjf;
-    }
-
-    // Transform according to the transformation tensors
-    this->transformCoupleField(pnf, cmpt);
-
-    // Multiply the field by coefficients and add into the result
-    const labelUList& faceCells = this->cyclicAMIPatch().faceCells();
-    forAll(faceCells, elemI)
-    {
-        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
-    }
+    notImplemented
+    (
+        "void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix"
+        "("
+            "scalarField&, "
+            "const scalarField&, "
+            "const scalarField& coeffs,"
+            "const direction, "
+            "const Pstream::commsTypes"
+        ") const"
+    );
 }
 
 
@@ -175,15 +159,16 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 
     pnf = this->cyclicAMIPatch().interpolate(pnf);
 
-    // for AMG solve - only apply jump to finest level
-    if (psiInternal.size() == this->internalField().size())
+    // only apply jump to original field
+    if (&psiInternal == &this->internalField())
     {
-        tmp<Field<Type> > tjf = jump();
+        Field<Type> jf(this->jump());
         if (!this->cyclicAMIPatch().owner())
         {
-            tjf = -tjf;
+            jf *= -1.0;
         }
-        pnf -= tjf;
+
+        pnf -= jf;
     }
 
     // Transform according to the transformation tensors
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H
index 6759b09a36..a2e77f0396 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H
@@ -148,6 +148,18 @@ public:
 };
 
 
+//- Update result field based on interface functionality
+template<>
+void jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
+(
+    scalarField& result,
+    const scalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes commsType
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C
index 323f460d9f..b7349812da 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C
@@ -36,6 +36,50 @@ namespace Foam
 
 makePatchFieldsTypeName(jumpCyclicAMI);
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<>
+void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
+(
+    scalarField& result,
+    const scalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes
+) const
+{
+    const labelUList& nbrFaceCells =
+        this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
+
+    scalarField pnf(psiInternal, nbrFaceCells);
+
+    pnf = this->cyclicAMIPatch().interpolate(pnf);
+
+    // only apply jump to original field
+    if (&psiInternal == &this->internalField())
+    {
+        Field<scalar> jf(this->jump());
+
+        if (!this->cyclicAMIPatch().owner())
+        {
+            jf *= -1.0;
+        }
+
+        pnf -= jf;
+    }
+
+    // Transform according to the transformation tensors
+    this->transformCoupleField(pnf, cmpt);
+
+    // Multiply the field by coefficients and add into the result
+    const labelUList& faceCells = this->cyclicAMIPatch().faceCells();
+    forAll(faceCells, elemI)
+    {
+        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
index d397dbbc5e..775c0f1ced 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
@@ -98,13 +98,6 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::tmp<Foam::Field<Type> > Foam::fanFvPatchField<Type>::jump() const
-{
-    return this->jump_;
-}
-
-
 template<class Type>
 void Foam::fanFvPatchField<Type>::updateCoeffs()
 {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
index a38a5d9cf8..218b4b1862 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@@ -39,7 +39,6 @@ Description
     \table
         Property     | Description             | Required    | Default value
         patchType    | underlying patch type should be \c cyclic| yes |
-        jump         | current jump value      | yes         |
         jumpTable    | jump data, e.g. \c csvFile | yes      |
     \endtable
 
@@ -49,7 +48,6 @@ Description
     {
         type            fan;
         patchType       cyclic;
-        jump            uniform 0;
         jumpTable       csvFile;
         csvFileCoeffs
         {
@@ -176,13 +174,8 @@ public:
 
     // Member functions
 
-        // Evaluation functions
-
-            //- Return the "jump" across the patch
-            tmp<Field<Type> > jump() const;
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
 };
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
index 5d99097dc1..8535c57127 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
@@ -116,11 +116,6 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
             this->jumpTable_ = DataEntry<scalar>::New("jumpTable", dict);
         }
     }
-    else
-    {
-        // Dummy jump table
-        this->jumpTable_.reset(new DataEntry<scalar>("jumpTable"));
-    }
 
     if (dict.found("value"))
     {
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C
index 1647c3b257..250442dad0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C
@@ -35,7 +35,7 @@ Foam::uniformJumpFvPatchField<Type>::uniformJumpFvPatchField
 )
 :
     fixedJumpFvPatchField<Type>(p, iF),
-    jumpTable_(0)
+    jumpTable_(new DataEntry<Type>("jumpTable"))
 {}
 
 

From 926b30e70d54501c71490af4734b19fde77bec57 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 29 Oct 2012 12:22:18 +0000
Subject: [PATCH 06/10] ENH: tutorial file update

---
 tutorials/mesh/cvMesh/simpleShapes/Allrun | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tutorials/mesh/cvMesh/simpleShapes/Allrun b/tutorials/mesh/cvMesh/simpleShapes/Allrun
index f1259d32d2..4c2ecc9473 100755
--- a/tutorials/mesh/cvMesh/simpleShapes/Allrun
+++ b/tutorials/mesh/cvMesh/simpleShapes/Allrun
@@ -12,11 +12,11 @@ runApplication surfaceClean \
     constant/triSurface/coneAndSphere_clean.obj
 mv log.surfaceClean log.surfaceClean.coneAndSphere
 # Orient so point to be meshed is inside surface
-surfaceOrient \
+runApplication surfaceOrient \
     constant/triSurface/coneAndSphere_clean.obj \
     -inside '(0 -0.5 0)' \
-    constant/triSurface/coneAndSphere_clean_orient.obj \
-    > log.surfaceOrient.coneAndSphere 2>&1
+    constant/triSurface/coneAndSphere_clean_orient.obj
+mv log.surfaceOrient log.surfaceOrient.coneAndSphere
 
 # Same for outside
 runApplication surfaceClean \
@@ -24,11 +24,11 @@ runApplication surfaceClean \
     1e-4 1e-6 \
     constant/triSurface/domain_clean.stl
 mv log.surfaceClean log.surfaceClean.domain
-surfaceOrient \
+runApplication surfaceOrient \
     constant/triSurface/domain_clean.stl \
     -inside '(0 -0.5 0)' \
     constant/triSurface/domain_clean_orient.stl
-    > log.surfaceOrient.domain 2>&1
+mv log.surfaceOrient log.surfaceOrient.domain
 
 runApplication surfaceFeatureExtract
 mv log.surfaceFeatureExtract log.surfaceFeatureExtract.coneAndSphere_clean

From ab6440105031f516c5499eaad7270edd88a88cc2 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 29 Oct 2012 12:26:39 +0000
Subject: [PATCH 07/10] ENH: put into local Foam namespace

---
 .../jumpCyclic/jumpCyclicFvPatchField.C       | 28 +++++++------------
 .../jumpCyclic/jumpCyclicFvPatchFields.C      | 14 +++-------
 2 files changed, 14 insertions(+), 28 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
index 2287477bdf..6313dab9c1 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
@@ -25,15 +25,10 @@ License
 
 #include "jumpCyclicFvPatchField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
-jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
+Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
@@ -44,7 +39,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 
 
 template<class Type>
-jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
+Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 (
     const jumpCyclicFvPatchField<Type>& ptf,
     const fvPatch& p,
@@ -57,7 +52,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 
 
 template<class Type>
-jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
+Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF,
@@ -72,7 +67,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 
 
 template<class Type>
-jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
+Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 (
     const jumpCyclicFvPatchField<Type>& ptf
 )
@@ -82,7 +77,7 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 
 
 template<class Type>
-jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
+Foam::jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 (
     const jumpCyclicFvPatchField<Type>& ptf,
     const DimensionedField<Type, volMesh>& iF
@@ -95,7 +90,8 @@ jumpCyclicFvPatchField<Type>::jumpCyclicFvPatchField
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-tmp<Field<Type> > jumpCyclicFvPatchField<Type>::patchNeighbourField() const
+Foam::tmp<Foam::Field<Type> >
+Foam::jumpCyclicFvPatchField<Type>::patchNeighbourField() const
 {
     const Field<Type>& iField = this->internalField();
     const labelUList& nbrFaceCells =
@@ -133,7 +129,7 @@ tmp<Field<Type> > jumpCyclicFvPatchField<Type>::patchNeighbourField() const
 
 
 template<class Type>
-void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
+void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 (
     scalarField& result,
     const scalarField& psiInternal,
@@ -144,7 +140,7 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 {
     notImplemented
     (
-        "void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix"
+        "void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix"
         "("
             "scalarField&, "
             "const scalarField&, "
@@ -157,7 +153,7 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 
 
 template<class Type>
-void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
+void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const Field<Type>& psiInternal,
@@ -205,8 +201,4 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
index adfc1326ff..2001e1e2f4 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
@@ -27,19 +27,17 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-makePatchFieldsTypeName(jumpCyclic);
+    makePatchFieldsTypeName(jumpCyclic);
+}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<>
-void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
+void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
 (
     scalarField& result,
     const scalarField& psiInternal,
@@ -88,8 +86,4 @@ void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

From 459aa318fbe4801f807259575295ff11f893ec57 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 29 Oct 2012 12:40:58 +0000
Subject: [PATCH 08/10] ENH: Updated temperature jump BCs

---
 .../energyJump/energyJump/energyJumpFvPatchScalarField.C  | 8 +++++++-
 .../energyJumpAMI/energyJumpAMIFvPatchScalarField.C       | 8 +++++++-
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
index 62ac6516ef..71a097f55f 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
@@ -113,9 +113,15 @@ void Foam::energyJumpFvPatchScalarField::updateCoeffs()
                 thermo.T().boundaryField()[patchID]
             );
 
+        fixedJumpFvPatchScalarField& Tbp =
+            const_cast<fixedJumpFvPatchScalarField&>(TbPatch);
+
+        // force update of jump
+        Tbp.updateCoeffs();
+
         const labelUList& faceCells = this->patch().faceCells();
 
-        jump_ = thermo.he(pp, TbPatch.jump(), faceCells);
+        jump_ = thermo.he(pp, Tbp.jump(), faceCells);
     }
 
     fixedJumpFvPatchField<scalar>::updateCoeffs();
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
index 9f596c3cde..81a992b9ab 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
@@ -113,9 +113,15 @@ void Foam::energyJumpAMIFvPatchScalarField::updateCoeffs()
                 thermo.T().boundaryField()[patchID]
             );
 
+        fixedJumpAMIFvPatchScalarField& Tbp =
+            const_cast<fixedJumpAMIFvPatchScalarField&>(TbPatch);
+
+        // force update of jump
+        Tbp.updateCoeffs();
+
         const labelUList& faceCells = this->patch().faceCells();
 
-        jump_ = thermo.he(pp, TbPatch.jump(), faceCells);
+        jump_ = thermo.he(pp, Tbp.jump(), faceCells);
     }
 
     fixedJumpAMIFvPatchField<scalar>::updateCoeffs();

From 357939a77052e6410c1a14688705ceead6a56486 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 29 Oct 2012 13:20:42 +0000
Subject: [PATCH 09/10] ENH: snappyHexMesh: refine by distance to feature

---
 .../snappyHexMesh/snappyHexMeshDict           |  13 +-
 .../autoHexMeshDriver/autoRefineDriver.C      |   3 +
 .../meshRefinement/meshRefinement.H           |   9 +
 .../meshRefinement/meshRefinementRefine.C     | 119 +++++++++++-
 .../refinementFeatures/refinementFeatures.C   | 180 +++++++++++++++++-
 .../refinementFeatures/refinementFeatures.H   |  34 +++-
 6 files changed, 336 insertions(+), 22 deletions(-)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index f0cc40a8e4..4bef1736f0 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -96,20 +96,29 @@ castellatedMeshControls
     // refinement.
     nCellsBetweenLevels 1;
 
+
     // Explicit feature edge refinement
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Specifies a level for any cell intersected by explicitly provided
     // edges.
     // This is a featureEdgeMesh, read from constant/triSurface for now.
+    // Specify 'levels' in the same way as the 'distance' mode in the
+    // refinementRegions (see below). The old specification
+    //      level   2;
+    // is equivalent to
+    //      levels  ((0 2));
+
     features
     (
         //{
         //    file "someLine.eMesh";
-        //    level 2;
+        //    //level 2;
+        //    levels ((0.0 2) (1.0 3));
         //}
     );
 
+
     // Surface based refinement
     // ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -178,7 +187,7 @@ castellatedMeshControls
     // three modes
     // - distance. 'levels' specifies per distance to the surface the
     //   wanted refinement level. The distances need to be specified in
-    //   descending order.
+    //   increasing order.
     // - inside. 'levels' is only one entry and only the level is used. All
     //   cells inside the surface get refined up to the level. The surface
     //   needs to be closed for this to be possible.
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 985bcfd103..dfd769fd15 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -97,6 +97,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                     refineParams.curvature(),
 
                     true,               // featureRefinement
+                    false,              // featureDistanceRefinement
                     false,              // internalRefinement
                     false,              // surfaceRefinement
                     false,              // curvatureRefinement
@@ -207,6 +208,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
                 refineParams.curvature(),
 
                 false,              // featureRefinement
+                false,              // featureDistanceRefinement
                 false,              // internalRefinement
                 true,               // surfaceRefinement
                 true,               // curvatureRefinement
@@ -368,6 +370,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 refineParams.curvature(),
 
                 false,              // featureRefinement
+                true,               // featureDistanceRefinement
                 true,               // internalRefinement
                 false,              // surfaceRefinement
                 false,              // curvatureRefinement
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 918aae9000..b98230e412 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -264,6 +264,14 @@ private:
                 label& nRefine
             ) const;
 
+            //- Mark cells for distance-to-feature based refinement.
+            label markInternalDistanceToFeatureRefinement
+            (
+                const label nAllowRefine,
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
             //- Mark cells for refinement-shells based refinement.
             label markInternalRefinement
             (
@@ -656,6 +664,7 @@ public:
                 const scalar curvature,
 
                 const bool featureRefinement,
+                const bool featureDistanceRefinement,
                 const bool internalRefinement,
                 const bool surfaceRefinement,
                 const bool curvatureRefinement,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 4d9d10ab80..2c15a847ea 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -290,7 +290,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
         forAll(features_, featI)
         {
             const featureEdgeMesh& featureMesh = features_[featI];
-            const label featureLevel = features_.levels()[featI];
+            const label featureLevel = features_.levels()[featI][0];
             const labelListList& pointEdges = featureMesh.pointEdges();
 
             // Find regions on edgeMesh
@@ -511,6 +511,90 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 }
 
 
+// Mark cells for distance-to-feature based refinement.
+Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
+(
+    const label nAllowRefine,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    // Detect if there are any distance shells
+    if (features_.maxDistance() <= 0.0)
+    {
+        return 0;
+    }
+
+    label oldNRefine = nRefine;
+
+    // Collect cells to test
+    pointField testCc(cellLevel.size()-nRefine);
+    labelList testLevels(cellLevel.size()-nRefine);
+    label testI = 0;
+
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] == -1)
+        {
+            testCc[testI] = cellCentres[cellI];
+            testLevels[testI] = cellLevel[cellI];
+            testI++;
+        }
+    }
+
+    // Do test to see whether cells is inside/outside shell with higher level
+    labelList maxLevel;
+    features_.findHigherLevel(testCc, testLevels, maxLevel);
+
+    // Mark for refinement. Note that we didn't store the original cellID so
+    // now just reloop in same order.
+    testI = 0;
+    forAll(cellLevel, cellI)
+    {
+        if (refineCell[cellI] == -1)
+        {
+            if (maxLevel[testI] > testLevels[testI])
+            {
+                bool reachedLimit = !markForRefine
+                (
+                    maxLevel[testI],    // mark with any positive value
+                    nAllowRefine,
+                    refineCell[cellI],
+                    nRefine
+                );
+
+                if (reachedLimit)
+                {
+                    if (debug)
+                    {
+                        Pout<< "Stopped refining internal cells"
+                            << " since reaching my cell limit of "
+                            << mesh_.nCells()+7*nRefine << endl;
+                    }
+                    break;
+                }
+            }
+            testI++;
+        }
+    }
+
+    if
+    (
+        returnReduce(nRefine, sumOp<label>())
+      > returnReduce(nAllowRefine, sumOp<label>())
+    )
+    {
+        Info<< "Reached refinement limit." << endl;
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
 // Mark cells for non-surface intersection based refinement.
 Foam::label Foam::meshRefinement::markInternalRefinement
 (
@@ -1128,6 +1212,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const scalar curvature,
 
     const bool featureRefinement,
+    const bool featureDistanceRefinement,
     const bool internalRefinement,
     const bool surfaceRefinement,
     const bool curvatureRefinement,
@@ -1196,8 +1281,24 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 nRefine
             );
 
-            Info<< "Marked for refinement due to explicit features    : "
-                << nFeatures << " cells."  << endl;
+            Info<< "Marked for refinement due to explicit features             "
+                << ": " << nFeatures << " cells."  << endl;
+        }
+
+        // Inside distance-to-feature shells
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if (featureDistanceRefinement)
+        {
+            label nShell = markInternalDistanceToFeatureRefinement
+            (
+                nAllowRefine,
+
+                refineCell,
+                nRefine
+            );
+            Info<< "Marked for refinement due to distance to explicit features "
+                ": " << nShell << " cells."  << endl;
         }
 
         // Inside refinement shells
@@ -1212,8 +1313,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 refineCell,
                 nRefine
             );
-            Info<< "Marked for refinement due to refinement shells    : "
-                << nShell << " cells."  << endl;
+            Info<< "Marked for refinement due to refinement shells             "
+                << ": " << nShell << " cells."  << endl;
         }
 
         // Refinement based on intersection of surface
@@ -1230,8 +1331,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 refineCell,
                 nRefine
             );
-            Info<< "Marked for refinement due to surface intersection : "
-                << nSurf << " cells."  << endl;
+            Info<< "Marked for refinement due to surface intersection          "
+                << ": " << nSurf << " cells."  << endl;
         }
 
         // Refinement based on curvature of surface
@@ -1254,8 +1355,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
                 refineCell,
                 nRefine
             );
-            Info<< "Marked for refinement due to curvature/regions    : "
-                << nCurv << " cells."  << endl;
+            Info<< "Marked for refinement due to curvature/regions             "
+                << ": " << nCurv << " cells."  << endl;
         }
 
         // Pack cells-to-refine
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
index 0cdd440d1a..e72f13d726 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -25,6 +25,7 @@ License
 
 #include "refinementFeatures.H"
 #include "Time.H"
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -34,15 +35,15 @@ void Foam::refinementFeatures::read
     const PtrList<dictionary>& featDicts
 )
 {
-    forAll(featDicts, i)
+    forAll(featDicts, featI)
     {
-        const dictionary& dict = featDicts[i];
+        const dictionary& dict = featDicts[featI];
 
         fileName featFileName(dict.lookup("file"));
 
         set
         (
-            i,
+            featI,
             new featureEdgeMesh
             (
                 IOobject
@@ -58,15 +59,74 @@ void Foam::refinementFeatures::read
             )
         );
 
-        const featureEdgeMesh& eMesh = operator[](i);
+        const featureEdgeMesh& eMesh = operator[](featI);
 
         //eMesh.mergePoints(meshRefiner_.mergeDistance());
-        levels_[i] = readLabel(dict.lookup("level"));
 
-        Info<< "Refinement level " << levels_[i]
-            << " for all cells crossed by feature " << featFileName
-            << " (" << eMesh.points().size() << " points, "
+        if (dict.found("levels"))
+        {
+            List<Tuple2<scalar, label> > distLevels(dict["levels"]);
+
+            if (dict.size() < 1)
+            {
+                FatalErrorIn
+                (
+                    "refinementFeatures::read"
+                    "(const objectRegistry&"
+                    ", const PtrList<dictionary>&)"
+                )   << " : levels should be at least size 1" << endl
+                    << "levels : "  << dict["levels"]
+                    << exit(FatalError);
+            }
+
+            distances_[featI].setSize(distLevels.size());
+            levels_[featI].setSize(distLevels.size());
+
+            forAll(distLevels, j)
+            {
+                distances_[featI][j] = distLevels[j].first();
+                levels_[featI][j] = distLevels[j].second();
+
+                // Check in incremental order
+                if (j > 0)
+                {
+                    if
+                    (
+                        (distances_[featI][j] <= distances_[featI][j-1])
+                     || (levels_[featI][j] > levels_[featI][j-1])
+                    )
+                    {
+                        FatalErrorIn
+                        (
+                            "refinementFeatures::read"
+                            "(const objectRegistry&"
+                            ", const PtrList<dictionary>&)"
+                        )   << " : Refinement should be specified in order"
+                            << " of increasing distance"
+                            << " (and decreasing refinement level)." << endl
+                            << "Distance:" << distances_[featI][j]
+                            << " refinementLevel:" << levels_[featI][j]
+                            << exit(FatalError);
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Look up 'level' for single level
+            levels_[featI] = labelList(1, readLabel(dict.lookup("level")));
+            distances_[featI] = scalarField(1, 0.0);
+        }
+
+        Info<< "Refinement level according to distance to "
+            << featFileName << " (" << eMesh.points().size() << " points, "
             << eMesh.edges().size() << " edges)." << endl;
+        forAll(levels_[featI], j)
+        {
+            Info<< "    level " << levels_[featI][j]
+                << " for all cells within " << distances_[featI][j]
+                << " meter." << endl;
+        }
     }
 }
 
@@ -127,6 +187,80 @@ void Foam::refinementFeatures::buildTrees
 }
 
 
+
+// Find maximum level of a shell.
+void Foam::refinementFeatures::findHigherLevel
+(
+    const pointField& pt,
+    const label featI,
+    labelList& maxLevel
+) const
+{
+    const labelList& levels = levels_[featI];
+
+    const scalarField& distances = distances_[featI];
+
+    // Collect all those points that have a current maxLevel less than
+    // (any of) the shell. Also collect the furthest distance allowable
+    // to any shell with a higher level.
+
+    pointField candidates(pt.size());
+    labelList candidateMap(pt.size());
+    scalarField candidateDistSqr(pt.size());
+    label candidateI = 0;
+
+    forAll(maxLevel, pointI)
+    {
+        forAllReverse(levels, levelI)
+        {
+            if (levels[levelI] > maxLevel[pointI])
+            {
+                candidates[candidateI] = pt[pointI];
+                candidateMap[candidateI] = pointI;
+                candidateDistSqr[candidateI] = sqr(distances[levelI]);
+                candidateI++;
+                break;
+            }
+        }
+    }
+    candidates.setSize(candidateI);
+    candidateMap.setSize(candidateI);
+    candidateDistSqr.setSize(candidateI);
+
+    // Do the expensive nearest test only for the candidate points.
+    const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
+
+    List<pointIndexHit> nearInfo(candidates.size());
+    forAll(candidates, candidateI)
+    {
+        nearInfo[candidateI] = tree.findNearest
+        (
+            candidates[candidateI],
+            candidateDistSqr[candidateI]
+        );
+    }
+
+    // Update maxLevel
+    forAll(nearInfo, candidateI)
+    {
+        if (nearInfo[candidateI].hit())
+        {
+            // Check which level it actually is in.
+            label minDistI = findLower
+            (
+                distances,
+                mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
+            );
+
+            label pointI = candidateMap[candidateI];
+
+            // pt is inbetween shell[minDistI] and shell[minDistI+1]
+            maxLevel[pointI] = levels[minDistI+1];
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::refinementFeatures::refinementFeatures
@@ -136,6 +270,7 @@ Foam::refinementFeatures::refinementFeatures
 )
 :
     PtrList<featureEdgeMesh>(featDicts.size()),
+    distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
     pointTrees_(featDicts.size())
@@ -175,6 +310,7 @@ Foam::refinementFeatures::refinementFeatures
 )
 :
     PtrList<featureEdgeMesh>(featDicts.size()),
+    distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
     pointTrees_(featDicts.size())
@@ -336,4 +472,32 @@ void Foam::refinementFeatures::findNearestPoint
 }
 
 
+void Foam::refinementFeatures::findHigherLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    labelList& maxLevel
+) const
+{
+    // Maximum level of any shell. Start off with level of point.
+    maxLevel = ptLevel;
+
+    forAll(*this, featI)
+    {
+        findHigherLevel(pt, featI, maxLevel);
+    }
+}
+
+
+Foam::scalar Foam::refinementFeatures::maxDistance() const
+{
+    scalar overallMax = -GREAT;
+    forAll(distances_, featI)
+    {
+        overallMax = max(overallMax, max(distances_[featI]));
+    }
+    return overallMax;
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
index 5442e140dd..aa2807bac1 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
@@ -57,8 +57,11 @@ private:
 
     // Private data
 
-        //- Refinement levels
-        labelList levels_;
+        //- Per shell the list of ranges
+        List<scalarField> distances_;
+
+        //- Per shell per distance the refinement level
+        labelListList levels_;
 
         //- Edge
         PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
@@ -75,6 +78,13 @@ private:
         //- Build edge tree and feature point tree
         void buildTrees(const label, const labelList&);
 
+        //- Find shell level higher than ptLevel
+        void findHigherLevel
+        (
+            const pointField& pt,
+            const label featI,
+            labelList& maxLevel
+        ) const;
 
 public:
 
@@ -101,11 +111,18 @@ public:
 
         // Access
 
-            const labelList& levels() const
+            //- Per featureEdgeMesh the list of level
+            const labelListList& levels() const
             {
                 return levels_;
             }
 
+            //- Per featureEdgeMesh the list of ranges
+            const List<scalarField>& distances() const
+            {
+                return distances_;
+            }
+
             const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
             {
                 return edgeTrees_;
@@ -119,6 +136,9 @@ public:
 
         // Query
 
+            //- Highest distance of all features
+            scalar maxDistance() const;
+
             //- Find nearest point on nearest feature edge
             void findNearestEdge
             (
@@ -141,6 +161,14 @@ public:
                 labelList& nearIndex
             ) const;
 
+            //- Find shell level higher than ptLevel
+            void findHigherLevel
+            (
+                const pointField& pt,
+                const labelList& ptLevel,
+                labelList& maxLevel
+            ) const;
+
 };
 
 

From 49130c846782982093fedebef8d6e6ab8654a6e5 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 29 Oct 2012 14:31:23 +0000
Subject: [PATCH 10/10] ENH: Updated min/max function object to include
 boundary values

---
 .../field/fieldMinMax/fieldMinMaxTemplates.C  | 95 ++++++++++++++-----
 1 file changed, 71 insertions(+), 24 deletions(-)

diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
index 64a03c43b1..33564024b1 100644
--- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
@@ -45,31 +45,56 @@ void Foam::fieldMinMax::calcMinMaxFields
         const label procI = Pstream::myProcNo();
 
         const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
+        const fvMesh& mesh = field.mesh();
+
+        const volVectorField::GeometricBoundaryField& CfBoundary =
+            mesh.C().boundaryField();
+
         switch (mode)
         {
             case mdMag:
             {
-                const scalarField magField(mag(field));
+                const volScalarField magField(mag(field));
+                const volScalarField::GeometricBoundaryField& magFieldBoundary =
+                    magField.boundaryField();
 
-                labelList minIs(Pstream::nProcs());
                 scalarList minVs(Pstream::nProcs());
                 List<vector> minCs(Pstream::nProcs());
-                minIs[procI] = findMin(magField);
-                minVs[procI] = magField[minIs[procI]];
-                minCs[procI] = field.mesh().C()[minIs[procI]];
+                label minProcI = findMin(magField);
+                minVs[procI] = magField[minProcI];
+                minCs[procI] = field.mesh().C()[minProcI];
 
-                Pstream::gatherList(minIs);
-                Pstream::gatherList(minVs);
-                Pstream::gatherList(minCs);
 
                 labelList maxIs(Pstream::nProcs());
                 scalarList maxVs(Pstream::nProcs());
                 List<vector> maxCs(Pstream::nProcs());
-                maxIs[procI] = findMax(magField);
-                maxVs[procI] = magField[maxIs[procI]];
-                maxCs[procI] = field.mesh().C()[maxIs[procI]];
+                label maxProcI = findMax(magField);
+                maxVs[procI] = magField[maxProcI];
+                maxCs[procI] = field.mesh().C()[maxProcI];
+
+                forAll(magFieldBoundary, patchI)
+                {
+                    const scalarField& mfp = magFieldBoundary[patchI];
+                    const vectorField& Cfp = CfBoundary[patchI];
+
+                    label minPI = findMin(mfp);
+                    if (mfp[minPI] < minVs[procI])
+                    {
+                        minVs[procI] = mfp[minPI];
+                        minCs[procI] = Cfp[minPI];
+                    }
+
+                    label maxPI = findMax(mfp);
+                    if (mfp[maxPI] > maxVs[procI])
+                    {
+                        maxVs[procI] = mfp[maxPI];
+                        maxCs[procI] = Cfp[maxPI];
+                    }
+                }
+
+                Pstream::gatherList(minVs);
+                Pstream::gatherList(minCs);
 
-                Pstream::gatherList(maxIs);
                 Pstream::gatherList(maxVs);
                 Pstream::gatherList(maxCs);
 
@@ -127,25 +152,47 @@ void Foam::fieldMinMax::calcMinMaxFields
             }
             case mdCmpt:
             {
-                List<Type> minVs(Pstream::nProcs());
-                labelList minIs(Pstream::nProcs());
-                List<vector> minCs(Pstream::nProcs());
-                minIs[procI] = findMin(field);
-                minVs[procI] = field[minIs[procI]];
-                minCs[procI] = field.mesh().C()[minIs[procI]];
+                const typename fieldType::GeometricBoundaryField&
+                    fieldBoundary = field.boundaryField();
+
+                List<Type> minVs(Pstream::nProcs());
+                List<vector> minCs(Pstream::nProcs());
+                label minProcI = findMin(field);
+                minVs[procI] = field[minProcI];
+                minCs[procI] = field.mesh().C()[minProcI];
 
-                Pstream::gatherList(minIs);
                 Pstream::gatherList(minVs);
                 Pstream::gatherList(minCs);
 
                 List<Type> maxVs(Pstream::nProcs());
-                labelList maxIs(Pstream::nProcs());
                 List<vector> maxCs(Pstream::nProcs());
-                maxIs[procI] = findMax(field);
-                maxVs[procI] = field[maxIs[procI]];
-                maxCs[procI] = field.mesh().C()[maxIs[procI]];
+                label maxProcI = findMax(field);
+                maxVs[procI] = field[maxProcI];
+                maxCs[procI] = field.mesh().C()[maxProcI];
+
+                forAll(fieldBoundary, patchI)
+                {
+                    const Field<Type>& fp = fieldBoundary[patchI];
+                    const vectorField& Cfp = CfBoundary[patchI];
+
+                    label minPI = findMin(fp);
+                    if (fp[minPI] < minVs[procI])
+                    {
+                        minVs[procI] = fp[minPI];
+                        minCs[procI] = Cfp[minPI];
+                    }
+
+                    label maxPI = findMax(fp);
+                    if (fp[maxPI] > maxVs[procI])
+                    {
+                        maxVs[procI] = fp[maxPI];
+                        maxCs[procI] = Cfp[maxPI];
+                    }
+                }
+
+                Pstream::gatherList(minVs);
+                Pstream::gatherList(minCs);
 
-                Pstream::gatherList(maxIs);
                 Pstream::gatherList(maxVs);
                 Pstream::gatherList(maxCs);