merge into master

applications/solvers/Lagrangian/reactingParcelFoam/Make/files

@@ -0,0 +1,3 @@
+reactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/reactingParcelFoam

applications/solvers/combustion/reactingFoam/createFields.H

@@ -37,7 +37,6 @@ volVectorField U
 
 volScalarField& p = thermo->p();
 const volScalarField& psi = thermo->psi();
-const volScalarField& T = thermo->T();
 volScalarField& h = thermo->h();
 
 
@@ -99,5 +98,5 @@ volScalarField dQ
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
+    dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0.0)
 );

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files

@@ -1,10 +1,6 @@
 regionProperties/regionProperties.C
 
-coupleManager/coupleManager.C
-
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
-derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
-derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
 derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 fluid/compressibleCourantNo.C

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -2,9 +2,7 @@ EXE_INC = \
     -Ifluid \
     -Isolid \
     -IregionProperties \
-    -IcoupleManager \
+    -I$(LIB_SRC)/meshTools/lnInclude \
-    -IderivedFvPatchFields/solidWallTemperatureCoupled \
-    -IderivedFvPatchFields/solidWallHeatFluxTemperatureCoupled \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -67,6 +67,7 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
+#       include "readPIMPLEControls.H"
 
         if (fluidRegions.size())
         {
@@ -78,22 +79,36 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        forAll(fluidRegions, i)
+        if (nOuterCorr != 1)
         {
-            Info<< "\nSolving for fluid region "
+            forAll(fluidRegions, i)
-                << fluidRegions[i].name() << endl;
+            {
-#           include "setRegionFluidFields.H"
+#               include "setRegionFluidFields.H"
-#           include "readFluidMultiRegionPISOControls.H"
+#               include "storeOldFluidFields.H"
-#           include "solveFluid.H"
+            }
         }
 
-        forAll(solidRegions, i)
+
+        // --- PIMPLE loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
-            Info<< "\nSolving for solid region "
+            forAll(fluidRegions, i)
-                << solidRegions[i].name() << endl;
+            {
-#           include "setRegionSolidFields.H"
+                Info<< "\nSolving for fluid region "
-#           include "readSolidMultiRegionPISOControls.H"
+                    << fluidRegions[i].name() << endl;
-#           include "solveSolid.H"
+#               include "setRegionFluidFields.H"
+#               include "readFluidMultiRegionPIMPLEControls.H"
+#               include "solveFluid.H"
+            }
+
+            forAll(solidRegions, i)
+            {
+                Info<< "\nSolving for solid region "
+                    << solidRegions[i].name() << endl;
+#               include "setRegionSolidFields.H"
+#               include "readSolidMultiRegionPIMPLEControls.H"
+#               include "solveSolid.H"
+            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C

@@ -1,186 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "coupleManager.H"
-#include "OFstream.H"
-#include "regionProperties.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::coupleManager::coupleManager(const fvPatch& patch)
-:
-    patch_(patch),
-    neighbourRegionName_("undefined-neighbourRegionName"),
-    neighbourPatchName_("undefined-neighbourPatchName"),
-    neighbourFieldName_("undefined-neighbourFieldName"),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-Foam::coupleManager::coupleManager
-(
-    const fvPatch& patch,
-    const dictionary& dict
-)
-:
-    patch_(patch),
-    neighbourRegionName_(dict.lookup("neighbourRegionName")),
-    neighbourPatchName_(dict.lookup("neighbourPatchName")),
-    neighbourFieldName_(dict.lookup("neighbourFieldName")),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-Foam::coupleManager::coupleManager
-(
-    const coupleManager& cm
-)
-:
-    patch_(cm.patch()),
-    neighbourRegionName_(cm.neighbourRegionName()),
-    neighbourPatchName_(cm.neighbourPatchName()),
-    neighbourFieldName_(cm.neighbourFieldName()),
-    localRegion_(patch_.boundaryMesh().mesh())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::coupleManager::~coupleManager()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::coupleManager::regionOwner() const
-{
-    const fvMesh& nbrRegion = neighbourRegion();
-
-    const regionProperties& props =
-        localRegion_.objectRegistry::parent().lookupObject<regionProperties>
-        (
-            "regionProperties"
-        );
-
-    label myIndex = findIndex(props.fluidRegionNames(), localRegion_.name());
-    if (myIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), localRegion_.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::regionOwner() const")
-                << "Cannot find region " << localRegion_.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        myIndex = props.fluidRegionNames().size() + i;
-    }
-    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
-    if (nbrIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::regionOwner() const")
-                << "Cannot find region " << nbrRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        nbrIndex = props.fluidRegionNames().size() + i;
-    }
-
-    return myIndex < nbrIndex;
-}
-
-
-void Foam::coupleManager::checkCouple() const
-{
-    Info<< "neighbourRegionName_ = " << neighbourRegionName_ << endl;
-    Info<< "neighbourPatchName_ = " << neighbourPatchName_ << endl;
-    Info<< "neighbourFieldName_ = " << neighbourFieldName_ << endl;
-
-    const fvPatch& nPatch = neighbourPatch();
-
-    if (patch_.size() != nPatch.size())
-    {
-        FatalErrorIn("Foam::coupleManager::checkCouple()")
-            << "Unequal patch sizes:" << nl
-            << "    patch name (size) = " << patch_.name()
-            << "(" << patch_.size() << ")" << nl
-            << "    neighbour patch name (size) = "
-            << nPatch.name() << "(" << patch_.size() << ")" << nl
-            << abort(FatalError);
-    }
-}
-
-
-void Foam::coupleManager::coupleToObj() const
-{
-    const fvPatch& nPatch = neighbourPatch();
-
-    OFstream obj
-        (
-             patch_.name() + "_to_" + nPatch.name() + "_couple.obj"
-        );
-    const vectorField& c1 = patch_.Cf();
-    const vectorField& c2 = neighbourPatch().Cf();
-
-    if (c1.size() != c2.size())
-    {
-        FatalErrorIn("coupleManager::coupleToObj() const")
-            << "Coupled patches are of unequal size:" << nl
-            << "    patch0 = " << patch_.name()
-            << "(" << patch_.size() <<  ")" << nl
-            << "    patch1 = " << nPatch.name()
-            << "(" << nPatch.size() <<  ")" << nl
-            << abort(FatalError);
-    }
-
-    forAll(c1, i)
-    {
-        obj << "v " << c1[i].x() << " " << c1[i].y() << " " << c1[i].z() << nl
-            << "v " << c2[i].x() << " " << c2[i].y() << " " << c2[i].z() << nl
-            << "l " << (2*i + 1) << " " << (2*i + 2) << endl;
-    }
-}
-
-
-void Foam::coupleManager::writeEntries(Ostream& os) const
-{
-    os.writeKeyword("neighbourRegionName");
-    os << neighbourRegionName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("neighbourPatchName");
-    os << neighbourPatchName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("neighbourFieldName");
-    os << neighbourFieldName_ << token::END_STATEMENT << nl;
-}
-
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H

@@ -1,170 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    coupleManager
-
-Description
-    Object to handle the coupling of region patches. It can be queried to
-    return the neighbour information.
-
-SourceFiles
-    coupleManager.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef coupleManager_H
-#define coupleManager_H
-
-#include "Ostream.H"
-#include "dictionary.H"
-#include "fvPatch.H"
-#include "fvMesh.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class coupleManager Declaration
-\*---------------------------------------------------------------------------*/
-
-class coupleManager
-{
-    // Private data
-
-        //- Reference to the local fvPatch
-        const fvPatch& patch_;
-
-        //- Name of neighbour region
-        word neighbourRegionName_;
-
-        //- Name of patch on the neighbour region
-        word neighbourPatchName_;
-
-        //- Name of field on the neighbour region
-        word neighbourFieldName_;
-
-        //- Reference to the local region
-        const fvMesh& localRegion_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise assignment
-        void operator=(const coupleManager&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from fvPatch
-        coupleManager(const fvPatch& patch);
-
-        //- Construct from fvPatch and dictionary
-        coupleManager(const fvPatch& patch, const dictionary& dict);
-
-        //- Copy constructor
-        coupleManager(const coupleManager& cm);
-
-
-    // Destructor
-
-        ~coupleManager();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return a reference to the local patch
-            inline const fvPatch& patch() const;
-
-            //- Return the name of the neighbour region
-            inline const word& neighbourRegionName() const;
-
-            //- Return the name of the patch on the neighbour region
-            inline const word& neighbourPatchName() const;
-
-            //- Return the name of the field on the neighbour region
-            inline const word& neighbourFieldName() const;
-
-            //- Return a reference to the neighbour mesh
-            inline const fvMesh& neighbourRegion() const;
-
-            //- Return the neighbour patch ID
-            inline label neighbourPatchID() const;
-
-            //- Return a reference to the neighbour patch
-            inline const fvPatch& neighbourPatch() const;
-
-            //- Return a reference to the neighbour patch field
-            template<class Type>
-            inline const fvPatchField<Type>& neighbourPatchField() const;
-
-            //- Am I owner (= first to evaluate) of this region interface? 
-            bool regionOwner() const;
-
-            //- Check that the couple is valid
-            void checkCouple() const;
-
-
-        // Edit
-
-            //- Return the name of the neighbour region
-            word& neighbourRegionName();
-
-            //- Return the name of the patch on the neighbour region
-            word& neighbourPatchName();
-
-            //- Return the name of the field on the neighbour region
-            word& neighbourFieldName();
-
-
-        // Write
-
-           //- Write couple to obj file
-           void coupleToObj() const;
-
-           //- Write entries for patches
-           void writeEntries(Ostream& os) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "coupleManagerI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H

@@ -1,98 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-inline const Foam::fvPatch& Foam::coupleManager::patch() const
-{
-    return patch_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourRegionName() const
-{
-    return neighbourRegionName_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourPatchName() const
-{
-    return neighbourPatchName_;
-};
-
-
-inline const Foam::word& Foam::coupleManager::neighbourFieldName() const
-{
-     return neighbourFieldName_;
-};
-
-
-inline const Foam::fvMesh& Foam::coupleManager::neighbourRegion() const
-{
-    return localRegion_.objectRegistry::parent()
-        .lookupObject<fvMesh>(neighbourRegionName_);
-}
-
-
-inline Foam::label Foam::coupleManager::neighbourPatchID() const
-{
-    return neighbourRegion().boundaryMesh().findPatchID(neighbourPatchName_);
-}
-
-
-inline const Foam::fvPatch& Foam::coupleManager::neighbourPatch() const
-{
-    return neighbourRegion().boundary()[neighbourPatchID()];
-}
-
-
-template<class Type>
-inline const Foam::fvPatchField<Type>&
-Foam::coupleManager::neighbourPatchField() const
-{
-    return neighbourPatch().lookupPatchField
-        <GeometricField<Type, fvPatchField, volMesh>, Type>
-            (neighbourFieldName_);
-}
-
-
-inline Foam::word& Foam::coupleManager::neighbourRegionName()
-{
-    return neighbourRegionName_;
-};
-
-
-inline Foam::word& Foam::coupleManager::neighbourPatchName()
-{
-    return neighbourPatchName_;
-};
-
-
-inline Foam::word& Foam::coupleManager::neighbourFieldName()
-{
-    return neighbourFieldName_;
-};
-
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C

@@ -1,150 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-#include "solidWallTemperatureCoupledFvPatchScalarField.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    coupleManager_(p),
-    KName_("undefined-K")
-{}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
-    KName_(ptf.KName_)
-{}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
-    KName_(dict.lookup("K"))
-{
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
-solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-(
-    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& whftcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(whftcsf, iF),
-    coupleManager_(whftcsf.coupleManager_),
-    KName_(whftcsf.KName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<scalar>& neighbourField =
-        coupleManager_.neighbourPatchField<scalar>();
-
-    const fvPatchField<scalar>& K =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
-
-    gradient() = -refCast<const solidWallTemperatureCoupledFvPatchScalarField>
-        (neighbourField).flux()/K;
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H

@@ -1,165 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-
-Description
-    Fixed heat-flux boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallHeatFluxTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-SourceFiles
-    solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
-#define solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
-#include "coupleManager.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-   Class solidWallHeatFluxTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-:
-    public fixedGradientFvPatchScalarField
-{
-    // Private data
-
-        //- Couple manager object
-        coupleManager coupleManager_;
-
-        //- Name of thermal conductivity field
-        word KName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallHeatFluxTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        // solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-                    (
-                        *this
-                    )
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
-                    (
-                        *this,
-                        iF
-                    )
-            );
-        }
-
-
-    // Member functions
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C

@@ -28,8 +28,62 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
+#include "directMappedPatchBase.H"
 #include "regionProperties.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
+(
+    const polyMesh& nbrRegion
+) const
+{
+    const fvMesh& myRegion = patch().boundaryMesh().mesh();
+
+    const regionProperties& props =
+        myRegion.objectRegistry::parent().lookupObject<regionProperties>
+        (
+            "regionProperties"
+        );
+
+    label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
+    if (myIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), myRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn
+            (
+                "solidWallMixedTemperatureCoupledFvPatchScalarField"
+                "::interfaceOwner(const polyMesh&) const"
+            )   << "Cannot find region " << myRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        myIndex = props.fluidRegionNames().size() + i;
+    }
+    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
+    if (nbrIndex == -1)
+    {
+        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
+
+        if (i == -1)
+        {
+            FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
+                << "Cannot find region " << nbrRegion.name()
+                << " neither in fluids " << props.fluidRegionNames()
+                << " nor in solids " << props.solidRegionNames()
+                << exit(FatalError);
+        }
+        nbrIndex = props.fluidRegionNames().size() + i;
+    }
+
+    return myIndex < nbrIndex;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
@@ -40,7 +94,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    coupleManager_(p),
+    neighbourFieldName_("undefined-neighbourFieldName"),
     KName_("undefined-K")
 {
     this->refValue() = 0.0;
@@ -60,7 +114,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
+    neighbourFieldName_(ptf.neighbourFieldName_),
     KName_(ptf.KName_),
     fixesValue_(ptf.fixesValue_)
 {}
@@ -75,14 +129,46 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
+    neighbourFieldName_(dict.lookup("neighbourFieldName")),
     KName_(dict.lookup("K"))
 {
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+            "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
+            "(\n"
+            "    const fvPatch& p,\n"
+            "    const DimensionedField<scalar, volMesh>& iF,\n"
+            "    const dictionary& dict\n"
+            ")\n"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << directMappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << dimensionedInternalField().name()
+            << " in file " << dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
     fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
-    refValue() = static_cast<scalarField>(*this);
+
-    refGrad() = 0.0;
+    if (dict.found("refValue"))
-    valueFraction() = 1.0;
+    {
-    fixesValue_ = true;
+        // Full restart
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
+        fixesValue_ = readBool(dict.lookup("fixesValue"));
+    }
+    else
+    {
+        // Start from user entered data. Assume fixedValue.
+        refValue() = *this;
+        refGrad() = 0.0;
+        valueFraction() = 1.0;
+        fixesValue_ = true;
+    }
 }
 
 
@@ -94,7 +180,7 @@ solidWallMixedTemperatureCoupledFvPatchScalarField
 )
 :
     mixedFvPatchScalarField(wtcsf, iF),
-    coupleManager_(wtcsf.coupleManager_),
+    neighbourFieldName_(wtcsf.neighbourFieldName_),
     KName_(wtcsf.KName_),
     fixesValue_(wtcsf.fixesValue_)
 {}
@@ -116,33 +202,111 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
         return;
     }
 
+    // Get the coupling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        patch().patch()
+    );
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+
     tmp<scalarField> intFld = patchInternalField();
 
+    if (interfaceOwner(nbrMesh))
+    {
+        // Note: other side information could be cached - it only needs
+        // to be updated the first time round the iteration (i.e. when
+        // switching regions) but unfortunately we don't have this information.
+
+        const mapDistribute& distMap = mpp.map();
+        const fvPatch& nbrPatch = refCast<const fvMesh>
+        (
+            nbrMesh
+        ).boundary()[mpp.samplePolyPatch().index()];
+
+
+        // Calculate the temperature by harmonic averaging
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
+        refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
+        (
+            nbrPatch.lookupPatchField<volScalarField, scalar>
+            (
+                neighbourFieldName_
+            )
+        );
+
+        // Swap to obtain full local values of neighbour internal field
+        scalarField nbrIntFld = nbrField.patchInternalField();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrIntFld
+        );
+
+        // Swap to obtain full local values of neighbour K*delta
+        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            distMap.constructSize(),
+            distMap.subMap(),           // what to send
+            distMap.constructMap(),     // what to receive
+            nbrKDelta
+        );
+
+
+        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
+
+        // Calculate common wall temperature. Reuse *this to store common value.
+        scalarField Twall
+        (
+            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
+          / (myKDelta() + nbrKDelta)
+        );
+        // Assign to me
+        fvPatchScalarField::operator=(Twall);
+        // Distribute back and assign to neighbour
+        mapDistribute::distribute
+        (
+            Pstream::defaultCommsType,
+            distMap.schedule(),
+            nbrField.size(),
+            distMap.constructMap(),     // reverse : what to send
+            distMap.subMap(),
+            Twall
+        );
+        const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
+        (
+            nbrField
+        ).fvPatchScalarField::operator=(Twall);
+    }
+
+
+    // Switch between fixed value (of harmonic avg) or gradient
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     label nFixed = 0;
 
     // Like snGrad but bypass switching on refValue/refGrad.
-    tmp<scalarField> normalGradient =
+    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
-        (*this-intFld())
-      * patch().deltaCoeffs();
 
     if (debug)
     {
+        scalar Q = gSum(K()*patch().magSf()*normalGradient());
+
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
             << "updateCoeffs() :"
             << " patch:" << patch().name()
+            << " heatFlux:" << Q
             << " walltemperature "
-            << " min:"
+            << " min:" << gMin(*this)
-            <<  returnReduce
+            << " max:" << gMax(*this)
-                (
-                    (this->size() > 0 ? min(*this) : VGREAT),
-                    minOp<scalar>()
-                )
-            << " max:"
-            <<  returnReduce
-                (
-                    (this->size() > 0 ? max(*this) : -VGREAT),
-                    maxOp<scalar>()
-                )
             << " avg:" << gAverage(*this)
             << endl;
     }
@@ -150,7 +314,7 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
     forAll(*this, i)
     {
         // if outgoing flux use fixed value.
-        if (intFld()[i] > operator[](i))
+        if (normalGradient()[i] < 0.0)
         {
             this->refValue()[i] = operator[](i);
             this->refGrad()[i] = 0.0;   // not used
@@ -185,80 +349,16 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
 }
 
 
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
-(
-    const Pstream::commsTypes
-)
-{
-    if (!this->updated())
-    {
-        this->updateCoeffs();
-    }
-
-    if (!coupleManager_.regionOwner())
-    {
-        // I am the last one to evaluate.
-
-        tmp<scalarField> intFld = patchInternalField();
-
-        const fvPatch& nbrPatch = coupleManager_.neighbourPatch();
-
-        solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
-        refCast<solidWallMixedTemperatureCoupledFvPatchScalarField>
-        (
-            const_cast<fvPatchField<scalar>&>
-            (
-                coupleManager_.neighbourPatchField<scalar>()
-            )
-        );
-        tmp<scalarField> nbrIntFld = nbrField.patchInternalField();
-        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-        tmp<scalarField> nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
-
-        // Calculate common wall temperature and assign to both sides
-        scalarField::operator=
-        (
-            (myKDelta()*intFld + nbrKDelta()*nbrIntFld)
-          / (myKDelta() + nbrKDelta())
-        );
-
-        nbrField.scalarField::operator=(*this);
-
-        if (debug)
-        {
-            Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-                << "updateCoeffs() :"
-                << " patch:" << patch().name()
-                << " setting master and slave to wall temperature "
-                << " min:"
-                <<  returnReduce
-                    (
-                        (this->size() > 0 ? min(*this) : VGREAT),
-                        minOp<scalar>()
-                    )
-                << " max:"
-                <<  returnReduce
-                    (
-                        (this->size() > 0 ? max(*this) : -VGREAT),
-                        maxOp<scalar>()
-                    )
-                << " avg:" << gAverage(*this)
-                << endl;
-        }
-    }
-
-    fvPatchScalarField::evaluate();
-}
-
-
 void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
 (
     Ostream& os
 ) const
 {
     mixedFvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
+    os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
+        << token::END_STATEMENT << nl;
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
 }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

@@ -37,13 +37,22 @@ Description
         myInterfacePatchName
         {
             type                solidWallMixedTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
             neighbourFieldName  T;
             K                   K;
             value               uniform 300;
         }
 
+    Needs to be on underlying directMapped(Wall)FvPatch.
+
+    Note: runs in parallel with arbitrary decomposition. Uses directMapped
+    functionality to calculate exchange.
+
+    Note: lags interface data so both sides use same data.
+    - problem: schedule to calculate average would interfere
+    with standard processor swaps.
+    - so: updateCoeffs sets both to same Twall. Only need to do
+    this for last outer iteration but don't have access to this.
+
 SourceFiles
     solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
@@ -54,7 +63,6 @@ SourceFiles
 
 #include "fvPatchFields.H"
 #include "mixedFvPatchFields.H"
-#include "coupleManager.H"
 #include "fvPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,14 +80,21 @@ class solidWallMixedTemperatureCoupledFvPatchScalarField
 {
     // Private data
 
-        //- Couple manager object
+        //- Name of field on the neighbour region
-        coupleManager coupleManager_;
+        const word neighbourFieldName_;
-
+        
         //- Name of thermal conductivity field
-        word KName_;
+        const word KName_;
 
         bool fixesValue_;
 
+
+    // Private Member Functions
+
+        //- Am I or neighbour owner of interface
+        bool interfaceOwner(const polyMesh& nbrRegion) const;
+
+
 public:
 
     //- Runtime type information
@@ -162,12 +177,6 @@ public:
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
-        //- Evaluate the patch field
-        virtual void evaluate
-        (
-            const Pstream::commsTypes commsType=Pstream::blocking
-        );
-
         //- Write
         virtual void write(Ostream&) const;
 };

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C

@@ -1,156 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    coupleManager_(p),
-    KName_("undefined-K")
-{}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const solidWallTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    coupleManager_(ptf.coupleManager_),
-    KName_(ptf.KName_)
-{}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    coupleManager_(p, dict),
-    KName_(dict.lookup("K"))
-{
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::solidWallTemperatureCoupledFvPatchScalarField::
-solidWallTemperatureCoupledFvPatchScalarField
-(
-    const solidWallTemperatureCoupledFvPatchScalarField& wtcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(wtcsf, iF),
-    coupleManager_(wtcsf.coupleManager_),
-    KName_(wtcsf.KName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::solidWallTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const fvPatchField<scalar>& neighbourField =
-        coupleManager_.neighbourPatchField<scalar>();
-
-    operator==(neighbourField);
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    coupleManager_.writeEntries(os);
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::solidWallTemperatureCoupledFvPatchScalarField::flux() const
-{
-    const fvPatchScalarField& Kw =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
-
-    const fvPatchScalarField& Tw = *this;
-
-    return Tw.snGrad()*Kw;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H

@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    solidWallHeatFluxCoupledFvPatchScalarField
-
-Description
-    Fixed value boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallTemperatureCoupled;
-            neighbourRegionName fluid;
-            neighbourPatchName  fluidSolidInterface;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-SourceFiles
-    solidWallTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallTemperatureCoupledFvPatchScalarField_H
-#define solidWallTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
-#include "coupleManager.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-        Class solidWallTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallTemperatureCoupledFvPatchScalarField
-:
-    public fixedValueFvPatchScalarField
-{
-    // Private data
-
-        //- Couple manager object
-        coupleManager coupleManager_;
-
-        //- Name of thermal conductivity field
-        word KName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given solidWallTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallTemperatureCoupledFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallTemperatureCoupledFvPatchScalarField(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        //- (intensive) flux
-        tmp<scalarField> flux() const;
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H

@@ -14,12 +14,10 @@
          (
             UEqn()
          ==
-           -fvc::reconstruct
+            fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
-                    fvc::snGrad(pd)
+              - fvc::snGrad(p)*mesh.magSf()
-                  + ghf*fvc::snGrad(rho)
-                ) * mesh.magSf()
             )
         );
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -6,43 +6,15 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
     PtrList<volScalarField> DpDtFluid(fluidRegions.size());
-    PtrList<volScalarField> ghFluid(fluidRegions.size());
-    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
-    PtrList<volScalarField> pdFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
 
-    dimensionedScalar pRef
-    (
-        "pRef",
-        dimensionSet(1, -1, -2, 0, 0),
-        rp.lookup("pRef")
-    );
-
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
     {
         Info<< "*** Reading fluid mesh thermophysical properties for region "
             << fluidRegions[i].name() << nl << endl;
 
-        Info<< "    Adding to pdFluid\n" << endl;
-        pdFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "pd",
-                    runTime.timeName(),
-                    fluidRegions[i],
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                fluidRegions[i]
-            )
-        );
-
         Info<< "    Adding to thermoFluid\n" << endl;
         thermoFluid.set
         (
@@ -145,6 +117,7 @@
             i,
             new volScalarField
             (
+                "DpDt",
                 fvc::DDt
                 (
                     surfaceScalarField
@@ -157,36 +130,5 @@
             )
         );
 
-        const dictionary& environmentalProperties =
-            fluidRegions[i].lookupObject<IOdictionary>
-                ("environmentalProperties");
-        dimensionedVector g(environmentalProperties.lookup("g"));
-
-        Info<< "    Adding to ghFluid\n" << endl;
-        ghFluid.set
-        (
-            i,
-            new volScalarField
-            (
-                "gh",
-                 g & fluidRegions[i].C()
-            )
-        );
-        ghfFluid.set
-        (
-            i,
-            new surfaceScalarField
-            (
-                "ghf",
-                 g & fluidRegions[i].Cf()
-            )
-        );
-
-        Info<< "    Updating p from pd\n" << endl;
-        thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef;
-        thermoFluid[i].correct();
-
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
     }
-
-

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H

@@ -1,5 +1,5 @@
 {
-    tmp<fvScalarMatrix> hEqn
+    fvScalarMatrix hEqn
     (
         fvm::ddt(rho, h)
       + fvm::div(phi, h)
@@ -7,8 +7,16 @@
      ==
         DpDt
     );
-    hEqn().relax();
+    if (oCorr == nOuterCorr-1)
-    hEqn().solve();
+    {
+        hEqn.relax();
+        hEqn.solve(mesh.solver("hFinal"));
+    }
+    else
+    {
+        hEqn.relax();
+        hEqn.solve();
+    }
 
     thermo.correct();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -1,5 +1,5 @@
 {
-    bool closedVolume = pd.needReference();
+    bool closedVolume = p.needReference();
 
     rho = thermo.rho();
 
@@ -17,31 +17,34 @@
         )
     );
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+    phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rhorUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, pd)
+            fvm::ddt(psi, p)
-          + fvc::ddt(psi)*pRef
-          + fvc::ddt(psi, rho)*gh
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, pd)
+          - fvm::laplacian(rhorUAf, p)
         );
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        if
+        (
+            oCorr == nOuterCorr-1
+         && corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
@@ -49,27 +52,24 @@
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
-    // Update pressure field (including bc)
+    // Update pressure substantive derivative
-    p == pd + rho*gh + pRef;
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     // Solve continuity
-#   include "rhoEqn.H"
+    #include "rhoEqn.H"
 
     // Update continuity errors
-#   include "compressibleContinuityErrors.H"
+    #include "compressibleContinuityErrors.H"
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi);
+        p += (massIni - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
         rho = thermo.rho();
     }
 
     // Update thermal conductivity
     K = thermoFluid[i].Cp()*turb.alphaEff();
-
-    // Update pd (including bc)
-    pd == p - (rho*gh + pRef);
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPIMPLEControls.H

@@ -0,0 +1,9 @@
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+
+    int nCorr(readInt(pimple.lookup("nCorrectors")));
+
+    int nNonOrthCorr =
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
+
+    bool momentumPredictor =
+        pimple.lookupOrDefault<Switch>("momentumPredictor", true);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H

@@ -1 +0,0 @@
-    solve(fvm::ddt(rho) + fvc::div(phi));

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -7,12 +7,15 @@
     surfaceScalarField& phi = phiFluid[i];
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
-    const volScalarField& gh = ghFluid[i];
-    const surfaceScalarField& ghf = ghfFluid[i];
-    volScalarField& pd = pdFluid[i];
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
     const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
+
+    const dictionary& environmentalProperties =
+        fluidRegions[i].lookupObject<IOdictionary>
+        ("environmentalProperties");
+
+    const dimensionedVector g(environmentalProperties.lookup("g"));

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H

@@ -1,15 +1,18 @@
-#   include "rhoEqn.H"
+if (oCorr == 0)
-    for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
+{
-    {
+    #include "rhoEqn.H"
-    #   include "UEqn.H"
+}
 
-    #   include "hEqn.H"
+#include "UEqn.H"
 
-     // --- PISO loop
+#include "hEqn.H"
 
-        for (int corr=0; corr<nCorr; corr++)
+// --- PISO loop
-        {
+for (int corr=0; corr<nCorr; corr++)
-    #       include "pEqn.H"
+{
-        }
+    #include "pEqn.H"
-    }
+}
-    turb.correct();
+
+turb.correct();
+
+rho = thermo.rho();

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H

@@ -0,0 +1,2 @@
+    p.storePrevIter();
+    rho.storePrevIter();

applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H

@@ -0,0 +1,7 @@
+    // We do not have a top-level mesh. Construct the fvSolution for
+    // the runTime instead.
+    fvSolution solutionDict(runTime);
+
+    const dictionary& pimple = solutionDict.subDict("PIMPLE");
+
+    int nOuterCorr(readInt(pimple.lookup("nOuterCorrectors")));

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPIMPLEControls.H

@@ -0,0 +1,4 @@
+    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
+
+    int nNonOrthCorr =
+        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H

@@ -1,4 +1,4 @@
-    dictionary piso = solidRegions[i].solutionDict().subDict("PISO");
+    const dictionary& piso = solidRegions[i].solutionDict().subDict("PISO");
 
     int nNonOrthCorr = 0;
     if (piso.found("nNonOrthogonalCorrectors"))

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -1,4 +1,4 @@
-//    fvMesh& mesh = solidRegions[i];
+    fvMesh& mesh = solidRegions[i];
 
     volScalarField& rho = rhos[i];
     volScalarField& cp = cps[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -1,10 +1,13 @@
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        solve
+        tmp<fvScalarMatrix> TEqn
         (
-            fvm::ddt(rho*cp, T) - fvm::laplacian(K, T)
+            fvm::ddt(rho*cp, T)
+          - fvm::laplacian(K, T)
         );
+        TEqn().relax();
+        TEqn().solve();
     }
 
     Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;

applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p)
-                  - fvc::snGrad(pd)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H

@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 

applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -88,24 +88,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
@@ -152,11 +134,11 @@
     );
 
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }

applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H

@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
+        pEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
+        pEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
-          - ghf*fvc::snGrad(rho)
+          + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         if
@@ -51,9 +51,9 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
+               *pEqnComp()
-              + pdEqnIncomp,
+              + pEqnIncomp,
-                mesh.solver(pd.name() + "Final")
+                mesh.solver(p.name() + "Final")
             );
         }
         else
@@ -64,8 +64,8 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
+               *pEqnComp()
-              + pdEqnIncomp
+              + pEqnIncomp
             );
         }
 
@@ -73,26 +73,21 @@
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
+               *(pEqnComp & p);
-            phi += pdEqnIncomp.flux();
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
+    p.max(pMin);
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p)
-                  - fvc::snGrad(pd)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -46,11 +46,6 @@
     #include "createPhi.H"
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
@@ -88,24 +83,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
+        pEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
+        pEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
-          - ghf*fvc::snGrad(rho)
+          + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         solve
@@ -44,31 +44,27 @@
                 max(alpha1, scalar(0))*(psi1/rho1)
               + max(alpha2, scalar(0))*(psi2/rho2)
             )
-           *pdEqnComp()
+           *pEqnComp()
-          + pdEqnIncomp
+          + pEqnIncomp
         );
 
         if (nonOrth == nNonOrthCorr)
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
+               *(pEqnComp & p);
-            phi += pdEqnIncomp.flux();
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
+    p.max(pMin);
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 }

