Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -176,6 +176,10 @@ public:
         //- Find patch index given a name
         label findPatchID(const word& patchName) const;
 
+        //- Find patch indices for a given polyPatch type
+        template<class Type>
+        labelHashSet findPatchIDs() const;
+
         //- Return patch index for a given face label
         label whichPatch(const label faceIndex) const;
 
@@ -235,6 +239,7 @@ public:
             IOstream::compressionType cmp
         ) const;
 
+
     // Member Operators
 
         //- Return const and non-const reference to polyPatch by index.
@@ -259,6 +264,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "polyBoundaryMeshTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMeshTemplates.C

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::labelHashSet Foam::polyBoundaryMesh::findPatchIDs() const
+{
+    const polyBoundaryMesh& bm = *this;
+
+    labelHashSet patchIDs(bm.size());
+
+    forAll(bm, patchI)
+    {
+        if (isA<Type>(bm[patchI]))
+        {
+            patchIDs.insert(patchI);
+        }
+    }
+    return patchIDs;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/Make/files

@@ -36,6 +36,7 @@ $(derivedFvPatches)/regionCoupled/regionCoupledWallFvPatch.C
 
 
 wallDist = fvMesh/wallDist
+$(wallDist)/patchDist.C
 $(wallDist)/wallPointYPlus/wallPointYPlus.C
 $(wallDist)/nearWallDistNoSearch.C
 $(wallDist)/nearWallDist.C

src/finiteVolume/fvMesh/wallDist/patchDist.C

@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchDist.H"
+#include "patchWave.H"
+#include "fvMesh.H"
+#include "emptyFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchDist::patchDist
+(
+    const fvMesh& mesh,
+    const labelHashSet& patchIDs,
+    const bool correctWalls
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            "y",
+            mesh.time().timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("y", dimLength, GREAT)
+    ),
+    patchIDs_(patchIDs),
+    correctWalls_(correctWalls),
+    nUnset_(0)
+{
+    patchDist::correct();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchDist::~patchDist()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::patchDist::correct()
+{
+    // Calculate distance starting from patch faces
+    patchWave wave(mesh(), patchIDs_, correctWalls_);
+
+    // Transfer cell values from wave into *this
+    transfer(wave.distance());
+
+    // Transfer values on patches into boundaryField of *this
+    forAll(boundaryField(), patchI)
+    {
+        if (!isA<emptyFvPatchScalarField>(boundaryField()[patchI]))
+        {
+            scalarField& waveFld = wave.patchDistance()[patchI];
+
+            boundaryField()[patchI].transfer(waveFld);
+        }
+    }
+
+    // Transfer number of unset values
+    nUnset_ = wave.nUnset();
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/wallDist/patchDist.H

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::patchDist
+
+Description
+    Calculation of distance to nearest patch for all cells and boundary.
+    Uses meshWave to do actual calculation.
+
+    Distance correction:
+
+    if correctWalls = true:
+    For each cell with face on wall calculate the true nearest point
+    (by triangle decomposition) on that face and do the same for that face's
+    pointNeighbours. This will find the true nearest distance in almost all
+    cases. Only very skewed cells or cells close to another wall might be
+    missed.
+
+    For each cell with only one point on wall the same is done except now it
+    takes the pointFaces() of the wall point to look for the nearest point.
+
+Note
+
+    correct() : for now does complete recalculation. (which usually is
+    ok since mesh is smoothed). However for topology change where geometry
+    in most of domain does not change you could think of starting from the
+    old cell values. Tried but not done since:
+    - meshWave would have to be called with old cellInfo.
+      This is List\<wallInfo\> of nCells.
+    - cannot construct from distance (y_) only since we don't know a value
+      for origin_. (origin_ = GREAT already used to denote illegal value.)
+    - so we would have to store a List\<wallInfo\> which unfortunately does
+      not get resized/mapped automatically upon mesh changes.
+
+SourceFiles
+    patchDist.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchDist_H
+#define patchDist_H
+
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                          Class patchDist Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchDist
+:
+    public volScalarField
+{
+
+
+private:
+
+    // Private Member Data
+
+        //- Set of patch IDs
+        labelHashSet patchIDs_;
+
+        //- Do accurate distance calculation for near-wall cells.
+        bool correctWalls_;
+
+        //- Number of unset cells and faces.
+        label nUnset_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        patchDist(const patchDist&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const patchDist&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh and flag whether or not to correct wall.
+        //  Calculate for all cells. correctWalls : correct wall (face&point)
+        //  cells for correct distance, searching neighbours.
+        patchDist
+        (
+            const fvMesh& mesh,
+            const labelHashSet& patchIDs,
+            const bool correctWalls = true
+        );
+
+
+    //- Destructor
+    virtual ~patchDist();
+
+
+    // Member Functions
+
+        const volScalarField& y() const
+        {
+            return *this;
+        }
+
+        label nUnset() const
+        {
+            return nUnset_;
+        }
+
+        //- Correct for mesh geom/topo changes
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/wallDist/wallDist.C

@@ -24,34 +24,24 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "wallDist.H"
-#include "patchWave.H"
 #include "fvMesh.H"
 #include "wallPolyPatch.H"
-#include "fvPatchField.H"
-#include "Field.H"
-#include "emptyFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::wallDist::wallDist(const fvMesh& mesh, const bool correctWalls)
+Foam::wallDist::wallDist
+(
+    const fvMesh& mesh,
+    const bool correctWalls
+)
 :
-    volScalarField
+    patchDist
     (
-        IOobject
-        (
-            "y",
-            mesh.time().timeName(),
-            mesh
-        ),
         mesh,
-        dimensionedScalar("y", dimLength, GREAT)
+        mesh.boundaryMesh().findPatchIDs<wallPolyPatch>(),
-    ),
+        correctWalls
-    cellDistFuncs(mesh),
+    )
-    correctWalls_(correctWalls),
+{}
-    nUnset_(0)
-{
-    wallDist::correct();
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -60,33 +50,4 @@ Foam::wallDist::~wallDist()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::wallDist::correct()
-{
-    // Get patchids of walls
-    labelHashSet wallPatchIDs(getPatchIDs<wallPolyPatch>());
-
-    // Calculate distance starting from wallPatch faces.
-    patchWave wave(cellDistFuncs::mesh(), wallPatchIDs, correctWalls_);
-
-    // Transfer cell values from wave into *this
-    transfer(wave.distance());
-
-    // Transfer values on patches into boundaryField of *this
-    forAll(boundaryField(), patchI)
-    {
-        if (!isA<emptyFvPatchScalarField>(boundaryField()[patchI]))
-        {
-            scalarField& waveFld = wave.patchDistance()[patchI];
-
-            boundaryField()[patchI].transfer(waveFld);
-        }
-    }
-
-    // Transfer number of unset values
-    nUnset_ = wave.nUnset();
-}
-
-
 // ************************************************************************* //

src/finiteVolume/fvMesh/wallDist/wallDist.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,33 +25,7 @@ Class
     Foam::wallDist
 
 Description
-    Calculation of distance to nearest wall for all cells and boundary.
+    Specialisation of patchDist for wall distance calculation
-    Uses meshWave to do actual calculation.
-
-    Distance correction:
-
-    if correctWalls = true:
-    For each cell with face on wall calculate the true nearest point
-    (by triangle decomposition) on that face and do that same for that face's
-    pointNeighbours. This will find the true nearest distance in almost all
-    cases. Only very skewed cells or cells close to another wall might be
-    missed.
-
-    For each cell with only point on wall the same is done except now it takes
-    the pointFaces() of the wall point to look for the nearest point.
-
-Note
-
-    correct() : for now does complete recalculation. (which usually is
-    ok since mesh is smoothed). However for topology change where geometry
-    in most of domain does not change you could think of starting from the
-    old cell values. Tried but not done since:
-    - meshWave would have to be called with old cellInfo.
-      This is List\<wallInfo\> of nCells.
-    - cannot construct from distance (y_) only since we don't know a value
-      for origin_. (origin_ = GREAT already used to denote illegal value.)
-    - so we would have to store a List\<wallInfo\> which unfortunately does
-      not get resized/mapped automatically upon mesh changes.
 
 SourceFiles
     wallDist.C
@@ -61,9 +35,7 @@ SourceFiles
 #ifndef wallDist_H
 #define wallDist_H
 
-#include "volFields.H"
+#include "patchDist.H"
-#include "cellDistFuncs.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,27 +45,16 @@ namespace Foam
 class fvMesh;
 
 /*---------------------------------------------------------------------------*\
-                           Class wallDist Declaration
+                          Class wallDist Declaration
 \*---------------------------------------------------------------------------*/
 
 class wallDist
 :
-    public volScalarField,
+    public patchDist
-    public cellDistFuncs
 {
 
-
 private:
 
-    // Private Member Data
-
-        //- Do accurate distance calculation for near-wall cells.
-        bool correctWalls_;
-
-        //- Number of unset cells and faces.
-        label nUnset_;
-
-
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -110,27 +71,15 @@ public:
         //- Construct from mesh and flag whether or not to correct wall.
         //  Calculate for all cells. correctWalls : correct wall (face&point)
         //  cells for correct distance, searching neighbours.
-        wallDist(const fvMesh& mesh, bool correctWalls = true);
+        wallDist
+        (
+            const fvMesh& mesh,
+            const bool correctWalls = true
+        );
 
 
     //- Destructor
     virtual ~wallDist();
-
-
-    // Member Functions
-
-        const volScalarField& y() const
-        {
-            return *this;
-        }
-
-        label nUnset() const
-        {
-            return nUnset_;
-        }
-
-        //- Correct for mesh geom/topo changes
-        virtual void correct();
 };
 
 

src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,7 @@ License
 #include "fvcGrad.H"
 #include "unitConversion.H"
 #include "fvPatchField.H"
+#include "patchDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,6 +44,35 @@ namespace surfaceFilmModels
 defineTypeNameAndDebug(contactAngleForce, 0);
 addToRunTimeSelectionTable(force, contactAngleForce, dictionary);
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void contactAngleForce::initialise()
+{
+    const wordReList zeroForcePatches(coeffs_.lookup("zeroForcePatches"));
+
+    if (zeroForcePatches.size())
+    {
+        const polyBoundaryMesh& pbm = owner_.regionMesh().boundaryMesh();
+        scalar dLim = readScalar(coeffs_.lookup("zeroForceDistance"));
+
+        Info<< "        Assigning zero contact force within " << dLim
+            << " of patches:" << endl;
+
+        labelHashSet patchIDs = pbm.patchSet(zeroForcePatches);
+
+        forAllConstIter(labelHashSet, patchIDs, iter)
+        {
+            label patchI = iter.key();
+            Info<< "            " << pbm[patchI].name() << endl;
+        }
+
+        patchDist dist(owner_.regionMesh(), patchIDs);
+
+        mask_ = pos(dist - dimensionedScalar("dLim", dimLength, dLim));
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 contactAngleForce::contactAngleForce
@@ -61,8 +91,23 @@ contactAngleForce::contactAngleForce
             coeffs_.subDict("contactAngleDistribution"),
             rndGen_
         )
+    ),
+    mask_
+    (
+        IOobject
+        (
+            typeName + ".contactForceMask",
+            owner_.time().timeName(),
+            owner_.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        owner_.regionMesh(),
+        dimensionedScalar("mask", dimless, 0.0)
     )
-{}
+{
+    initialise();
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -121,7 +166,7 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
             cellI = cellN;
         }
 
-        if (cellI != -1)
+        if (cellI != -1 && mask_[cellI] > 0)
         {
             const scalar dx = owner_.regionMesh().deltaCoeffs()[faceI];
             const vector n =
@@ -137,20 +182,26 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         if (!owner().isCoupledPatch(patchI))
         {
             const fvPatchField<scalar>& alphaf = alpha.boundaryField()[patchI];
+            const fvPatchField<scalar>& maskf = mask_.boundaryField()[patchI];
             const scalarField& dx = alphaf.patch().deltaCoeffs();
             const labelUList& faceCells = alphaf.patch().faceCells();
 
             forAll(alphaf, faceI)
             {
-                label cellO = faceCells[faceI];
+                if (maskf[faceI] > 0)
-
-                if ((alpha[cellO] > 0.5) && (alphaf[faceI] < 0.5))
                 {
-                    const vector n =
+                    label cellO = faceCells[faceI];
-                        gradAlpha[cellO]/(mag(gradAlpha[cellO]) + ROOTVSMALL);
+
-                    scalar theta = cos(degToRad(distribution_->sample()));
+                    if ((alpha[cellO] > 0.5) && (alphaf[faceI] < 0.5))
-                    force[cellO] += Ccf_*n*sigma[cellO]*(1.0 - theta)/dx[faceI];
+                    {
-                    nHits[cellO]++;
+                        const vector n =
+                            gradAlpha[cellO]
+                           /(mag(gradAlpha[cellO]) + ROOTVSMALL);
+                        scalar theta = cos(degToRad(distribution_->sample()));
+                        force[cellO] +=
+                            Ccf_*n*sigma[cellO]*(1.0 - theta)/dx[faceI];
+                        nHits[cellO]++;
+                    }
                 }
             }
         }

src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,9 @@ Class
 Description
     Film contact angle force
 
+    The effect of the contact angle force can be ignored over a specified
+    distance from patches.
+
 SourceFiles
     contactAngleForce.C
 
@@ -69,10 +72,15 @@ private:
         //- Parcel size PDF model
         const autoPtr<distributionModels::distributionModel> distribution_;
 
+        //- Mask over which force is applied
+        volScalarField mask_;
 
 
     // Private member functions
 
+        //- Initialise
+        void initialise();
+
         //- Disallow default bitwise copy construct
         contactAngleForce(const contactAngleForce&);