applications/solvers/multiphase/interDyMFoam/correctPhi.H

@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -27,7 +27,7 @@
             fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)

applications/solvers/multiphase/interDyMFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -92,47 +92,17 @@
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
     }
 
 
-    volScalarField p
+    label pRefCell = 0;
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*(g & mesh.C())
-    );
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
     scalar pRefValue = 0.0;
-
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }

applications/solvers/multiphase/interDyMFoam/interDyMFoam.C

@@ -114,18 +114,6 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -7,38 +7,38 @@
 
     phi = phiU +
     (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
-      - ghf*fvc::snGrad(rho)
+      + fvc::interpolate(rho)*(g & mesh.Sf())
-    )*rAUf*mesh.magSf();
+    )*rAUf;
 
-    if (pd.needReference())
+    if (p.needReference())
     {
         fvc::makeRelative(phi, U);
-        adjustPhi(phi, U, pd);
+        adjustPhi(phi, U, p);
         fvc::makeAbsolute(phi, U);
     }
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUf, pd) == fvc::div(phi)
+            fvm::laplacian(rAUf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 

applications/solvers/multiphase/interFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p)
-                  - fvc::snGrad(pd)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -1,11 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
@@ -22,7 +22,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -37,7 +37,7 @@
             fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)

applications/solvers/multiphase/interFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -83,45 +83,9 @@
     );
 
 
-    Info<< "Calculating field g.h\n" << endl;
+    label pRefCell = 0;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
     scalar pRefValue = 0.0;
-
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
 
 
     // Construct interface from alpha1 distribution

applications/solvers/multiphase/interFoam/interFoam.C

@@ -89,18 +89,6 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/interFoam/pEqn.H

@@ -12,33 +12,33 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
-          - ghf*fvc::snGrad(rho)
+          + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
-    adjustPhi(phi, U, pd);
+    adjustPhi(phi, U, p);
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -25,10 +25,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p)
-                  - fvc::snGrad(pd)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H

@@ -1,11 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
@@ -22,7 +22,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -37,7 +37,7 @@
             fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -66,26 +66,9 @@
     rho.oldTime();
 
 
-    label pdRefCell = 0;
+    label pRefCell = 0;
-    scalar pdRefValue = 0.0;
+    scalar pRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-    Info<< "Calculating field g.h" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
 
 
     // Construct interface from alpha1 distribution

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -13,11 +13,11 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
-          - ghf*fvc::snGrad(rho)
+          + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
-    adjustPhi(phi, U, pd);
+    adjustPhi(phi, U, p);
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -25,31 +25,29 @@
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvc::div(phi) - fvm::laplacian(rUAf, pd)
+            fvc::div(phi) - fvm::laplacian(rUAf, p)
-          + (vDotvP - vDotcP)*(rho*gh - pSat) + fvm::Sp(vDotvP - vDotcP, pd)
+          - (vDotvP - vDotcP)*pSat + fvm::Sp(vDotvP - vDotcP, p)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
-    p = pd + rho*gh;
-
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 }

applications/solvers/multiphase/multiphaseInterFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     mixture.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p)
-                  - fvc::snGrad(pd)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/multiphaseInterFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -45,45 +45,9 @@
     rho.oldTime();
 
 
-    Info<< "Calculating field g.h\n" << endl;
+    label pRefCell = 0;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
     scalar pRefValue = 0.0;
-
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
 
 
     // Construct incompressible turbulence model

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C

@@ -81,18 +81,6 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -12,33 +12,33 @@
 
     phi = phiU +
         (
-            mixture.surfaceTensionForce()
+            mixture.surfaceTensionForce()*mesh.magSf()
-          - ghf*fvc::snGrad(rho)
+          + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
-    adjustPhi(phi, U, pd);
+    adjustPhi(phi, U, p);
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 

applications/solvers/multiphase/settlingFoam/UEqn.H

@@ -22,10 +22,8 @@
           ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
-                  - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p)*mesh.magSf()
-                  - fvc::snGrad(p)
-                )*mesh.magSf()
             )
         );
     }

applications/solvers/multiphase/settlingFoam/createFields.H

@@ -337,6 +337,3 @@
         ),
         mut + mul
     );
-
-    Info<< "Calculating field (g.h)f\n" << endl;
-    surfaceScalarField ghf = surfaceScalarField("ghf", g & mesh.Cf());

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -15,7 +15,7 @@ phi =
     );
 
 surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf();
+phi += fvc::interpolate(rho)*(g & mesh.Sf())*rUAf;
 
 for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {

applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H

@@ -14,7 +14,14 @@
     //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
-    UEqn.relax();
+    if (oCorr == nOuterCorr-1)
+    {
+        UEqn.relax(1);
+    }
+    else
+    {
+        UEqn.relax();
+    }
 
     if (momentumPredictor)
     {
@@ -22,9 +29,11 @@
         (
             UEqn
          ==
-           -fvc::reconstruct
+            fvc::reconstruct
             (
-                mesh.magSf()*(fvc::snGrad(pd) + ghf*fvc::snGrad(rho))
+                fvc::interpolate(rho)*(g & mesh.Sf())
-            )
+              - mesh.magSf()*fvc::snGrad(p)
+            ),
+            mesh.solver(oCorr == nOuterCorr-1 ? "UFinal" : "U")
         );
     }

applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H

@@ -3,7 +3,12 @@
     (
         fvm::ddt(alpha1)
       + fvm::div(phi, alpha1)
-      - fvm::laplacian(Dab, alpha1)
+      //- fvm::Sp(fvc::div(phi), alpha1)
+      - fvm::laplacian
+        (
+            Dab + alphatab*turbulence->nut(), alpha1,
+            "laplacian(Dab,alpha1)"
+        )
     );
 
     alpha1Eqn.solve();

applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
-    volScalarField pd
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -50,6 +50,9 @@
 
     dimensionedScalar Dab(twoPhaseProperties.lookup("Dab"));
 
+    // Read the reciprocal of the turbulent Schmidt number
+    dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab"));
+
     // Need to store rho for ddt(rho, U)
     volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2);
     rho.oldTime();
@@ -72,45 +75,9 @@
     );
 
 
-    Info<< "Calculating field g.h\n" << endl;
+    label pRefCell = 0;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("gh", g & mesh.Cf());
-
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
     scalar pRefValue = 0.0;
-
+    setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
 
 
     // Construct incompressible turbulence model

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -7,26 +7,37 @@
     surfaceScalarField phiU
     (
         "phiU",
-        (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rUAf*mesh.magSf();
+    phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rUAf;
-
-    adjustPhi(phi, U, pd);
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
-        pdEqn.solve();
+
+        if
+        (
+            corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver(p.name() + "Final"));
+        }
+        else
+        {
+            pEqn.solve(mesh.solver(p.name()));
+        }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 

applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C

@@ -40,53 +40,52 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPIMPLEControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
-#   include "setRootCase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (runTime.run())
     {
+        #include "readPIMPLEControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
+        // --- Pressure-velocity PIMPLE corrector loop
-#       include "CourantNo.H"
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
-
-        twoPhaseProperties.correct();
-
-#       include "alphaEqn.H"
-
-#       include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            twoPhaseProperties.correct();
+
+            #include "alphaEqn.H"
+
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            #include "continuityErrs.H"
+
+            turbulence->correct();
         }
 
-#       include "continuityErrs.H"
-
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/test/extendedStencil/Make/files

@@ -0,0 +1,4 @@
+
+testExtendedStencil.C
+
+EXE = $(FOAM_USER_APPBIN)/testExtendedStencil

applications/test/extendedStencil/Make/options

@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume

applications/test/extendedStencil/testExtendedStencil.C

@@ -0,0 +1,499 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    testExtendedStencil
+
+Description
+    Test app for determining extended stencil.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "Time.H"
+#include "mapDistribute.H"
+#include "OFstream.H"
+#include "meshTools.H"
+//#include "FECCellToFaceStencil.H"
+//#include "CFCCellToFaceStencil.H"
+//#include "CPCCellToFaceStencil.H"
+//#include "CECCellToFaceStencil.H"
+//#include "extendedCentredCellToFaceStencil.H"
+//#include "extendedUpwindCellToFaceStencil.H"
+
+//#include "centredCFCCellToFaceStencilObject.H"
+//#include "centredFECCellToFaceStencilObject.H"
+//#include "centredCPCCellToFaceStencilObject.H"
+//#include "centredCECCellToFaceStencilObject.H"
+
+//#include "upwindFECCellToFaceStencilObject.H"
+//#include "upwindCPCCellToFaceStencilObject.H"
+//#include "upwindCECCellToFaceStencilObject.H"
+
+//#include "upwindCFCCellToFaceStencilObject.H"
+#include "centredCFCFaceToCellStencilObject.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void writeStencilOBJ
+(
+    const fileName& fName,
+    const point& fc,
+    const List<point>& stencilCc
+)
+{
+    OFstream str(fName);
+    label vertI = 0;
+
+    meshTools::writeOBJ(str, fc);
+    vertI++;
+
+    forAll(stencilCc, i)
+    {
+        meshTools::writeOBJ(str, stencilCc[i]);
+        vertI++;
+        str << "l 1 " << vertI << nl;
+    }
+}
+
+
+// Stats
+void writeStencilStats(const labelListList& stencil)
+{
+    label sumSize = 0;
+    label nSum = 0;
+    label minSize = labelMax;
+    label maxSize = labelMin;
+
+    forAll(stencil, i)
+    {
+        const labelList& sCells = stencil[i];
+
+        if (sCells.size() > 0)
+        {
+            sumSize += sCells.size();
+            nSum++;
+            minSize = min(minSize, sCells.size());
+            maxSize = max(maxSize, sCells.size());
+        }
+    }
+    reduce(sumSize, sumOp<label>());
+    reduce(nSum, sumOp<label>());
+    sumSize /= nSum;
+
+    reduce(minSize, minOp<label>());
+    reduce(maxSize, maxOp<label>());
+
+    Info<< "Stencil size :" << nl
+        << "    average : " << sumSize << nl
+        << "    min     : " << minSize << nl
+        << "    max     : " << maxSize << nl
+        << endl;
+}
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addTimeOptions.H"
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    // Get times list
+    instantList Times = runTime.times();
+#   include "checkTimeOptions.H"
+    runTime.setTime(Times[startTime], startTime);
+#   include "createMesh.H"
+
+    // Force calculation of extended edge addressing
+    const labelListList& edgeFaces = mesh.edgeFaces();
+    const labelListList& edgeCells = mesh.edgeCells();
+    const labelListList& pointCells = mesh.pointCells();
+    Info<< "dummy:" << edgeFaces.size() + edgeCells.size() + pointCells.size()
+        << endl;
+
+
+    // Centred, semi-extended stencil (edge cells only)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+//    {
+//        //const FECCellToFaceStencil cfcStencil(mesh);
+//        //const extendedCentredCellToFaceStencil addressing
+//        //(
+//        //    cfcStencil
+//        //);
+//        const extendedCentredStencil& addressing =
+//        centredFECCellToFaceStencilObject::New
+//        (
+//            mesh
+//        );
+//
+//        Info<< "faceEdgeCell:" << endl;
+//        writeStencilStats(addressing.stencil());
+//
+//        // Collect stencil cell centres
+//        List<List<point> > stencilPoints(mesh.nFaces());
+//        addressing.collectData
+//        (
+//            mesh.C(),
+//            stencilPoints
+//        );
+//
+//        forAll(stencilPoints, faceI)
+//        {
+//            writeStencilOBJ
+//            (
+//                runTime.path()/"faceEdgeCell" + Foam::name(faceI) + ".obj",
+//                mesh.faceCentres()[faceI],
+//                stencilPoints[faceI]
+//            );
+//        }
+//    }
+
+
+
+
+//    // Centred, face stencil
+//    // ~~~~~~~~~~~~~~~~~~~~~
+//
+//    {
+//        const extendedCentredCellToFaceStencil& addressing =
+//        centredCFCCellToFaceStencilObject::New
+//        (
+//            mesh
+//        );
+//        
+//        Info<< "cellFaceCell:" << endl;
+//        writeStencilStats(addressing.stencil());
+//        
+//        
+//        //// Do some interpolation.
+//        //{
+//        //    const labelListList& stencil = addressing.stencil();
+//        //    List<List<scalar> > stencilWeights(stencil.size());
+//        //    forAll(stencil, faceI)
+//        //    {
+//        //        const labelList& fStencil = stencil[faceI];
+//        //
+//        //        if (fStencil.size() > 0)
+//        //        {
+//        //            // Uniform weights
+//        //            stencilWeights[faceI] = scalarList
+//        //            (
+//        //                fStencil.size(),
+//        //                1.0/fStencil.size()
+//        //            );
+//        //        }
+//        //    }
+//        //
+//        //    tmp<surfaceVectorField> tfc
+//        //    (
+//        //        addressing.weightedSum(mesh.C(), stencilWeights)
+//        //    );
+//        //}
+//
+//
+//        // Collect stencil cell centres
+//        List<List<point> > stencilPoints(mesh.nFaces());
+//        addressing.collectData
+//        (
+//            mesh.C(),
+//            stencilPoints
+//        );
+//        
+//        forAll(stencilPoints, faceI)
+//        {
+//            if (stencilPoints[faceI].size() >= 15)
+//            {
+//                writeStencilOBJ
+//                (
+//                    runTime.path()/"centredFace" + Foam::name(faceI) + ".obj",
+//                    mesh.faceCentres()[faceI],
+//                    stencilPoints[faceI]
+//                );
+//            }
+//        }
+//    }
+
+
+//    // Centred, point stencil
+//    // ~~~~~~~~~~~~~~~~~~~~~~
+//
+//    {
+//        //const extendedCentredCellToFaceStencil& addressing =
+//        //centredCPCStencilObject::New
+//        //(
+//        //    mesh
+//        //);
+//        //
+//        //Info<< "cellPointCell:" << endl;
+//        //writeStencilStats(addressing.stencil());
+//        //
+//        //
+//        //// Collect stencil cell centres
+//        //List<List<point> > stencilPoints(mesh.nFaces());
+//        //addressing.collectData
+//        //(
+//        //    mesh.C(),
+//        //    stencilPoints
+//        //);
+//        //
+//        //forAll(stencilPoints, faceI)
+//        //{
+//        //    writeStencilOBJ
+//        //    (
+//        //        runTime.path()/"centredPoint" + Foam::name(faceI) + ".obj",
+//        //        mesh.faceCentres()[faceI],
+//        //        stencilPoints[faceI]
+//        //    );
+//        //}
+//    }
+
+
+
+//    // Centred, edge stencil
+//    // ~~~~~~~~~~~~~~~~~~~~~~
+//
+//    {
+//        //const extendedCentredCellToFaceStencil& addressing =
+//        //centredCECStencilObject::New
+//        //(
+//        //    mesh
+//        //);
+//        //
+//        //Info<< "cellEdgeCell:" << endl;
+//        //writeStencilStats(addressing.stencil());
+//        //
+//        //
+//        //// Collect stencil cell centres
+//        //List<List<point> > stencilPoints(mesh.nFaces());
+//        //addressing.collectData
+//        //(
+//        //    mesh.C(),
+//        //    stencilPoints
+//        //);
+//        //
+//        //forAll(stencilPoints, faceI)
+//        //{
+//        //    writeStencilOBJ
+//        //    (
+//        //        runTime.path()/"centredEdge" + Foam::name(faceI) + ".obj",
+//        //        mesh.faceCentres()[faceI],
+//        //        stencilPoints[faceI]
+//        //    );
+//        //}
+//    }
+
+
+
+    // Upwind, semi-extended stencil
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //{
+    //    const extendedUpwindCellToFaceStencil& addressing =
+    //    upwindFECCellToFaceStencilObject::New
+    //    (
+    //        mesh,
+    //        0.5
+    //    );
+    //
+    //    Info<< "upwind-faceEdgeCell:" << endl;
+    //    writeStencilStats(addressing.ownStencil());
+    //
+    //    {
+    //        // Collect stencil cell centres
+    //        List<List<point> > ownPoints(mesh.nFaces());
+    //        addressing.collectData
+    //        (
+    //            addressing.ownMap(),
+    //            addressing.ownStencil(),
+    //            mesh.C(),
+    //            ownPoints
+    //        );
+    //
+    //        forAll(ownPoints, faceI)
+    //        {
+    //            writeStencilOBJ
+    //            (
+    //                runTime.path()/"ownFEC" + Foam::name(faceI) + ".obj",
+    //                mesh.faceCentres()[faceI],
+    //                ownPoints[faceI]
+    //            );
+    //        }
+    //    }
+    //    {
+    //        // Collect stencil cell centres
+    //        List<List<point> > neiPoints(mesh.nFaces());
+    //        addressing.collectData
+    //        (
+    //            addressing.neiMap(),
+    //            addressing.neiStencil(),
+    //            mesh.C(),
+    //            neiPoints
+    //        );
+    //
+    //        forAll(neiPoints, faceI)
+    //        {
+    //            writeStencilOBJ
+    //            (
+    //                runTime.path()/"neiFEC" + Foam::name(faceI) + ".obj",
+    //                mesh.faceCentres()[faceI],
+    //                neiPoints[faceI]
+    //            );
+    //        }
+    //    }
+    //}
+
+
+
+    // Upwind, extended stencil
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //{
+    //    const extendedUpwindCellToFaceStencil& addressing =
+    //    upwindCFCCellToFaceStencilObject::New
+    //    (
+    //        mesh,
+    //        0.5
+    //    );
+    //
+    //    Info<< "upwind-cellFaceCell:" << endl;
+    //    writeStencilStats(addressing.ownStencil());
+    //
+    //    {
+    //        // Collect stencil cell centres
+    //        List<List<point> > ownPoints(mesh.nFaces());
+    //        addressing.collectData
+    //        (
+    //            addressing.ownMap(),
+    //            addressing.ownStencil(),
+    //            mesh.C(),
+    //            ownPoints
+    //        );
+    //
+    //        forAll(ownPoints, faceI)
+    //        {
+    //            writeStencilOBJ
+    //            (
+    //                runTime.path()/"ownCFC" + Foam::name(faceI) + ".obj",
+    //                mesh.faceCentres()[faceI],
+    //                ownPoints[faceI]
+    //            );
+    //        }
+    //    }
+    //    {
+    //        // Collect stencil cell centres
+    //        List<List<point> > neiPoints(mesh.nFaces());
+    //        addressing.collectData
+    //        (
+    //            addressing.neiMap(),
+    //            addressing.neiStencil(),
+    //            mesh.C(),
+    //            neiPoints
+    //        );
+    //
+    //        forAll(neiPoints, faceI)
+    //        {
+    //            writeStencilOBJ
+    //            (
+    //                runTime.path()/"neiCFC" + Foam::name(faceI) + ".obj",
+    //                mesh.faceCentres()[faceI],
+    //                neiPoints[faceI]
+    //            );
+    //        }
+    //    }
+    //}
+
+
+
+    //---- CELL CENTRED STENCIL -----
+
+    // Centred, cell stencil
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        const extendedCentredFaceToCellStencil& addressing =
+        centredCFCFaceToCellStencilObject::New
+        (
+            mesh
+        );
+        
+        Info<< "cellFaceCell:" << endl;
+        writeStencilStats(addressing.stencil());
+        
+        // Collect stencil face centres
+        List<List<point> > stencilPoints(mesh.nCells());
+        addressing.collectData
+        (
+            mesh.Cf(),
+            stencilPoints
+        );
+        
+        forAll(stencilPoints, cellI)
+        {
+            writeStencilOBJ
+            (
+                runTime.path()/"centredCell" + Foam::name(cellI) + ".obj",
+                mesh.cellCentres()[cellI],
+                stencilPoints[cellI]
+            );
+        }
+    }
+
+
+//XXXXXX
+//    // Evaluate
+//    List<List<scalar> > stencilData(faceStencils.size());
+//    collectStencilData
+//    (
+//        distMap,
+//        faceStencils,
+//        vf,
+//        stencilData
+//    );
+//    for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
+//    {
+//        const scalarList& stData = stencilData[faceI];
+//        const scalarList& stWeight = fit[faceI];
+//
+//        forAll(stData, i)
+//        {
+//            sf[faceI] += stWeight[i]*stData[i];
+//        }
+//    }
+//    See finiteVolume/lnInclude/leastSquaresGrad.C
+//XXXXXX
+
+    Pout<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -566,7 +566,7 @@ int main(int argc, char *argv[])
             IOobject
             (
                 "cellProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::MUST_READ,
@@ -579,7 +579,7 @@ int main(int argc, char *argv[])
             IOobject
             (
                 "boundaryProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::MUST_READ,
@@ -603,7 +603,7 @@ int main(int argc, char *argv[])
                 IOobject
                 (
                     "faceProcAddressing",
-                    "constant",
+                    procMesh.facesInstance(),
                     procMesh.meshSubDir,
                     procMesh,
                     IOobject::MUST_READ,
@@ -645,7 +645,7 @@ int main(int argc, char *argv[])
                 IOobject
                 (
                     "pointProcAddressing",
-                    "constant",
+                    procMesh.facesInstance(),
                     procMesh.meshSubDir,
                     procMesh,
                     IOobject::MUST_READ,

applications/utilities/parallelProcessing/decomposePar/decomposeParDict

@@ -55,7 +55,17 @@ metisCoeffs
 }
 
 scotchCoeffs
-{}
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
 
 manualCoeffs
 {

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -269,7 +269,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 this->polyMesh::name(),  // region name of undecomposed mesh
-                "constant",
+                pointsInstance(),
                 processorDb
             ),
             xferMove(procPoints),
@@ -620,7 +620,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "pointProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,
@@ -635,7 +635,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "faceProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,
@@ -650,7 +650,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "cellProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,
@@ -665,7 +665,7 @@ bool domainDecomposition::writeDecomposition()
             IOobject
             (
                 "boundaryProcAddressing",
-                "constant",
+                procMesh.facesInstance(),
                 procMesh.meshSubDir,
                 procMesh,
                 IOobject::NO_READ,

applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C

@@ -56,13 +56,12 @@ processorVolPatchFieldDecomposer
     const unallocLabelList& addressingSlice
 )
 :
-    addressing_(addressingSlice.size()),
+    directAddressing_(addressingSlice.size())
-    weights_(addressingSlice.size())
 {
     const labelList& own = mesh.faceOwner();
     const labelList& neighb = mesh.faceNeighbour();
 
-    forAll (addressing_, i)
+    forAll (directAddressing_, i)
     {
         // Subtract one to align addressing.  
         label ai = mag(addressingSlice[i]) - 1;
@@ -74,18 +73,14 @@ processorVolPatchFieldDecomposer
             // on the parallel boundary.
             // Give face the value of the neighbour.
 
-            addressing_[i].setSize(1);
-            weights_[i].setSize(1);
-            weights_[i][0] = 1.0;
-
             if (addressingSlice[i] >= 0)
             {
                 // I have the owner so use the neighbour value
-                addressing_[i][0] = neighb[ai];
+                directAddressing_[i] = neighb[ai];
             }
             else
             {
-                addressing_[i][0] = own[ai];
+                directAddressing_[i] = own[ai];
             }
         }
         else
@@ -96,12 +91,7 @@ processorVolPatchFieldDecomposer
             // up the different (face) list of data), so I will
             // just grab the value from the owner cell
 
-            addressing_[i].setSize(1);
+            directAddressing_[i] = own[ai];
-            weights_[i].setSize(1);
-
-            addressing_[i][0] = own[ai];
-
-            weights_[i][0] = 1.0;
         }
     }
 }

applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H

@@ -96,15 +96,16 @@ public:
         };
 
 
-        //- Processor patch field decomposer class
+        //- Processor patch field decomposer class. Maps either owner or
+        //  neighbour data (no interpolate anymore - processorFvPatchField
+        //  holds neighbour data)
         class processorVolPatchFieldDecomposer
         :
             public fvPatchFieldMapper
         {
             // Private data
 
-                labelListList addressing_;
+                labelList directAddressing_;
-                scalarListList weights_;
 
         public:
 
@@ -120,27 +121,23 @@ public:
 
                 label size() const
                 {
-                    return addressing_.size();
+                    return directAddressing_.size();
                 }
 
                 bool direct() const
                 {
-                    return false;
+                    return true;
                 }
 
-                const labelListList& addressing() const
+                const unallocLabelList& directAddressing() const
                 {
-                    return addressing_;
+                    return directAddressing_;
-                }
-
-                const scalarListList& weights() const
-                {
-                    return weights_;
                 }
         };
 
 
-        //- Processor patch field decomposer class
+        //- Processor patch field decomposer class. Surface field is assumed
+        //  to have direction (so manipulates sign when mapping)
         class processorSurfacePatchFieldDecomposer
         :
             public fvPatchFieldMapper

applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C

@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
             // Give patch area
             if (isType<cyclicPolyPatch>(mesh.boundaryMesh()[patchi]))
             {
-                Info<< "    Cyclic patch area: " << nl;
+                Info<< "    Cyclic patch vector area: " << nl;
                 label nFaces = mesh.boundaryMesh()[patchi].size();
                 vector sum1 = vector::zero;
                 vector sum2 = vector::zero;
@@ -92,12 +92,18 @@ int main(int argc, char *argv[])
                 Info<< "    - half 1 = " << sum1 << ", " << mag(sum1) << nl
                     << "    - half 2 = " << sum2 << ", " << mag(sum2) << nl
                     << "    - total  = " << (sum1 + sum2) << ", "
-                    << mag(sum1 + sum2) << endl;;
+                    << mag(sum1 + sum2) << endl;
+                Info<< "    Cyclic patch area magnitude = "
+                    << gSum(mesh.magSf().boundaryField()[patchi])/2.0 << endl;
             }
             else
             {
-                Info<< "    Patch area = "
+                Info<< "    Area vector of patch "
+                    << patchName << '[' << patchi << ']' << " = "
                     << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
+                Info<< "    Area magnitude of patch "
+                    << patchName << '[' << patchi << ']' << " = "
+                    << gSum(mesh.magSf().boundaryField()[patchi]) << endl;
             }
 
             // Read field and calc integral
@@ -107,15 +113,26 @@ int main(int argc, char *argv[])
                     << fieldName << endl;
 
                 volScalarField field(fieldHeader, mesh);
-                vector sumField = gSum
-                (
-                    mesh.Sf().boundaryField()[patchi]
-                   *field.boundaryField()[patchi]
-                );
 
-                Info<< "    Integral of " << fieldName << " over patch "
+                Info<< "    Integral of " << fieldName
+                    << " over vector area of patch "
                     << patchName << '[' << patchi << ']' << " = "
-                    << sumField << nl;
+                    << gSum
+                       (
+                           mesh.Sf().boundaryField()[patchi]
+                          *field.boundaryField()[patchi]
+                       )
+                    << nl;
+
+                Info<< "    Integral of " << fieldName
+                    << " over area magnitude of patch "
+                    << patchName << '[' << patchi << ']' << " = "
+                    << gSum
+                       (
+                           mesh.magSf().boundaryField()[patchi]
+                          *field.boundaryField()[patchi]
+                       )
+                    << nl;
             }
             else if
             (

applications/utilities/postProcessing/wall/yPlusLES/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/LES/LESModel \
     -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \

applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C

@@ -34,6 +34,7 @@ Description
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 #include "LESModel.H"
 #include "nearWallDist.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,18 @@ int main(int argc, char *argv[])
     {
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
-        mesh.readUpdate();
+        fvMesh::readUpdateState state = mesh.readUpdate();
+
+        // Wall distance
+        if (timeI == 0 || state != fvMesh::UNCHANGED)
+        {
+            Info<< "Calculating wall distance\n" << endl;
+            wallDist y(mesh, true);
+            Info<< "Writing wall distance to field "
+                << y.name() << nl << endl;
+            y.write();
+        }
+
 
         volScalarField yPlus
         (
@@ -116,6 +128,9 @@ int main(int argc, char *argv[])
             }
         }
 
+        Info<< "Writing yPlus to field "
+            << yPlus.name() << nl << endl;
+
         yPlus.write();
     }
 

applications/utilities/postProcessing/wall/yPlusRAS/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \

applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.C

@@ -34,6 +34,7 @@ Description
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "wallFvPatch.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,17 @@ int main(int argc, char *argv[])
     {
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
-        mesh.readUpdate();
+        fvMesh::readUpdateState state = mesh.readUpdate();
+
+        // Wall distance
+        if (timeI == 0 || state != fvMesh::UNCHANGED)
+        {
+            Info<< "Calculating wall distance\n" << endl;
+            wallDist y(mesh, true);
+            Info<< "Writing wall distance to field "
+                << y.name() << nl << endl;
+            y.write();
+        }
 
         volScalarField yPlus
         (
@@ -106,6 +117,9 @@ int main(int argc, char *argv[])
             }
         }
 
+        Info<< "Writing yPlus to field "
+            << yPlus.name() << nl << endl;
+
         yPlus.write();
     }
 

applications/utilities/surface/surfaceRedistributePar/Make/files

@@ -0,0 +1,3 @@
+surfaceRedistributePar.C
+
+EXE = $(FOAM_APPBIN)/surfaceRedistributePar

applications/utilities/surface/surfaceRedistributePar/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -ltriSurface

applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C

@@ -0,0 +1,295 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    surfaceRedistributePar
+
+Description
+    (Re)distribution of triSurface. Either takes an undecomposed surface
+    or an already decomposed surface and redistribute it so each processor
+    has all triangles that overlap its mesh.
+
+Note
+    - best decomposition option is hierarchGeomDecomp since
+    guarantees square decompositions.
+    - triangles might be present on multiple processors.
+    - merging uses geometric tolerance so take care with writing precision.
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeBoundBox.H"
+#include "FixedList.H"
+#include "argList.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "distributedTriSurfaceMesh.H"
+#include "mapDistribute.H"
+#include "triSurfaceFields.H"
+#include "Pair.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Print on master all the per-processor surface stats.
+void writeProcStats
+(
+    const triSurface& s,
+    const List<List<treeBoundBox> >& meshBb
+)
+{
+    // Determine surface bounding boxes, faces, points
+    List<treeBoundBox> surfBb(Pstream::nProcs());
+    {
+        surfBb[Pstream::myProcNo()] = boundBox(s.points(), false);
+        Pstream::gatherList(surfBb);
+        Pstream::scatterList(surfBb);
+    }
+
+    labelList nPoints(Pstream::nProcs());
+    nPoints[Pstream::myProcNo()] = s.points().size();
+    Pstream::gatherList(nPoints);
+    Pstream::scatterList(nPoints);
+
+    labelList nFaces(Pstream::nProcs());
+    nFaces[Pstream::myProcNo()] = s.size();
+    Pstream::gatherList(nFaces);
+    Pstream::scatterList(nFaces);
+
+    forAll(surfBb, procI)
+    {
+        const List<treeBoundBox>& bbs = meshBb[procI];
+
+        Info<< "processor" << procI << endl
+            << "\tMesh bounds          : " << bbs[0] << nl;
+        for (label i = 1; i < bbs.size(); i++)
+        {
+            Info<< "\t                       " << bbs[i]<< nl;
+        }
+        Info<< "\tSurface bounding box : " << surfBb[procI] << nl
+            << "\tTriangles            : " << nFaces[procI] << nl
+            << "\tVertices             : " << nPoints[procI]
+            << endl;
+    }
+    Info<< endl;
+}
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::validArgs.append("triSurfaceMesh");
+    argList::validArgs.append("distributionType");
+
+    argList::validOptions.insert("keepNonMapped", "");
+#   include "setRootCase.H"
+#   include "createTime.H"
+    runTime.functionObjects().off();
+
+    fileName surfFileName(args.additionalArgs()[0]);
+    Info<< "Reading surface from " << surfFileName << nl << endl;
+
+    const word distType(args.additionalArgs()[1]);
+
+    Info<< "Using distribution method "
+        << distributedTriSurfaceMesh::distributionTypeNames_[distType]
+        << " " << distType << nl << endl;
+
+    bool keepNonMapped = args.options().found("keepNonMapped");
+
+    if (keepNonMapped)
+    {
+        Info<< "Preserving surface outside of mesh bounds." << nl << endl;
+    }
+    else
+    {
+        Info<< "Removing surface outside of mesh bounds." << nl << endl;
+    }
+
+
+    if (!Pstream::parRun())
+    {
+        FatalErrorIn(args.executable())
+            << "Please run this program on the decomposed case."
+            << " It will read surface " << surfFileName
+            << " and decompose it such that it overlaps the mesh bounding box."
+            << exit(FatalError);
+    }
+
+
+#   include "createPolyMesh.H"
+
+    Random rndGen(653213);
+
+    // Determine mesh bounding boxes:
+    List<List<treeBoundBox> > meshBb(Pstream::nProcs());
+    {
+        meshBb[Pstream::myProcNo()] = List<treeBoundBox>
+        (
+            1,
+            treeBoundBox
+            (
+                boundBox(mesh.points(), false)
+            ).extend(rndGen, 1E-3)
+        );
+        Pstream::gatherList(meshBb);
+        Pstream::scatterList(meshBb);
+    }
+
+    IOobject io
+    (
+        surfFileName,         // name
+        //runTime.findInstance("triSurface", surfFileName),   // instance
+        runTime.constant(),   // instance
+        "triSurface",         // local
+        runTime,              // registry
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    const fileName actualPath(io.filePath());
+    fileName localPath(actualPath);
+    localPath.replace(runTime.rootPath() + '/', "");
+
+    if (actualPath == io.objectPath())
+    {
+        Info<< "Loading local (decomposed) surface " << localPath << nl <<endl;
+    }
+    else
+    {
+        Info<< "Loading undecomposed surface " << localPath << nl << endl;
+    }
+
+
+    // Create dummy dictionary for bounding boxes if does not exist.
+    if (!isFile(actualPath / "Dict"))
+    {
+        dictionary dict;
+        dict.add("bounds", meshBb[Pstream::myProcNo()]);
+        dict.add("distributionType", distType);
+        dict.add("mergeDistance", SMALL);
+
+        IOdictionary ioDict
+        (
+            IOobject
+            (
+                io.name() + "Dict",
+                io.instance(),
+                io.local(),
+                io.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            dict
+        );
+
+        Info<< "Writing dummy bounds dictionary to " << ioDict.name()
+            << nl << endl;
+
+        ioDict.regIOobject::writeObject
+        (
+            IOstream::ASCII,
+            IOstream::currentVersion,
+            ioDict.time().writeCompression()        
+        );
+    }
+
+
+    // Load surface
+    distributedTriSurfaceMesh surfMesh(io);
+    Info<< "Loaded surface" << nl << endl;
+
+
+    // Generate a test field
+    {
+        const triSurface& s = static_cast<const triSurface&>(surfMesh);
+
+        autoPtr<triSurfaceVectorField> fcPtr
+        (
+            new triSurfaceVectorField
+            (
+                IOobject
+                (
+                    surfMesh.searchableSurface::name(),     // name
+                    surfMesh.searchableSurface::instance(), // instance
+                    surfMesh.searchableSurface::local(),    // local
+                    surfMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                surfMesh,
+                dimLength
+            )
+        );
+        triSurfaceVectorField& fc = fcPtr();
+
+        forAll(fc, triI)
+        {
+            fc[triI] = s[triI].centre(s.points());
+        }
+
+        // Steal pointer and store object on surfMesh
+        fcPtr.ptr()->store();
+    }
+
+
+    // Write per-processor stats
+    Info<< "Before redistribution:" << endl;
+    writeProcStats(surfMesh, meshBb);
+
+
+    // Do redistribution
+    Info<< "Redistributing surface" << nl << endl;
+    autoPtr<mapDistribute> faceMap;
+    autoPtr<mapDistribute> pointMap;
+    surfMesh.distribute
+    (
+        meshBb[Pstream::myProcNo()],
+        keepNonMapped,
+        faceMap,
+        pointMap
+    );
+    faceMap.clear();
+    pointMap.clear();
+
+    Info<< endl;
+
+
+    // Write per-processor stats
+    Info<< "After redistribution:" << endl;
+    writeProcStats(surfMesh, meshBb);
+
+
+    Info<< "Writing surface." << nl << endl;
+    surfMesh.searchableSurface::write();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

etc/apps/paraview3/bashrc

@@ -46,8 +46,8 @@ do
     fi
 done
 
-paraviewMajor=paraview-3.5
+paraviewMajor=paraview-3.6
-export ParaView_VERSION=3.5-cvs
+export ParaView_VERSION=3.6
 
 export ParaView_INST_DIR=$WM_THIRD_PARTY_DIR/paraview-$ParaView_VERSION
 export ParaView_DIR=$ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER

etc/apps/paraview3/cshrc

@@ -44,8 +44,8 @@ foreach cmake ( cmake-2.6.4 cmake-2.6.2 cmake-2.4.6 )
     endif
 end
 
-set paraviewMajor=paraview-3.5
+set paraviewMajor=paraview-3.6
-setenv ParaView_VERSION 3.5-cvs
+setenv ParaView_VERSION 3.6
 
 setenv ParaView_INST_DIR $WM_THIRD_PARTY_DIR/paraview-$ParaView_VERSION
 setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H

@@ -84,6 +84,9 @@ public:
         //- Start of procI+1 data
         inline const labelList& offsets() const;
 
+        //- my local size
+        inline label localSize() const;
+
         //- Global sum of localSizes
         inline label size() const;
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H

@@ -34,6 +34,17 @@ inline const Foam::labelList& Foam::globalIndex::offsets() const
 }
 
 
+inline Foam::label Foam::globalIndex::localSize() const
+{
+    return
+    (
+        Pstream::myProcNo() == 0
+      ? offsets_[Pstream::myProcNo()]
+      : offsets_[Pstream::myProcNo()] - offsets_[Pstream::myProcNo()-1]
+    );
+}
+
+
 inline Foam::label Foam::globalIndex::size() const
 {
     return offsets_[Pstream::nProcs()-1];

src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C

@@ -75,6 +75,24 @@ extern "C"
 }
 
 
+// Hack: scotch generates floating point errors so need to switch of error
+//       trapping!
+#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
+#    define LINUX
+#endif
+
+#if defined(LINUX) && defined(__GNUC__)
+#    define LINUX_GNUC
+#endif
+
+#ifdef LINUX_GNUC
+#   ifndef __USE_GNU
+#       define __USE_GNU
+#   endif
+#   include <fenv.h>
+#endif
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -113,13 +131,30 @@ Foam::label Foam::scotchDecomp::decompose
 {
     // Strategy
     // ~~~~~~~~
+
     // Default.
     SCOTCH_Strat stradat;
     check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
-    //SCOTCH_stratGraphMap(&stradat, &argv[i][2]);
+
-    //fprintf(stdout, "S\tStrat=");
+    if (decompositionDict_.found("scotchCoeffs"))
-    //SCOTCH_stratSave(&stradat, stdout);
+    {
-    //fprintf(stdout, "\n");
+        const dictionary& scotchCoeffs =
+            decompositionDict_.subDict("scotchCoeffs");
+
+
+        string strategy;
+        if (scotchCoeffs.readIfPresent("strategy", strategy))
+        {
+            if (debug)
+            {
+                Info<< "scotchDecomp : Using strategy " << strategy << endl;
+            }
+            SCOTCH_stratGraphMap(&stradat, strategy.c_str());
+            //fprintf(stdout, "S\tStrat=");
+            //SCOTCH_stratSave(&stradat, stdout);
+            //fprintf(stdout, "\n");
+        }
+    }
 
 
     // Graph
@@ -153,37 +188,40 @@ Foam::label Foam::scotchDecomp::decompose
         const dictionary& scotchCoeffs =
             decompositionDict_.subDict("scotchCoeffs");
 
-        Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
+        if (scotchCoeffs.found("writeGraph"))
-
-        if (writeGraph)
         {
-            OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
+            Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
 
-            Info<< "Dumping Scotch graph file to " << str.name() << endl
+            if (writeGraph)
-                << "Use this in combination with gpart." << endl;
-
-            label version = 0;
-            str << version << nl;
-            // Numer of vertices
-            str << xadj.size()-1 << ' ' << adjncy.size() << nl;
-            // Numbering starts from 0
-            label baseval = 0;
-            // Has weights?
-            label hasEdgeWeights = 0;
-            label hasVertexWeights = 0;
-            label numericflag = 10*hasEdgeWeights+hasVertexWeights;
-            str << baseval << ' ' << numericflag << nl;
-            for (label cellI = 0; cellI < xadj.size()-1; cellI++)
             {
-                label start = xadj[cellI];
+                OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
-                label end = xadj[cellI+1];
-                str << end-start;
 
-                for (label i = start; i < end; i++)
+                Info<< "Dumping Scotch graph file to " << str.name() << endl
+                    << "Use this in combination with gpart." << endl;
+
+                label version = 0;
+                str << version << nl;
+                // Numer of vertices
+                str << xadj.size()-1 << ' ' << adjncy.size() << nl;
+                // Numbering starts from 0
+                label baseval = 0;
+                // Has weights?
+                label hasEdgeWeights = 0;
+                label hasVertexWeights = 0;
+                label numericflag = 10*hasEdgeWeights+hasVertexWeights;
+                str << baseval << ' ' << numericflag << nl;
+                for (label cellI = 0; cellI < xadj.size()-1; cellI++)
                 {
-                    str << ' ' << adjncy[i];
+                    label start = xadj[cellI];
+                    label end = xadj[cellI+1];
+                    str << end-start;
+
+                    for (label i = start; i < end; i++)
+                    {
+                        str << ' ' << adjncy[i];
+                    }
+                    str << nl;
                 }
-                str << nl;
             }
         }
     }
@@ -195,12 +233,36 @@ Foam::label Foam::scotchDecomp::decompose
 
     SCOTCH_Arch archdat;
     check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
-    check
+
-    (
+    List<label> processorWeights;
-        // SCOTCH_archCmpltw for weighted.
+    if (decompositionDict_.found("scotchCoeffs"))
-        SCOTCH_archCmplt(&archdat, nProcessors_),
+    {
-        "SCOTCH_archCmplt"
+        const dictionary& scotchCoeffs =
-    );
+            decompositionDict_.subDict("scotchCoeffs");
+
+        scotchCoeffs.readIfPresent("processorWeights", processorWeights);
+    }
+    if (processorWeights.size())
+    {
+        if (debug)
+        {
+            Info<< "scotchDecomp : Using procesor weights " << processorWeights
+                << endl;
+        }
+        check
+        (
+            SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
+            "SCOTCH_archCmpltw"
+        );
+    }
+    else
+    {
+        check
+        (
+            SCOTCH_archCmplt(&archdat, nProcessors_),
+            "SCOTCH_archCmplt"
+        );
+    }
 
 
     //SCOTCH_Mapping mapdat;
@@ -209,6 +271,16 @@ Foam::label Foam::scotchDecomp::decompose
     //SCOTCH_graphMapExit(&grafdat, &mapdat);
 
 
+    // Hack:switch off fpu error trapping
+#   ifdef LINUX_GNUC
+    int oldExcepts = fedisableexcept
+    (
+        FE_DIVBYZERO
+      | FE_INVALID
+      | FE_OVERFLOW
+    );
+#   endif
+
     finalDecomp.setSize(xadj.size()-1);
     finalDecomp = 0;
     check
@@ -223,6 +295,11 @@ Foam::label Foam::scotchDecomp::decompose
         "SCOTCH_graphMap"
     );
 
+#   ifdef LINUX_GNUC
+    feenableexcept(oldExcepts);
+#   endif
+
+
 
     //finalDecomp.setSize(xadj.size()-1);
     //check

src/finiteVolume/Make/files

@@ -39,26 +39,40 @@ fvMeshMapper = fvMesh/fvMeshMapper
 $(fvMeshMapper)/fvPatchMapper.C
 $(fvMeshMapper)/fvSurfaceMapper.C
 
+
 extendedStencil = fvMesh/extendedStencil
-$(extendedStencil)/extendedStencil.C
-$(extendedStencil)/extendedUpwindStencil.C
-$(extendedStencil)/extendedCentredStencil.C
 
-$(extendedStencil)/faceStencil/faceStencil.C
+cellToCell = $(extendedStencil)/cellToCell
-$(extendedStencil)/faceStencil/faceEdgeCellStencil.C
+$(cellToCell)/fullStencils/cellToCellStencil.C
-$(extendedStencil)/faceStencil/cellFaceCellStencil.C
+$(cellToCell)/fullStencils/CFCCellToCellStencil.C
-$(extendedStencil)/faceStencil/cellPointCellStencil.C
+$(cellToCell)/fullStencils/CPCCellToCellStencil.C
-$(extendedStencil)/faceStencil/cellEdgeCellStencil.C
+$(cellToCell)/fullStencils/CECCellToCellStencil.C
 
-$(extendedStencil)/extendedStencilMeshObjects/centredCECStencilObject.C
+cellToFace = $(extendedStencil)/cellToFace
-$(extendedStencil)/extendedStencilMeshObjects/centredCFCStencilObject.C
+$(cellToFace)/fullStencils/cellToFaceStencil.C
-$(extendedStencil)/extendedStencilMeshObjects/centredCPCStencilObject.C
+$(cellToFace)/fullStencils/CFCCellToFaceStencil.C
-$(extendedStencil)/extendedStencilMeshObjects/centredFECStencilObject.C
+$(cellToFace)/fullStencils/CECCellToFaceStencil.C
+$(cellToFace)/fullStencils/CPCCellToFaceStencil.C
+$(cellToFace)/fullStencils/FECCellToFaceStencil.C
+$(cellToFace)/extendedCellToFaceStencil.C
+$(cellToFace)/extendedCentredCellToFaceStencil.C
+$(cellToFace)/extendedUpwindCellToFaceStencil.C
+$(cellToFace)/MeshObjects/centredCECCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/centredCFCCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/centredCPCCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/centredFECCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/upwindCECCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/upwindCFCCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/upwindCPCCellToFaceStencilObject.C
+$(cellToFace)/MeshObjects/upwindFECCellToFaceStencilObject.C
+
+faceToCell = $(extendedStencil)/faceToCell
+$(faceToCell)/fullStencils/faceToCellStencil.C
+$(faceToCell)/fullStencils/CFCFaceToCellStencil.C
+$(faceToCell)/extendedFaceToCellStencil.C
+$(faceToCell)/extendedCentredFaceToCellStencil.C
+$(faceToCell)/MeshObjects/centredCFCFaceToCellStencilObject.C
 
-$(extendedStencil)/extendedStencilMeshObjects/upwindCECStencilObject.C
-$(extendedStencil)/extendedStencilMeshObjects/upwindCFCStencilObject.C
-$(extendedStencil)/extendedStencilMeshObjects/upwindCPCStencilObject.C
-$(extendedStencil)/extendedStencilMeshObjects/upwindFECStencilObject.C
 
 fvPatchFields = fields/fvPatchFields
 $(fvPatchFields)/fvPatchField/fvPatchFields.C

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C

@@ -108,7 +108,7 @@ void Foam::setRefCell
                  "    bool\n"
                  ")",
                 dict
-            )   << "Unable to set reference cell for field" << field.name()
+            )   << "Unable to set reference cell for field " << field.name()
                 << nl
                 << "    Please supply either " << refCellName
                 << " or " << refPointName << nl << exit(FatalIOError);

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C

@@ -113,7 +113,17 @@ void fixedFluxBuoyantPressureFvPatchScalarField::updateCoeffs()
     const fvPatchField<scalar>& rho =
         patch().lookupPatchField<volScalarField, scalar>("rho");
 
-    gradient() = -rho.snGrad()*(g.value() & patch().Cf());
+    // If the variable name is "pd" assume it is p - rho*g.h
+    // and set the gradient appropriately.
+    // Otherwise assume the variable is the static pressure.
+    if (dimensionedInternalField().name() == "pd")
+    {
+        gradient() = -rho.snGrad()*(g.value() & patch().Cf());
+    }
+    else
+    {
+        gradient() = rho*(g.value() & patch().nf());
+    }
 
     fixedGradientFvPatchScalarField::updateCoeffs();
 }

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H

@@ -26,7 +26,10 @@ Class
     Foam::fixedFluxBuoyantPressureFvPatchScalarField
 
 Description
-    Foam::fixedFluxBuoyantPressureFvPatchScalarField
+    Set the pressure gradient boundary condition appropriately for buoyant flow.
+
+    If the variable name is "pd" assume it is p - rho*g.h and set the gradient
+    appropriately.  Otherwise assume the variable is the static pressure.
 
 SourceFiles
     fixedFluxBuoyantPressureFvPatchScalarField.C

src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H

@@ -39,7 +39,7 @@ SourceFiles
 
 #include "snGradScheme.H"
 #include "quadraticFitSnGradData.H"
-#include "extendedStencil.H"
+#include "extendedCellToFaceStencil.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -133,7 +133,7 @@ public:
                 centralWeight_
             );
 
-            const extendedStencil& stencil = qfd.stencil();
+            const extendedCellToFaceStencil& stencil = qfd.stencil();
             const List<scalarList>& f = qfd.fit();
 
             tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > sft

src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H

@@ -38,7 +38,7 @@ SourceFiles
 
 #include "MeshObject.H"
 #include "fvMesh.H"
-#include "extendedStencil.H"
+#include "extendedCellToFaceStencil.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,7 +62,7 @@ class quadraticFitSnGradData
         const label minSize_;
 
         //- Extended stencil addressing
-        extendedStencil stencil_;
+        extendedCellToFaceStencil stencil_;
 
         //- For each cell in the mesh store the values which multiply the
         //  values of the stencil to obtain the gradient for each direction
@@ -107,7 +107,7 @@ public:
     // Member functions
 
         //- Return reference to the stencil
-        const extendedStencil& stencil() const
+        const extendedCellToFaceStencil& stencil() const
         {
             return stencil_;
         }

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -480,7 +480,7 @@ void Foam::fvMatrix<Type>::setReference
     const bool forceReference
 )
 {
-    if (celli >= 0 && (psi_.needReference() || forceReference))
+    if ((forceReference || psi_.needReference()) && celli >= 0)
     {
         source()[celli] += diag()[celli]*value;
         diag()[celli] += diag()[celli];

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CECCellToCellStencil.C

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cellEdgeCellStencil.H"
+#include "CECCellToCellStencil.H"
 #include "syncTools.H"
 //#include "meshTools.H"
 //#include "OFstream.H"
@@ -33,7 +33,7 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Calculates per edge the neighbour data (= edgeCells)
-void Foam::cellEdgeCellStencil::calcEdgeBoundaryData
+void Foam::CECCellToCellStencil::calcEdgeBoundaryData
 (
     const boolList& isValidBFace,
     const labelList& boundaryEdges,
@@ -72,7 +72,7 @@ void Foam::cellEdgeCellStencil::calcEdgeBoundaryData
 
 // Calculates per cell the neighbour data (= cell or boundary in global
 // numbering). First element is always cell itself!
-void Foam::cellEdgeCellStencil::calcCellStencil
+void Foam::CECCellToCellStencil::calcCellStencil
 (
     labelListList& globalCellCells
 ) const
@@ -189,20 +189,12 @@ void Foam::cellEdgeCellStencil::calcCellStencil
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::cellEdgeCellStencil::cellEdgeCellStencil(const polyMesh& mesh)
+Foam::CECCellToCellStencil::CECCellToCellStencil(const polyMesh& mesh)
 :
-    faceStencil(mesh)
+    cellToCellStencil(mesh)
 {
     // Calculate per cell the (edge) connected cells (in global numbering)
-    labelListList globalCellCells;
+    calcCellStencil(*this);
-    calcCellStencil(globalCellCells);
-
-    // Add stencils of neighbouring cells to create faceStencil
-    labelListList faceStencil;
-    calcFaceStencil(globalCellCells, faceStencil);
-
-    // Transfer to *this
-    transfer(faceStencil);
 }
 
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CECCellToCellStencil.H

@@ -23,19 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::cellEdgeCellStencil
+    Foam::CECCellToCellStencil
 
 Description
 
 SourceFiles
-    cellEdgeCellStencil.C
+    CECCellToCellStencil.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cellEdgeCellStencil_H
+#ifndef CECCellToCellStencil_H
-#define cellEdgeCellStencil_H
+#define CECCellToCellStencil_H
 
-#include "faceStencil.H"
+#include "cellToCellStencil.H"
 #include "boolList.H"
 #include "HashSet.H"
 #include "Map.H"
@@ -47,12 +47,12 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class cellEdgeCellStencil Declaration
+                           Class CECCellToCellStencil Declaration
 \*---------------------------------------------------------------------------*/
 
-class cellEdgeCellStencil
+class CECCellToCellStencil
 :
-    public faceStencil
+    public cellToCellStencil
 {
     // Private Member Functions
 
@@ -68,10 +68,10 @@ class cellEdgeCellStencil
 
 
         //- Disallow default bitwise copy construct
-        cellEdgeCellStencil(const cellEdgeCellStencil&);
+        CECCellToCellStencil(const CECCellToCellStencil&);
 
         //- Disallow default bitwise assignment
-        void operator=(const cellEdgeCellStencil&);
+        void operator=(const CECCellToCellStencil&);
 
 
 public:
@@ -79,7 +79,7 @@ public:
     // Constructors
 
         //- Construct from all cells and boundary faces
-        explicit cellEdgeCellStencil(const polyMesh&);
+        explicit CECCellToCellStencil(const polyMesh&);
 
 };
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CFCCellToCellStencil.C

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cellFaceCellStencil.H"
+#include "CFCCellToCellStencil.H"
 #include "syncTools.H"
 #include "SortableList.H"
 #include "emptyPolyPatch.H"
@@ -32,7 +32,7 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Calculates per face the neighbour data (= cell or boundary face)
-void Foam::cellFaceCellStencil::calcFaceBoundaryData
+void Foam::CFCCellToCellStencil::calcFaceBoundaryData
 (
     labelList& neiGlobal
 ) const
@@ -85,7 +85,7 @@ void Foam::cellFaceCellStencil::calcFaceBoundaryData
 
 // Calculates per cell the neighbour data (= cell or boundary in global
 // numbering). First element is always cell itself!
-void Foam::cellFaceCellStencil::calcCellStencil(labelListList& globalCellCells)
+void Foam::CFCCellToCellStencil::calcCellStencil(labelListList& globalCellCells)
  const
 {
     const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
@@ -147,20 +147,12 @@ void Foam::cellFaceCellStencil::calcCellStencil(labelListList& globalCellCells)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::cellFaceCellStencil::cellFaceCellStencil(const polyMesh& mesh)
+Foam::CFCCellToCellStencil::CFCCellToCellStencil(const polyMesh& mesh)
 :
-    faceStencil(mesh)
+    cellToCellStencil(mesh)
 {
     // Calculate per cell the (face) connected cells (in global numbering)
-    labelListList globalCellCells;
+    calcCellStencil(*this);
-    calcCellStencil(globalCellCells);
-
-    // Add stencils of neighbouring cells to create faceStencil
-    labelListList faceStencil;
-    calcFaceStencil(globalCellCells, faceStencil);
-
-    // Transfer to *this
-    transfer(faceStencil);
 }
 
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CFCCellToCellStencil.H

@@ -23,19 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::cellFaceCellStencil
+    Foam::CFCCellToCellStencil
 
 Description
 
 SourceFiles
-    cellFaceCellStencil.C
+    CFCCellToCellStencil.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cellFaceCellStencil_H
+#ifndef CFCCellToCellStencil_H
-#define cellFaceCellStencil_H
+#define CFCCellToCellStencil_H
 
-#include "faceStencil.H"
+#include "cellToCellStencil.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,12 +43,12 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class cellFaceCellStencil Declaration
+                           Class CFCCellToCellStencil Declaration
 \*---------------------------------------------------------------------------*/
 
-class cellFaceCellStencil
+class CFCCellToCellStencil
 :
-    public faceStencil
+    public cellToCellStencil
 {
     // Private Member Functions
 
@@ -57,17 +57,17 @@ class cellFaceCellStencil
         void calcCellStencil(labelListList& globalCellCells) const;
 
         //- Disallow default bitwise copy construct
-        cellFaceCellStencil(const cellFaceCellStencil&);
+        CFCCellToCellStencil(const CFCCellToCellStencil&);
 
         //- Disallow default bitwise assignment
-        void operator=(const cellFaceCellStencil&);
+        void operator=(const CFCCellToCellStencil&);
 
 public:
 
     // Constructors
 
         //- Construct from mesh
-        explicit cellFaceCellStencil(const polyMesh& mesh);
+        explicit CFCCellToCellStencil(const polyMesh& mesh);
 };
 
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CPCCellToCellStencil.C

@@ -24,13 +24,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cellPointCellStencil.H"
+#include "CPCCellToCellStencil.H"
 #include "syncTools.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Calculates per point the neighbour data (= pointCells)
-void Foam::cellPointCellStencil::calcPointBoundaryData
+void Foam::CPCCellToCellStencil::calcPointBoundaryData
 (
     const boolList& isValidBFace,
     const labelList& boundaryPoints,
@@ -69,7 +69,7 @@ void Foam::cellPointCellStencil::calcPointBoundaryData
 
 // Calculates per cell the neighbour data (= cell or boundary in global
 // numbering). First element is always cell itself!
-void Foam::cellPointCellStencil::calcCellStencil
+void Foam::CPCCellToCellStencil::calcCellStencil
 (
     labelListList& globalCellCells
 ) const
@@ -154,20 +154,13 @@ void Foam::cellPointCellStencil::calcCellStencil
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::cellPointCellStencil::cellPointCellStencil(const polyMesh& mesh)
+Foam::CPCCellToCellStencil::CPCCellToCellStencil(const polyMesh& mesh)
 :
-    faceStencil(mesh)
+    cellToCellStencil(mesh)
 {
     // Calculate per cell the (point) connected cells (in global numbering)
     labelListList globalCellCells;
-    calcCellStencil(globalCellCells);
+    calcCellStencil(*this);
-
-    // Add stencils of neighbouring cells to create faceStencil
-    labelListList faceStencil;
-    calcFaceStencil(globalCellCells, faceStencil);
-
-    // Transfer to *this
-    transfer(faceStencil);
 }
 
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/CPCCellToCellStencil.H

@@ -23,19 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::cellPointCellStencil
+    Foam::CPCCellToCellStencil
 
 Description
 
 SourceFiles
-    cellPointCellStencil.C
+    CPCCellToCellStencil.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cellPointCellStencil_H
+#ifndef CPCCellToCellStencil_H
-#define cellPointCellStencil_H
+#define CPCCellToCellStencil_H
 
-#include "faceStencil.H"
+#include "cellToCellStencil.H"
 #include "boolList.H"
 #include "HashSet.H"
 #include "Map.H"
@@ -46,12 +46,12 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class cellPointCellStencil Declaration
+                           Class CPCCellToCellStencil Declaration
 \*---------------------------------------------------------------------------*/
 
-class cellPointCellStencil
+class CPCCellToCellStencil
 :
-    public faceStencil
+    public cellToCellStencil
 {
     // Private Member Functions
 
@@ -67,10 +67,10 @@ class cellPointCellStencil
 
 
         //- Disallow default bitwise copy construct
-        cellPointCellStencil(const cellPointCellStencil&);
+        CPCCellToCellStencil(const CPCCellToCellStencil&);
 
         //- Disallow default bitwise assignment
-        void operator=(const cellPointCellStencil&);
+        void operator=(const CPCCellToCellStencil&);
 
 
 public:
@@ -78,7 +78,7 @@ public:
     // Constructors
 
         //- Construct from all cells and boundary faces
-        explicit cellPointCellStencil(const polyMesh&);
+        explicit CPCCellToCellStencil(const polyMesh&);
 
 };
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.C

@@ -0,0 +1,350 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "cellToCellStencil.H"
+#include "syncTools.H"
+#include "SortableList.H"
+#include "emptyPolyPatch.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Merge two list and guarantee global0,global1 are first.
+void Foam::cellToCellStencil::merge
+(
+    const label global0,
+    const label global1,
+    const labelList& listA,
+    labelList& listB
+)
+{
+    sort(listB);
+
+    // See if global0, global1 already present in listB
+    label nGlobalInsert = 0;
+
+    if (global0 != -1)
+    {
+        label index0 = findSortedIndex(listB, global0);
+        if (index0 == -1)
+        {
+            nGlobalInsert++;
+        }
+    }
+
+    if (global1 != -1)
+    {
+        label index1 = findSortedIndex(listB, global1);
+        if (index1 == -1)
+        {
+            nGlobalInsert++;
+        }
+    }
+
+
+    // For all in listA see if they are present
+    label nInsert = 0;
+
+    forAll(listA, i)
+    {
+        label elem = listA[i];
+
+        if (elem != global0 && elem != global1)
+        {
+            if (findSortedIndex(listB, elem) == -1)
+            {
+                nInsert++;
+            }
+        }
+    }
+
+    // Extend B with nInsert and whether global0,global1 need to be inserted.
+    labelList result(listB.size() + nGlobalInsert + nInsert);
+
+    label resultI = 0;
+
+    // Insert global0,1 first
+    if (global0 != -1)
+    {
+        result[resultI++] = global0;
+    }
+    if (global1 != -1)
+    {
+        result[resultI++] = global1;
+    }
+
+
+    // Insert listB
+    forAll(listB, i)
+    {
+        label elem = listB[i];
+
+        if (elem != global0 && elem != global1)
+        {
+            result[resultI++] = elem;
+        }
+    }
+
+
+    // Insert listA
+    forAll(listA, i)
+    {
+        label elem = listA[i];
+
+        if (elem != global0 && elem != global1)
+        {
+            if (findSortedIndex(listB, elem) == -1)
+            {
+                result[resultI++] = elem;
+            }
+        }
+    }
+
+    if (resultI != result.size())
+    {
+        FatalErrorIn("cellToCellStencil::merge(..)")
+            << "problem" << abort(FatalError);
+    }
+
+    listB.transfer(result);
+}
+
+
+// Merge two list and guarantee globalI is first.
+void Foam::cellToCellStencil::merge
+(
+    const label globalI,
+    const labelList& pGlobals,
+    labelList& cCells
+)
+{
+    labelHashSet set;
+    forAll(cCells, i)
+    {
+        if (cCells[i] != globalI)
+        {
+            set.insert(cCells[i]);
+        }
+    }
+
+    forAll(pGlobals, i)
+    {
+        if (pGlobals[i] != globalI)
+        {
+            set.insert(pGlobals[i]);
+        }
+    }
+
+    cCells.setSize(set.size()+1);
+    label n = 0;
+    cCells[n++] = globalI;
+
+    forAllConstIter(labelHashSet, set, iter)
+    {
+        cCells[n++] = iter.key();
+    }
+}
+
+
+void Foam::cellToCellStencil::validBoundaryFaces(boolList& isValidBFace) const
+{
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+    isValidBFace.setSize(mesh().nFaces()-mesh().nInternalFaces(), true);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled() || isA<emptyPolyPatch>(pp))
+        {
+            label bFaceI = pp.start()-mesh().nInternalFaces();
+            forAll(pp, i)
+            {
+                isValidBFace[bFaceI++] = false;
+            }
+        }
+    }
+}
+
+
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::cellToCellStencil::allCoupledFacesPatch() const
+{
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+    label nCoupled = 0;
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            nCoupled += pp.size();
+        }
+    }
+    labelList coupledFaces(nCoupled);
+    nCoupled = 0;
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                coupledFaces[nCoupled++] = faceI++;
+            }
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>
+            (
+                mesh().faces(),
+                coupledFaces
+            ),
+            mesh().points()
+        )
+    );
+}
+
+
+void Foam::cellToCellStencil::unionEqOp::operator()
+(
+    labelList& x,
+    const labelList& y
+) const
+{
+    if (y.size())
+    {
+        if (x.empty())
+        {
+            x = y;
+        }
+        else
+        {
+            labelHashSet set(x);
+            forAll(y, i)
+            {
+                set.insert(y[i]);
+            }
+            x = set.toc();
+        }
+    }
+}
+
+
+void Foam::cellToCellStencil::insertFaceCells
+(
+    const label exclude0,
+    const label exclude1,
+    const boolList& isValidBFace,
+    const labelList& faceLabels,
+    labelHashSet& globals
+) const
+{
+    const labelList& own = mesh().faceOwner();
+    const labelList& nei = mesh().faceNeighbour();
+
+    forAll(faceLabels, i)
+    {
+        label faceI = faceLabels[i];
+
+        label globalOwn = globalNumbering().toGlobal(own[faceI]);
+        if (globalOwn != exclude0 && globalOwn != exclude1)
+        {
+            globals.insert(globalOwn);
+        }
+
+        if (mesh().isInternalFace(faceI))
+        {
+            label globalNei = globalNumbering().toGlobal(nei[faceI]);
+            if (globalNei != exclude0 && globalNei != exclude1)
+            {
+                globals.insert(globalNei);
+            }
+        }
+        else
+        {
+            label bFaceI = faceI-mesh().nInternalFaces();
+
+            if (isValidBFace[bFaceI])
+            {
+                label globalI = globalNumbering().toGlobal
+                (
+                    mesh().nCells()
+                  + bFaceI
+                );
+
+                if (globalI != exclude0 && globalI != exclude1)
+                {
+                    globals.insert(globalI);
+                }
+            }
+        }
+    }
+}
+
+
+Foam::labelList Foam::cellToCellStencil::calcFaceCells
+(
+    const boolList& isValidBFace,
+    const labelList& faceLabels,
+    labelHashSet& globals
+) const
+{
+    globals.clear();
+
+    insertFaceCells
+    (
+        -1,
+        -1,
+        isValidBFace,
+        faceLabels,
+        globals
+    );
+
+    return globals.toc();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cellToCellStencil::cellToCellStencil(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    globalNumbering_(mesh_.nCells()+mesh_.nFaces()-mesh_.nInternalFaces())
+{}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::cellToCellStencil
+
+Description
+    baseclass for extended cell centred addressing. Contains per cell a
+    list of neighbouring cells and/or boundaryfaces in global addressing.
+
+SourceFiles
+    cellToCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellToCellStencil_H
+#define cellToCellStencil_H
+
+#include "globalIndex.H"
+#include "boolList.H"
+#include "HashSet.H"
+#include "indirectPrimitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class cellToCellStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellToCellStencil
+:
+    public labelListList
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        //- Global numbering for cells and boundary faces
+        const globalIndex globalNumbering_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const cellToCellStencil&);
+
+        
+protected:
+
+        //- Merge two lists.
+        static void merge
+        (
+            const label,
+            const label,
+            const labelList&,
+            labelList&
+        );
+
+        //- Merge two lists.
+        static void merge(const label, const labelList&, labelList&);
+
+        //- Valid boundary faces (not empty and not coupled)
+        void validBoundaryFaces(boolList& isValidBFace) const;
+
+        //- Return patch of all coupled faces.
+        autoPtr<indirectPrimitivePatch> allCoupledFacesPatch() const;
+
+        //- Combine operator for labelLists
+        class unionEqOp
+        {
+            public:
+            void operator()( labelList& x, const labelList& y ) const;
+        };
+
+        //- Collect cell neighbours of faces in global numbering
+        void insertFaceCells
+        (
+            const label exclude0,
+            const label exclude1,
+            const boolList& nonEmptyFace,
+            const labelList& faceLabels,
+            labelHashSet& globals
+        ) const;
+
+        //- Collect cell neighbours of faces in global numbering
+        labelList calcFaceCells
+        (
+            const boolList& nonEmptyFace,
+            const labelList& faceLabels,
+            labelHashSet& globals
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh
+        explicit cellToCellStencil(const polyMesh&);
+
+
+    // Member Functions
+
+        const polyMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Global numbering for cells and boundary faces
+        const globalIndex& globalNumbering() const
+        {
+            return globalNumbering_;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "centredCECCellToFaceStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(centredCECCellToFaceStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::centredCECStencilObject
+    Foam::centredCECCellToFaceStencilObject
 
 Description
 
@@ -31,11 +31,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef centredCECStencilObject_H
+#ifndef centredCECCellToFaceStencilObject_H
-#define centredCECStencilObject_H
+#define centredCECCellToFaceStencilObject_H
 
-#include "extendedCentredStencil.H"
+#include "extendedCentredCellToFaceStencil.H"
-#include "cellEdgeCellStencil.H"
+#include "CECCellToFaceStencil.H"
 #include "MeshObject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,35 +44,35 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class centredCECStencilObject Declaration
+                           Class centredCECCellToFaceStencilObject Declaration
 \*---------------------------------------------------------------------------*/
 
-class centredCECStencilObject
+class centredCECCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, centredCECStencilObject>,
+    public MeshObject<fvMesh, centredCECCellToFaceStencilObject>,
-    public extendedCentredStencil
+    public extendedCentredCellToFaceStencil
 {
 
 public:
 
-    TypeName("centredCECStencil");
+    TypeName("centredCECCellToFaceStencil");
 
     // Constructors
 
         //- Construct from uncompacted face stencil
-        explicit centredCECStencilObject
+        explicit centredCECCellToFaceStencilObject
         (
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCECStencilObject>(mesh),
+            MeshObject<fvMesh, centredCECCellToFaceStencilObject>(mesh),
-            extendedCentredStencil(cellEdgeCellStencil(mesh))
+            extendedCentredCellToFaceStencil(CECCellToFaceStencil(mesh))
         {}
 
 
     // Destructor
 
-        virtual ~centredCECStencilObject()
+        virtual ~centredCECCellToFaceStencilObject()
         {}
 };
 

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "centredCFCCellToFaceStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(centredCFCCellToFaceStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::centredCFCStencilObject
+    Foam::centredCFCCellToFaceStencilObject
 
 Description
 
@@ -31,11 +31,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef centredCFCStencilObject_H
+#ifndef centredCFCCellToFaceStencilObject_H
-#define centredCFCStencilObject_H
+#define centredCFCCellToFaceStencilObject_H
 
-#include "extendedCentredStencil.H"
+#include "extendedCentredCellToFaceStencil.H"
-#include "cellFaceCellStencil.H"
+#include "CFCCellToFaceStencil.H"
 #include "MeshObject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,34 +44,34 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class centredCFCStencilObject Declaration
+                           Class centredCFCCellToFaceStencilObject Declaration
 \*---------------------------------------------------------------------------*/
 
-class centredCFCStencilObject
+class centredCFCCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, centredCFCStencilObject>,
+    public MeshObject<fvMesh, centredCFCCellToFaceStencilObject>,
-    public extendedCentredStencil
+    public extendedCentredCellToFaceStencil
 {
 
 public:
 
-    TypeName("centredCFCStencil");
+    TypeName("centredCFCCellToFaceStencil");
 
     // Constructors
 
         //- Construct from uncompacted face stencil
-        explicit centredCFCStencilObject
+        explicit centredCFCCellToFaceStencilObject
         (
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCFCStencilObject>(mesh),
+            MeshObject<fvMesh, centredCFCCellToFaceStencilObject>(mesh),
-            extendedCentredStencil(cellFaceCellStencil(mesh))
+            extendedCentredCellToFaceStencil(CFCCellToFaceStencil(mesh))
         {}
 
 
     //- Destructor
-    virtual ~centredCFCStencilObject()
+    virtual ~centredCFCCellToFaceStencilObject()
     {}
 };