From 3c060421798a05412c4e982e2facee9c2192e61a Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Tue, 13 Apr 2010 18:24:54 +0100
Subject: [PATCH 01/12] ENH: particleForces.  Adding paramagnetic force. 
 Adding d to calcCoupled/calcNonCoupled arguments and call.

---
 .../KinematicParcel/KinematicParcel.C         | 25 +++++-
 .../particleForces/particleForces.C           | 90 +++++++++++++++++--
 .../particleForces/particleForces.H           | 31 ++++++-
 3 files changed, 133 insertions(+), 13 deletions(-)

diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index a381c9d431..e17e7138d6 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -160,10 +160,27 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const scalar utc = td.cloud().drag().utc(Re, d, mu) + ROOTVSMALL;
 
     // Momentum source due to particle forces
-    const vector FCoupled =
-        mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U);
-    const vector FNonCoupled =
-        mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U);
+    const vector FCoupled = mass*td.cloud().forces().calcCoupled
+    (
+        cellI,
+        dt,
+        rhoc_,
+        rho,
+        Uc_,
+        U,
+        d
+    );
+
+    const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled
+    (
+        cellI,
+        dt,
+        rhoc_,
+        rho,
+        Uc_,
+        U,
+        d
+    );
 
 
     // New particle velocity
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.C b/src/lagrangian/intermediate/particleForces/particleForces.C
index 265f9b59ae..938549a734 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.C
+++ b/src/lagrangian/intermediate/particleForces/particleForces.C
@@ -27,6 +27,8 @@ License
 #include "fvMesh.H"
 #include "volFields.H"
 #include "fvcGrad.H"
+#include "mathematicalConstants.H"
+#include "electromagneticConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -41,16 +43,25 @@ Foam::particleForces::particleForces
     dict_(dict.subDict("particleForces")),
     g_(g),
     gradUPtr_(NULL),
+    gradHPtr_(NULL),
     gravity_(dict_.lookup("gravity")),
     virtualMass_(dict_.lookup("virtualMass")),
     Cvm_(0.0),
     pressureGradient_(dict_.lookup("pressureGradient")),
-    UName_(dict_.lookupOrDefault<word>("U", "U"))
+    paramagnetic_(dict_.lookup("paramagnetic")),
+    chi_(0.0),
+    UName_(dict_.lookupOrDefault<word>("U", "U")),
+    HName_(dict_.lookupOrDefault<word>("H", "H"))
 {
     if (virtualMass_)
     {
         dict_.lookup("Cvm") >> Cvm_;
     }
+
+    if (paramagnetic_)
+    {
+        dict_.lookup("chi") >> chi_;
+    }
 }
 
 
@@ -60,11 +71,15 @@ Foam::particleForces::particleForces(const particleForces& f)
     dict_(f.dict_),
     g_(f.g_),
     gradUPtr_(f.gradUPtr_),
+    gradHPtr_(f.gradHPtr_),
     gravity_(f.gravity_),
     virtualMass_(f.virtualMass_),
     Cvm_(f.Cvm_),
     pressureGradient_(f.pressureGradient_),
-    UName_(f.UName_)
+    paramagnetic_(f.paramagnetic_),
+    chi_(f.chi_),
+    UName_(f.UName_),
+    HName_(f.HName_)
 {}
 
 
@@ -102,18 +117,42 @@ Foam::Switch Foam::particleForces::virtualMass() const
 }
 
 
+Foam::scalar Foam::particleForces::Cvm() const
+{
+    return Cvm_;
+}
+
+
 Foam::Switch Foam::particleForces::pressureGradient() const
 {
     return pressureGradient_;
 }
 
 
+Foam::Switch Foam::particleForces::paramagnetic() const
+{
+    return paramagnetic_;
+}
+
+
+Foam::scalar Foam::particleForces::chi() const
+{
+    return chi_;
+}
+
+
 const Foam::word& Foam::particleForces::UName() const
 {
     return UName_;
 }
 
 
+const Foam::word& Foam::particleForces::HName() const
+{
+    return HName_;
+}
+
+
 void Foam::particleForces::cacheFields(const bool store)
 {
     if (store && pressureGradient_)
@@ -129,6 +168,20 @@ void Foam::particleForces::cacheFields(const bool store)
             gradUPtr_ = NULL;
         }
     }
+
+    if (store && paramagnetic_)
+    {
+        const volVectorField& H = mesh_.lookupObject<volVectorField>(HName_);
+        gradHPtr_ = fvc::grad(H).ptr();
+    }
+    else
+    {
+        if (gradHPtr_)
+        {
+            delete gradHPtr_;
+            gradHPtr_ = NULL;
+        }
+    }
 }
 
 
@@ -139,7 +192,8 @@ Foam::vector Foam::particleForces::calcCoupled
     const scalar rhoc,
     const scalar rho,
     const vector& Uc,
-    const vector& U
+    const vector& U,
+    const scalar d
 ) const
 {
     vector Ftot = vector::zero;
@@ -172,7 +226,8 @@ Foam::vector Foam::particleForces::calcNonCoupled
     const scalar rhoc,
     const scalar rho,
     const vector& Uc,
-    const vector& U
+    const vector& U,
+    const scalar d
 ) const
 {
     vector Ftot = vector::zero;
@@ -183,9 +238,34 @@ Foam::vector Foam::particleForces::calcNonCoupled
         Ftot += g_*(1.0 - rhoc/rho);
     }
 
+    // Magnetic field force
+
+    if (paramagnetic_)
+    {
+        const volVectorField& H = mesh_.lookupObject<volVectorField>(HName_);
+
+        const volTensorField& gradH = *gradHPtr_;
+
+        Ftot +=
+            3.0*constant::electromagnetic::mu0.value()/rho
+           *chi_/(chi_ + 3)
+           *(H[cellI] & gradH[cellI]);
+
+        // force is:
+
+        // 4.0
+        // *constant::mathematical::pi
+        // *constant::electromagnetic::mu0.value()
+        // *pow3(d/2)
+        // *chi/(chi + 3)
+        // *(H[cellI] & gradH[cellI]);
+
+        // which is divided by mass ((4/3)*pi*r^3*rho) to produce
+        // acceleration
+    }
+
     return Ftot;
 }
 
 
 // ************************************************************************* //
-
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H
index 87abf7a7d0..19c7b1b784 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.H
+++ b/src/lagrangian/intermediate/particleForces/particleForces.H
@@ -69,6 +69,9 @@ class particleForces
         //- Velocity gradient field
         const volTensorField* gradUPtr_;
 
+        //- Magnetic field strength gradient field
+        const volTensorField* gradHPtr_;
+
 
         // Forces to include in particle motion evaluation
 
@@ -84,12 +87,21 @@ class particleForces
             //- Pressure gradient
             Switch pressureGradient_;
 
+            //- Paramagnetic force
+            Switch paramagnetic_;
+
+            //- Magnetic susceptibility of particle
+            scalar chi_;
+
 
         // Additional info
 
-            //- Name of velucity field - default = "U"
+            //- Name of velocity field - default = "U"
             const word UName_;
 
+            //- Name of magnetic field strength field - default = "H"
+            const word HName_;
+
 
 public:
 
@@ -128,7 +140,7 @@ public:
             Switch virtualMass() const;
 
             //- Return virtual mass force coefficient
-            Switch Cvm() const;
+            scalar Cvm() const;
 
             //- Return pressure gradient force activate switch
             Switch pressureGradient() const;
@@ -136,6 +148,15 @@ public:
             //- Return name of velocity field
             const word& UName() const;
 
+            //- Return paramagnetic force activate switch
+            Switch paramagnetic() const;
+
+            //- Return magnetic susceptibility
+            scalar chi() const;
+
+            //- Return name of velocity field
+            const word& HName() const;
+
 
        // Evaluation
 
@@ -150,7 +171,8 @@ public:
                 const scalar rhoc,
                 const scalar rho,
                 const vector& Uc,
-                const vector& U
+                const vector& U,
+                const scalar d
             ) const;
 
             //- Calculate external forces applied to the particles
@@ -161,7 +183,8 @@ public:
                 const scalar rhoc,
                 const scalar rho,
                 const vector& Uc,
-                const vector& U
+                const vector& U,
+                const scalar d
             ) const;
 };
 

From 55cd7409d8d3565df68651546c3bae86ae49b7b5 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 14 Apr 2010 09:51:46 +0100
Subject: [PATCH 02/12] STYLE:  particleForces.  Rename variable chi ->
 magneticSusceptibility.

---
 .../intermediate/particleForces/particleForces.C   | 14 +++++++-------
 .../intermediate/particleForces/particleForces.H   |  4 ++--
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/lagrangian/intermediate/particleForces/particleForces.C b/src/lagrangian/intermediate/particleForces/particleForces.C
index 938549a734..a9c11aa81f 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.C
+++ b/src/lagrangian/intermediate/particleForces/particleForces.C
@@ -49,7 +49,7 @@ Foam::particleForces::particleForces
     Cvm_(0.0),
     pressureGradient_(dict_.lookup("pressureGradient")),
     paramagnetic_(dict_.lookup("paramagnetic")),
-    chi_(0.0),
+    magneticSusceptibility_(0.0),
     UName_(dict_.lookupOrDefault<word>("U", "U")),
     HName_(dict_.lookupOrDefault<word>("H", "H"))
 {
@@ -60,7 +60,7 @@ Foam::particleForces::particleForces
 
     if (paramagnetic_)
     {
-        dict_.lookup("chi") >> chi_;
+        dict_.lookup("magneticSusceptibility") >> magneticSusceptibility_;
     }
 }
 
@@ -77,7 +77,7 @@ Foam::particleForces::particleForces(const particleForces& f)
     Cvm_(f.Cvm_),
     pressureGradient_(f.pressureGradient_),
     paramagnetic_(f.paramagnetic_),
-    chi_(f.chi_),
+    magneticSusceptibility_(f.magneticSusceptibility_),
     UName_(f.UName_),
     HName_(f.HName_)
 {}
@@ -135,9 +135,9 @@ Foam::Switch Foam::particleForces::paramagnetic() const
 }
 
 
-Foam::scalar Foam::particleForces::chi() const
+Foam::scalar Foam::particleForces::magneticSusceptibility() const
 {
-    return chi_;
+    return magneticSusceptibility_;
 }
 
 
@@ -248,7 +248,7 @@ Foam::vector Foam::particleForces::calcNonCoupled
 
         Ftot +=
             3.0*constant::electromagnetic::mu0.value()/rho
-           *chi_/(chi_ + 3)
+           *magneticSusceptibility_/(magneticSusceptibility_ + 3)
            *(H[cellI] & gradH[cellI]);
 
         // force is:
@@ -257,7 +257,7 @@ Foam::vector Foam::particleForces::calcNonCoupled
         // *constant::mathematical::pi
         // *constant::electromagnetic::mu0.value()
         // *pow3(d/2)
-        // *chi/(chi + 3)
+        // *magneticSusceptibility/(magneticSusceptibility + 3)
         // *(H[cellI] & gradH[cellI]);
 
         // which is divided by mass ((4/3)*pi*r^3*rho) to produce
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H
index 19c7b1b784..d7e891ab65 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.H
+++ b/src/lagrangian/intermediate/particleForces/particleForces.H
@@ -91,7 +91,7 @@ class particleForces
             Switch paramagnetic_;
 
             //- Magnetic susceptibility of particle
-            scalar chi_;
+            scalar magneticSusceptibility_;
 
 
         // Additional info
@@ -152,7 +152,7 @@ public:
             Switch paramagnetic() const;
 
             //- Return magnetic susceptibility
-            scalar chi() const;
+            scalar magneticSusceptibility() const;
 
             //- Return name of velocity field
             const word& HName() const;

From c948fadb10bb06282d786f0528515056660d84cc Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 14 Apr 2010 13:59:58 +0100
Subject: [PATCH 03/12] ENH: particleForces.  Not caching gradH field, looking
 up HdotGradH field directly.

Renaming Ftot to accelTot to remind that the functions return
accelerations.
---
 .../particleForces/particleForces.C           | 53 +++++++------------
 .../particleForces/particleForces.H           | 12 ++---
 2 files changed, 24 insertions(+), 41 deletions(-)

diff --git a/src/lagrangian/intermediate/particleForces/particleForces.C b/src/lagrangian/intermediate/particleForces/particleForces.C
index a9c11aa81f..7778736683 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.C
+++ b/src/lagrangian/intermediate/particleForces/particleForces.C
@@ -43,7 +43,6 @@ Foam::particleForces::particleForces
     dict_(dict.subDict("particleForces")),
     g_(g),
     gradUPtr_(NULL),
-    gradHPtr_(NULL),
     gravity_(dict_.lookup("gravity")),
     virtualMass_(dict_.lookup("virtualMass")),
     Cvm_(0.0),
@@ -51,7 +50,7 @@ Foam::particleForces::particleForces
     paramagnetic_(dict_.lookup("paramagnetic")),
     magneticSusceptibility_(0.0),
     UName_(dict_.lookupOrDefault<word>("U", "U")),
-    HName_(dict_.lookupOrDefault<word>("H", "H"))
+    HdotGradHName_(dict_.lookupOrDefault<word>("HdotGradH", "HdotGradH"))
 {
     if (virtualMass_)
     {
@@ -71,7 +70,6 @@ Foam::particleForces::particleForces(const particleForces& f)
     dict_(f.dict_),
     g_(f.g_),
     gradUPtr_(f.gradUPtr_),
-    gradHPtr_(f.gradHPtr_),
     gravity_(f.gravity_),
     virtualMass_(f.virtualMass_),
     Cvm_(f.Cvm_),
@@ -79,7 +77,7 @@ Foam::particleForces::particleForces(const particleForces& f)
     paramagnetic_(f.paramagnetic_),
     magneticSusceptibility_(f.magneticSusceptibility_),
     UName_(f.UName_),
-    HName_(f.HName_)
+    HdotGradHName_(f.HdotGradHName_)
 {}
 
 
@@ -147,9 +145,9 @@ const Foam::word& Foam::particleForces::UName() const
 }
 
 
-const Foam::word& Foam::particleForces::HName() const
+const Foam::word& Foam::particleForces::HdotGradHName() const
 {
-    return HName_;
+    return HdotGradHName_;
 }
 
 
@@ -168,20 +166,6 @@ void Foam::particleForces::cacheFields(const bool store)
             gradUPtr_ = NULL;
         }
     }
-
-    if (store && paramagnetic_)
-    {
-        const volVectorField& H = mesh_.lookupObject<volVectorField>(HName_);
-        gradHPtr_ = fvc::grad(H).ptr();
-    }
-    else
-    {
-        if (gradHPtr_)
-        {
-            delete gradHPtr_;
-            gradHPtr_ = NULL;
-        }
-    }
 }
 
 
@@ -196,7 +180,7 @@ Foam::vector Foam::particleForces::calcCoupled
     const scalar d
 ) const
 {
-    vector Ftot = vector::zero;
+    vector accelTot = vector::zero;
 
     // Virtual mass force
     if (virtualMass_)
@@ -205,17 +189,17 @@ Foam::vector Foam::particleForces::calcCoupled
         (
             "Foam::particleForces::calcCoupled(...) - virtual mass force"
         );
-//        Ftot += Cvm_*rhoc/rho*d(Uc - U)/dt;
+//        accelTot += Cvm_*rhoc/rho*d(Uc - U)/dt;
     }
 
     // Pressure gradient force
     if (pressureGradient_)
     {
         const volTensorField& gradU = *gradUPtr_;
-        Ftot += rhoc/rho*(U & gradU[cellI]);
+        accelTot += rhoc/rho*(U & gradU[cellI]);
     }
 
-    return Ftot;
+    return accelTot;
 }
 
 
@@ -230,26 +214,27 @@ Foam::vector Foam::particleForces::calcNonCoupled
     const scalar d
 ) const
 {
-    vector Ftot = vector::zero;
+    vector accelTot = vector::zero;
 
     // Gravity force
     if (gravity_)
     {
-        Ftot += g_*(1.0 - rhoc/rho);
+        accelTot += g_*(1.0 - rhoc/rho);
     }
 
     // Magnetic field force
 
     if (paramagnetic_)
     {
-        const volVectorField& H = mesh_.lookupObject<volVectorField>(HName_);
+        const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
+        (
+            HdotGradHName_
+        );
 
-        const volTensorField& gradH = *gradHPtr_;
-
-        Ftot +=
+        accelTot +=
             3.0*constant::electromagnetic::mu0.value()/rho
            *magneticSusceptibility_/(magneticSusceptibility_ + 3)
-           *(H[cellI] & gradH[cellI]);
+           *HdotGradH[cellI];
 
         // force is:
 
@@ -257,14 +242,14 @@ Foam::vector Foam::particleForces::calcNonCoupled
         // *constant::mathematical::pi
         // *constant::electromagnetic::mu0.value()
         // *pow3(d/2)
-        // *magneticSusceptibility/(magneticSusceptibility + 3)
-        // *(H[cellI] & gradH[cellI]);
+        // *magneticSusceptibility_/(magneticSusceptibility_ + 3)
+        // *HdotGradH[cellI];
 
         // which is divided by mass ((4/3)*pi*r^3*rho) to produce
         // acceleration
     }
 
-    return Ftot;
+    return accelTot;
 }
 
 
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H
index d7e891ab65..cc47366533 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.H
+++ b/src/lagrangian/intermediate/particleForces/particleForces.H
@@ -69,9 +69,6 @@ class particleForces
         //- Velocity gradient field
         const volTensorField* gradUPtr_;
 
-        //- Magnetic field strength gradient field
-        const volTensorField* gradHPtr_;
-
 
         // Forces to include in particle motion evaluation
 
@@ -99,8 +96,9 @@ class particleForces
             //- Name of velocity field - default = "U"
             const word UName_;
 
-            //- Name of magnetic field strength field - default = "H"
-            const word HName_;
+            //- Name of paramagnetic field strength field - default =
+            //  "HdotGradH"
+            const word HdotGradHName_;
 
 
 public:
@@ -154,8 +152,8 @@ public:
             //- Return magnetic susceptibility
             scalar magneticSusceptibility() const;
 
-            //- Return name of velocity field
-            const word& HName() const;
+            //- Return name of paramagnetic field strength field
+            const word& HdotGradHName() const;
 
 
        // Evaluation

From 93a0172f038407709263c07be637f79ca2a374a5 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 14 Apr 2010 15:47:39 +0100
Subject: [PATCH 04/12] BUG: incorrect statistics in parallel. Also removed
 sign of magSf.

---
 .../field/fieldValues/cellSource/cellSource.C | 10 +++---
 .../field/fieldValues/cellSource/cellSource.H |  3 ++
 .../field/fieldValues/faceSource/faceSource.C | 31 +++++++++--------
 .../field/fieldValues/faceSource/faceSource.H | 25 ++++++++++----
 .../fieldValues/faceSource/faceSourceI.H      |  4 +--
 .../faceSource/faceSourceTemplates.C          | 33 ++++++++++++++-----
 6 files changed, 73 insertions(+), 33 deletions(-)

diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
index 9d65cefc57..97cd9393bd 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
@@ -80,10 +80,11 @@ void Foam::fieldValues::cellSource::setCellZoneCells()
     }
 
     cellId_.setSize(count);
+    nCells_ = returnReduce(cellId_.size(), sumOp<label>());
 
     if (debug)
     {
-        Info<< "Original cell zone size = " << cZone.size()
+        Pout<< "Original cell zone size = " << cZone.size()
             << ", new size = " << count << endl;
     }
 }
@@ -109,8 +110,8 @@ void Foam::fieldValues::cellSource::initialise(const dictionary& dict)
     }
 
     Info<< type() << " " << name_ << ":" << nl
-        << "    total cells  = " << cellId_.size() << nl
-        << "    total volume = " << sum(filterField(mesh().V()))
+        << "    total cells  = " << nCells_ << nl
+        << "    total volume = " << gSum(filterField(mesh().V()))
         << nl << endl;
 
     if (operation_ == opWeightedAverage)
@@ -144,7 +145,7 @@ void Foam::fieldValues::cellSource::writeFileHeader()
     {
         outputFilePtr_()
             << "# Source : " << sourceTypeNames_[source_] << " "
-            << sourceName_ <<  nl << "# Cells  : " << cellId_.size() << nl
+            << sourceName_ <<  nl << "# Cells  : " << nCells_ << nl
             << "# Time" << tab << "sum(V)";
 
         forAll(fields_, i)
@@ -172,6 +173,7 @@ Foam::fieldValues::cellSource::cellSource
     fieldValue(name, obr, dict, loadFromFiles),
     source_(sourceTypeNames_.read(dict.lookup("source"))),
     operation_(operationTypeNames_.read(dict.lookup("operation"))),
+    nCells_(0),
     cellId_()
 {
     read(dict);
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H
index 8108f7891c..c271a1c997 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H
@@ -133,6 +133,9 @@ protected:
         //- Operation to apply to values
         operationType operation_;
 
+        //- Global number of cells
+        label nCells_;
+
         //- Local list of cell IDs
         labelList cellId_;
 
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
index b6788e1af2..e19a80dc2a 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
@@ -78,7 +78,7 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
 
     faceId_.setSize(fZone.size());
     facePatchId_.setSize(fZone.size());
-    flipMap_.setSize(fZone.size());
+    faceSign_.setSize(fZone.size());
 
     label count = 0;
     forAll(fZone, i)
@@ -134,11 +134,11 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
         {
             if (fZone.flipMap()[i])
             {
-                flipMap_[count] = -1;
+                faceSign_[count] = -1;
             }
             else
             {
-                flipMap_[count] = 1;
+                faceSign_[count] = 1;
             }
             faceId_[count] = faceId;
             facePatchId_[count] = facePatchId;
@@ -148,11 +148,12 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
 
     faceId_.setSize(count);
     facePatchId_.setSize(count);
-    flipMap_.setSize(count);
+    faceSign_.setSize(count);
+    nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
 
     if (debug)
     {
-        Info<< "Original face zone size = " << fZone.size()
+        Pout<< "Original face zone size = " << fZone.size()
             << ", new size = " << count << endl;
     }
 }
@@ -187,13 +188,14 @@ void Foam::fieldValues::faceSource::setPatchFaces()
 
     faceId_.setSize(nFaces);
     facePatchId_.setSize(nFaces);
-    flipMap_.setSize(nFaces);
+    faceSign_.setSize(nFaces);
+    nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
 
     forAll(faceId_, faceI)
     {
         faceId_[faceI] = faceI;
         facePatchId_[faceI] = patchId;
-        flipMap_[faceI] = 1;
+        faceSign_[faceI] = 1;
     }
 }
 
@@ -216,7 +218,7 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
         }
         default:
         {
-            FatalErrorIn("faceSource::initiliase()")
+            FatalErrorIn("faceSource::initialise()")
                 << type() << " " << name_ << ": "
                 << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
                 << nl << "    Unknown source type. Valid source types are:"
@@ -225,8 +227,10 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
     }
 
     Info<< type() << " " << name_ << ":" << nl
-        << "    total faces  = " << faceId_.size() << nl
-        << "    total area   = " << sum(filterField(mesh().magSf())) << nl;
+        << "    total faces  = " << nFaces_
+        << nl
+        << "    total area   = " << gSum(filterField(mesh().magSf(), false))
+        << nl;
 
     if (operation_ == opWeightedAverage)
     {
@@ -260,7 +264,7 @@ void Foam::fieldValues::faceSource::writeFileHeader()
     {
         outputFilePtr_()
             << "# Source : " << sourceTypeNames_[source_] << " "
-            << sourceName_ <<  nl << "# Faces  : " << faceId_.size() << nl
+            << sourceName_ <<  nl << "# Faces  : " << nFaces_ << nl
             << "# Time" << tab << "sum(magSf)";
 
         forAll(fields_, i)
@@ -288,9 +292,10 @@ Foam::fieldValues::faceSource::faceSource
     fieldValue(name, obr, dict, loadFromFiles),
     source_(sourceTypeNames_.read(dict.lookup("source"))),
     operation_(operationTypeNames_.read(dict.lookup("operation"))),
+    nFaces_(0),
     faceId_(),
     facePatchId_(),
-    flipMap_(),
+    faceSign_(),
     weightFieldName_("undefinedWeightedFieldName")
 {
     read(dict);
@@ -326,7 +331,7 @@ void Foam::fieldValues::faceSource::write()
         {
             outputFilePtr_()
                 << obr_.time().value() << tab
-                << sum(filterField(mesh().magSf()));
+                << sum(filterField(mesh().magSf(), false));
         }
 
         forAll(fields_, i)
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
index e80e1b31ea..dc3d092334 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
@@ -53,6 +53,15 @@ Description
       - areaAverage
       - areaIntegrate
       - weightedAverage
+      - min
+      - max
+
+    Notes:
+    - faces on empty patches get ignored
+    - if the field is a volField the faceZone can only consist of boundary
+    faces.
+    - all fields get oriented according to the faceZone (so you might e.g. see
+    negative pressure)
 
 SourceFiles
     faceSource.C
@@ -64,7 +73,6 @@ SourceFiles
 
 #include "NamedEnum.H"
 #include "fieldValue.H"
-#include "labelList.H"
 #include "surfaceFieldsFwd.H"
 #include "volFieldsFwd.H"
 
@@ -136,14 +144,17 @@ protected:
         //- Operation to apply to values
         operationType operation_;
 
+        //- Global number of faces
+        label nFaces_;
+
         //- Local list of face IDs
         labelList faceId_;
 
         //- Local list of patch ID per face
         labelList facePatchId_;
 
-        //- List of +1/-1 representing face flip map
-        labelList flipMap_;
+        //- List of +1/-1 representing face flip map (1 use as is, -1 negate)
+        labelList faceSign_;
 
         //- Weight field name - only used for opWeightedAverage mode
         word weightFieldName_;
@@ -209,7 +220,7 @@ public:
             inline const labelList& facePatch() const;
 
             //- Return the list of +1/-1 representing face flip map
-            inline const labelList& flipMap() const;
+            inline const labelList& faceSign() const;
 
 
         // Function object functions
@@ -228,14 +239,16 @@ public:
             template<class Type>
             tmp<Field<Type> > filterField
             (
-                const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+                const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
+                const bool applyOrientation
             ) const;
 
             //- Filter a volume field according to faceIds
             template<class Type>
             tmp<Field<Type> > filterField
             (
-                const GeometricField<Type, fvPatchField, volMesh>& field
+                const GeometricField<Type, fvPatchField, volMesh>& field,
+                const bool applyOrientation
             ) const;
 };
 
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceI.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceI.H
index 18ab0766eb..2ef320febb 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceI.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceI.H
@@ -49,9 +49,9 @@ Foam::fieldValues::faceSource::facePatch() const
 
 
 inline const Foam::labelList&
-Foam::fieldValues::faceSource::flipMap() const
+Foam::fieldValues::faceSource::faceSign() const
 {
-    return flipMap_;
+    return faceSign_;
 }
 
 
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
index 7c9ab9dbec..dfa293a082 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
@@ -60,14 +60,14 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::setFieldValues
 
     if (obr_.foundObject<sf>(fieldName))
     {
-        return filterField(obr_.lookupObject<sf>(fieldName));
+        return filterField(obr_.lookupObject<sf>(fieldName), true);
     }
     else if (obr_.foundObject<vf>(fieldName))
     {
-        return filterField(obr_.lookupObject<vf>(fieldName));
+        return filterField(obr_.lookupObject<vf>(fieldName), true);
     }
 
-    return tmp<Field<Type> >(new Field<Type>(0.0));
+    return tmp<Field<Type> >(new Field<Type>(0));
 }
 
 
@@ -131,10 +131,13 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
 
     if (ok)
     {
+        // Get (correctly oriented) field
         Field<Type> values = combineFields(setFieldValues<Type>(fieldName));
 
-        scalarField magSf = combineFields(filterField(mesh().magSf()));
+        // Get unoriented magSf
+        scalarField magSf = combineFields(filterField(mesh().magSf(), false));
 
+        // Get (correctly oriented) weighting field
         scalarField weightField =
             combineFields(setFieldValues<scalar>(weightFieldName_));
 
@@ -177,7 +180,8 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
 (
-    const GeometricField<Type, fvPatchField, volMesh>& field
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const bool applyOrientation
 ) const
 {
     tmp<Field<Type> > tvalues(new Field<Type>(faceId_.size()));
@@ -205,8 +209,14 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
                 << "    Unable to process internal faces for volume field "
                 << field.name() << nl << abort(FatalError);
         }
+    }
 
-        values[i] *= flipMap_[i];
+    if (applyOrientation)
+    {
+        forAll(values, i)
+        {
+            values[i] *= faceSign_[i];
+        }
     }
 
     return tvalues;
@@ -216,7 +226,8 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
 (
-    const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
+    const bool applyOrientation
 ) const
 {
     tmp<Field<Type> > tvalues(new Field<Type>(faceId_.size()));
@@ -234,8 +245,14 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
         {
             values[i] = field[faceI];
         }
+    }
 
-        values[i] *= flipMap_[i];
+    if (applyOrientation)
+    {
+        forAll(values, i)
+        {
+            values[i] *= faceSign_[i];
+        }
     }
 
     return tvalues;

From 36e845a7a3a24ae5b9b8319f2cd6cc6f15191476 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 14 Apr 2010 16:10:31 +0100
Subject: [PATCH 05/12] ENH: Allow derivation of constraint types (for
 pointFields only). Allow usage of non constraint patchFields on constraint
 patch types (all fields)

---
 .../GeometricField/GeometricBoundaryField.C   | 53 +++++++++++----
 .../GeometricField/GeometricField.C           | 22 +++++--
 .../GeometricField/GeometricField.H           | 15 +++--
 .../pointPatchField/newPointPatchField.C      | 65 ++++++++++++-------
 .../pointPatchField/pointPatchField.H         | 14 +++-
 .../fvPatchFields/fvPatchField/fvPatchField.H | 12 ++++
 .../fvPatchField/newFvPatchField.C            | 41 +++++++++---
 .../fvsPatchField/fvsPatchField.H             | 12 ++++
 .../fvsPatchField/newFvsPatchField.C          | 42 +++++++++---
 9 files changed, 210 insertions(+), 66 deletions(-)

diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
index 705505a716..7fa8d4b8bc 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
@@ -72,7 +72,8 @@ GeometricBoundaryField
 (
     const BoundaryMesh& bmesh,
     const DimensionedField<Type, GeoMesh>& field,
-    const wordList& patchFieldTypes
+    const wordList& patchFieldTypes,
+    const wordList& constraintTypes
 )
 :
     FieldField<PatchField, Type>(bmesh.size()),
@@ -83,18 +84,22 @@ GeometricBoundaryField
         Info<< "GeometricField<Type, PatchField, GeoMesh>::"
                "GeometricBoundaryField::"
                "GeometricBoundaryField(const BoundaryMesh&, "
-               "const Field<Type>&, const wordList&)"
+               "const Field<Type>&, const wordList&, const wordList&)"
             << endl;
     }
 
-    if (patchFieldTypes.size() != this->size())
+    if
+    (
+        patchFieldTypes.size() != this->size()
+     || (constraintTypes.size() && (constraintTypes.size() != this->size()))
+    )
     {
         FatalErrorIn
         (
             "GeometricField<Type, PatchField, GeoMesh>::"
             "GeometricBoundaryField::"
             "GeometricBoundaryField(const BoundaryMesh&, "
-            "const Field<Type>&, const wordList&)"
+            "const Field<Type>&, const wordList&, const wordList&)"
         )   << "Incorrect number of patch type specifications given" << nl
             << "    Number of patches in mesh = " << bmesh.size()
             << " number of patch type specifications = "
@@ -102,18 +107,38 @@ GeometricBoundaryField
             << abort(FatalError);
     }
 
-    forAll(bmesh_, patchi)
+    if (constraintTypes.size())
     {
-        set
-        (
-            patchi,
-            PatchField<Type>::New
+        forAll(bmesh_, patchi)
+        {
+            set
             (
-                patchFieldTypes[patchi],
-                bmesh_[patchi],
-                field
-            )
-        );
+                patchi,
+                PatchField<Type>::New
+                (
+                    patchFieldTypes[patchi],
+                    constraintTypes[patchi],
+                    bmesh_[patchi],
+                    field
+                )
+            );
+        }
+    }
+    else
+    {
+        forAll(bmesh_, patchi)
+        {
+            set
+            (
+                patchi,
+                PatchField<Type>::New
+                (
+                    patchFieldTypes[patchi],
+                    bmesh_[patchi],
+                    field
+                )
+            );
+        }
     }
 }
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 8c17c3886b..b4960c341d 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -231,14 +231,15 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const IOobject& io,
     const Mesh& mesh,
     const dimensionSet& ds,
-    const wordList& patchFieldTypes
+    const wordList& patchFieldTypes,
+    const wordList& actualPatchTypes
 )
 :
     DimensionedField<Type, GeoMesh>(io, mesh, ds),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
-    boundaryField_(mesh.boundary(), *this, patchFieldTypes)
+    boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 {
     if (debug)
     {
@@ -287,14 +288,15 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const IOobject& io,
     const Mesh& mesh,
     const dimensioned<Type>& dt,
-    const wordList& patchFieldTypes
+    const wordList& patchFieldTypes,
+    const wordList& actualPatchTypes
 )
 :
     DimensionedField<Type, GeoMesh>(io, mesh, dt),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
-    boundaryField_(mesh.boundary(), *this, patchFieldTypes)
+    boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 {
     if (debug)
     {
@@ -653,14 +655,22 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 (
     const IOobject& io,
     const GeometricField<Type, PatchField, GeoMesh>& gf,
-    const wordList& patchFieldTypes
+    const wordList& patchFieldTypes,
+    const wordList& actualPatchTypes
+
 )
 :
     DimensionedField<Type, GeoMesh>(io, gf),
     timeIndex_(gf.timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
-    boundaryField_(this->mesh().boundary(), *this, patchFieldTypes)
+    boundaryField_
+    (
+        this->mesh().boundary(),
+        *this,
+        patchFieldTypes,
+        actualPatchTypes
+    )
 {
     if (debug)
     {
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 7d529ab40f..fa8ec179e4 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -128,12 +128,14 @@ public:
 
             //- Construct from a BoundaryMesh,
             //  reference to the internal field
-            //  and a wordList of patch types
+            //  and a wordList of patch types and optional the actual patch
+            //  types (to override constraint patches)
             GeometricBoundaryField
             (
                 const BoundaryMesh&,
                 const DimensionedInternalField&,
-                const wordList&
+                const wordList& wantedPatchTypes,
+                const wordList& actualPatchTypes = wordList()
             );
 
             //- Construct from a BoundaryMesh,
@@ -283,7 +285,8 @@ public:
             const IOobject&,
             const Mesh&,
             const dimensionSet&,
-            const wordList& patchFieldTypes
+            const wordList& wantedPatchTypes,
+            const wordList& actualPatchTypes = wordList()
         );
 
         //- Constructor given IOobject, mesh, dimensioned<Type> and patch type.
@@ -301,7 +304,8 @@ public:
             const IOobject&,
             const Mesh&,
             const dimensioned<Type>&,
-            const wordList& patchFieldTypes
+            const wordList& wantedPatchTypes,
+            const wordList& actualPatchTypes = wordList()
         );
 
         //- Constructor from components
@@ -387,7 +391,8 @@ public:
         (
             const IOobject&,
             const GeometricField<Type, PatchField, GeoMesh>&,
-            const wordList& patchFieldTypes
+            const wordList& patchFieldTypes,
+            const wordList& actualPatchTypes = wordList()
         );
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/newPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/newPointPatchField.C
index b14fbeac8f..749511b2cf 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/newPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/newPointPatchField.C
@@ -29,6 +29,7 @@ template<class Type>
 Foam::autoPtr<Foam::pointPatchField<Type> > Foam::pointPatchField<Type>::New
 (
     const word& patchFieldType,
+    const word& actualPatchType,
     const pointPatch& p,
     const DimensionedField<Type, pointMesh>& iF
 )
@@ -36,7 +37,8 @@ Foam::autoPtr<Foam::pointPatchField<Type> > Foam::pointPatchField<Type>::New
     if (debug)
     {
         Info<< "PointPatchField<Type>::"
-               "New(const word&, const pointPatch&, const Field<Type>&) : "
+               "New(const word&, const word&"
+                ", const pointPatch&, const Field<Type>&) : "
                "constructing pointPatchField<Type>"
             << endl;
     }
@@ -49,7 +51,7 @@ Foam::autoPtr<Foam::pointPatchField<Type> > Foam::pointPatchField<Type>::New
         FatalErrorIn
         (
             "PointPatchField<Type>::New"
-            "(const word&, const pointPatch&, const Field<Type>&)"
+            "(const word&, const word&, const pointPatch&, const Field<Type>&)"
         )   << "Unknown patchFieldType type "
             << patchFieldType
             << endl << endl
@@ -60,31 +62,48 @@ Foam::autoPtr<Foam::pointPatchField<Type> > Foam::pointPatchField<Type>::New
 
     autoPtr<pointPatchField<Type> > pfPtr(cstrIter()(p, iF));
 
-    if (pfPtr().constraintType() == p.constraintType())
+    if
+    (
+        actualPatchType == word::null
+     || actualPatchType != p.type()
+    )
     {
-        // Compatible (constraint-wise) with the patch type
-        return pfPtr;
-    }
-    else
-    {
-        // Use default constraint type
-        typename pointPatchConstructorTable::iterator patchTypeCstrIter =
-            pointPatchConstructorTablePtr_->find(p.type());
-
-        if (patchTypeCstrIter == pointPatchConstructorTablePtr_->end())
+        if (pfPtr().constraintType() != p.constraintType())
         {
-            FatalErrorIn
-            (
-                "PointPatchField<Type>::New"
-                "(const word&, const pointPatch&, const Field<Type>&)"
-            )   << "inconsistent patch and patchField types for \n"
-                << "    patch type " << p.type()
-                << " and patchField type " << patchFieldType
-                << exit(FatalError);
-        }
+            // Use default constraint type
+            typename pointPatchConstructorTable::iterator patchTypeCstrIter =
+                pointPatchConstructorTablePtr_->find(p.type());
 
-        return patchTypeCstrIter()(p, iF);
+            if (patchTypeCstrIter == pointPatchConstructorTablePtr_->end())
+            {
+                FatalErrorIn
+                (
+                    "PointPatchField<Type>::New"
+                    "(const word&, const word&"
+                    ", const pointPatch&, const Field<Type>&)"
+                )   << "inconsistent patch and patchField types for \n"
+                    << "    patch type " << p.type()
+                    << " and patchField type " << patchFieldType
+                    << exit(FatalError);
+            }
+
+            return patchTypeCstrIter()(p, iF);
+        }
     }
+
+    return pfPtr;
+}
+
+
+template<class Type>
+Foam::autoPtr<Foam::pointPatchField<Type> > Foam::pointPatchField<Type>::New
+(
+    const word& patchFieldType,
+    const pointPatch& p,
+    const DimensionedField<Type, pointMesh>& iF
+)
+{
+    return New(patchFieldType, word::null, p, iF);
 }
 
 
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
index 98604bb4bf..869c536365 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
@@ -35,7 +35,7 @@ Description
 
 SourceFiles
     pointPatchField.C
-    newpointPatchField.C
+    newPointPatchField.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -199,6 +199,18 @@ public:
             const DimensionedField<Type, pointMesh>&
         );
 
+        //- Return a pointer to a new patchField created on freestore given
+        //  patch and internal field
+        //  (does not set the patch field values).
+        //  Allows override of constraint type
+        static autoPtr<pointPatchField<Type> > New
+        (
+            const word&,
+            const word& actualPatchType,
+            const pointPatch&,
+            const DimensionedField<Type, pointMesh>&
+        );
+
         //- Return a pointer to a new patchField created on freestore from
         //  a given pointPatchField mapped onto a new patch
         static autoPtr<pointPatchField<Type> > New
diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
index 1a8720af06..b1df3530ee 100644
--- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
@@ -226,6 +226,18 @@ public:
             const DimensionedField<Type, volMesh>&
         );
 
+        //- Return a pointer to a new patchField created on freestore given
+        //  patch and internal field
+        //  (does not set the patch field values).
+        //  Allows override of constraint type
+        static tmp<fvPatchField<Type> > New
+        (
+            const word&,
+            const word& actualPatchType,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
         //- Return a pointer to a new patchField created on freestore from
         //  a given fvPatchField mapped onto a new patch
         static tmp<fvPatchField<Type> > New
diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/newFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/fvPatchField/newFvPatchField.C
index 8f89f202d9..86884666a9 100644
--- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/newFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/newFvPatchField.C
@@ -29,14 +29,16 @@ template<class Type>
 Foam::tmp<Foam::fvPatchField<Type> > Foam::fvPatchField<Type>::New
 (
     const word& patchFieldType,
+    const word& actualPatchType,
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
 )
 {
     if (debug)
     {
-        Info<< "fvPatchField<Type>::New(const word&, const fvPatch&, "
-               "const DimensionedField<Type, volMesh>&) : patchFieldType="
+        Info<< "fvPatchField<Type>::New(const word&, const word&, "
+               "const fvPatch&, const DimensionedField<Type, volMesh>&) :"
+               " patchFieldType="
             << patchFieldType
             << endl;
     }
@@ -48,7 +50,7 @@ Foam::tmp<Foam::fvPatchField<Type> > Foam::fvPatchField<Type>::New
     {
         FatalErrorIn
         (
-            "fvPatchField<Type>::New(const word&, const fvPatch&, "
+            "fvPatchField<Type>::New(const word&, const word&, const fvPatch&, "
             "const DimensionedField<Type, volMesh>&)"
         )   << "Unknown patchTypefield type " << patchFieldType
             << endl << endl
@@ -57,12 +59,23 @@ Foam::tmp<Foam::fvPatchField<Type> > Foam::fvPatchField<Type>::New
             << exit(FatalError);
     }
 
-    typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(p.type());
-
-    if (patchTypeCstrIter != patchConstructorTablePtr_->end())
+    if
+    (
+        actualPatchType == word::null
+     || actualPatchType != p.type()
+    )
     {
-        return patchTypeCstrIter()(p, iF);
+        typename patchConstructorTable::iterator patchTypeCstrIter =
+            patchConstructorTablePtr_->find(p.type());
+
+        if (patchTypeCstrIter != patchConstructorTablePtr_->end())
+        {
+            return patchTypeCstrIter()(p, iF);
+        }
+        else
+        {
+            return cstrIter()(p, iF);
+        }
     }
     else
     {
@@ -71,6 +84,18 @@ Foam::tmp<Foam::fvPatchField<Type> > Foam::fvPatchField<Type>::New
 }
 
 
+template<class Type>
+Foam::tmp<Foam::fvPatchField<Type> > Foam::fvPatchField<Type>::New
+(
+    const word& patchFieldType,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+{
+    return New(patchFieldType, word::null, p, iF);
+}
+
+
 template<class Type>
 Foam::tmp<Foam::fvPatchField<Type> > Foam::fvPatchField<Type>::New
 (
diff --git a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H
index b98fdcd453..c400d78d20 100644
--- a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H
+++ b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H
@@ -216,6 +216,18 @@ public:
             const DimensionedField<Type, surfaceMesh>&
         );
 
+        //- Return a pointer to a new patchField created on freestore given
+        //  patch and internal field
+        //  (does not set the patch field values)
+        //  Allows override of constraint type
+        static tmp<fvsPatchField<Type> > New
+        (
+            const word&,
+            const word& actualPatchType, 
+            const fvPatch&,
+            const DimensionedField<Type, surfaceMesh>&
+        );
+
         //- Return a pointer to a new patchField created on freestore from
         //  a given fvsPatchField mapped onto a new patch
         static tmp<fvsPatchField<Type> > New
diff --git a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/newFvsPatchField.C b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/newFvsPatchField.C
index b33ee35043..b8ad302fdd 100644
--- a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/newFvsPatchField.C
+++ b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/newFvsPatchField.C
@@ -34,14 +34,15 @@ template<class Type>
 tmp<fvsPatchField<Type> > fvsPatchField<Type>::New
 (
     const word& patchFieldType,
+    const word& actualPatchType,
     const fvPatch& p,
     const DimensionedField<Type, surfaceMesh>& iF
 )
 {
     if (debug)
     {
-        Info<< "fvsPatchField<Type>::New(const word&, const fvPatch&, "
-               "const Field<Type>&) : "
+        Info<< "fvsPatchField<Type>::New(const word&, const word&"
+               ", const fvPatch&, const Field<Type>&) : "
                "constructing fvsPatchField<Type>"
             << endl;
     }
@@ -53,8 +54,8 @@ tmp<fvsPatchField<Type> > fvsPatchField<Type>::New
     {
         FatalErrorIn
         (
-            "fvsPatchField<Type>::New(const word&, const fvPatch&, "
-            "const Field<Type>&)"
+            "fvsPatchField<Type>::New(const word&, const word&, const fvPatch&"
+            ", const Field<Type>&)"
         )   << "Unknown patchTypefield type " << patchFieldType
             << endl << endl
             << "Valid patchField types are :" << endl
@@ -62,12 +63,23 @@ tmp<fvsPatchField<Type> > fvsPatchField<Type>::New
             << exit(FatalError);
     }
 
-    typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(p.type());
-
-    if (patchTypeCstrIter != patchConstructorTablePtr_->end())
+    if
+    (
+        actualPatchType == word::null
+     || actualPatchType != p.type()
+    )
     {
-        return patchTypeCstrIter()(p, iF);
+        typename patchConstructorTable::iterator patchTypeCstrIter =
+            patchConstructorTablePtr_->find(p.type());
+
+        if (patchTypeCstrIter != patchConstructorTablePtr_->end())
+        {
+            return patchTypeCstrIter()(p, iF);
+        }
+        else
+        {
+            return cstrIter()(p, iF);
+        }
     }
     else
     {
@@ -76,6 +88,18 @@ tmp<fvsPatchField<Type> > fvsPatchField<Type>::New
 }
 
 
+template<class Type>
+tmp<fvsPatchField<Type> > fvsPatchField<Type>::New
+(
+    const word& patchFieldType,
+    const fvPatch& p,
+    const DimensionedField<Type, surfaceMesh>& iF
+)
+{
+    return New(patchFieldType, word::null, p, iF);
+}
+
+
 template<class Type>
 tmp<fvsPatchField<Type> > fvsPatchField<Type>::New
 (

From 6d71a38e6673ff1b234b71fa46dad900e2ea7902 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 11:34:47 +0100
Subject: [PATCH 06/12] ENH: Allow tensorial K

---
 ...allHeatFluxTemperatureFvPatchScalarField.C | 42 ++++++++++++++++---
 ...allHeatFluxTemperatureFvPatchScalarField.H | 12 ++++--
 2 files changed, 45 insertions(+), 9 deletions(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
index e2b1d536d2..1a5f871d96 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
@@ -124,6 +124,41 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
 }
 
 
+Foam::tmp<Foam::scalarField>
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::K() const
+{
+    const fvMesh& mesh = patch().boundaryMesh().mesh();
+
+    if (mesh.objectRegistry::foundObject<volScalarField>(KName_))
+    {
+        return patch().lookupPatchField<volScalarField, scalar>(KName_);
+    }
+    else if (mesh.objectRegistry::foundObject<volSymmTensorField>(KName_))
+    {
+        const symmTensorField& KWall =
+            patch().lookupPatchField<volSymmTensorField, scalar>(KName_);
+
+        vectorField n = patch().nf();
+
+        return n & KWall & n;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "solidWallHeatFluxTemperatureFvPatchScalarField::K()"
+            " const"
+        )   << "Did not find field " << KName_
+            << " on mesh " << mesh.name() << " patch " << patch().name()
+            << endl
+            << "Please set 'K' to a valid volScalarField"
+            << " or a valid volSymmTensorField." << exit(FatalError);
+
+        return scalarField(0);
+    }
+}
+
+
 void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
 {
     if (updated())
@@ -131,12 +166,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const scalarField& Kw = patch().lookupPatchField<volScalarField, scalar>
-    (
-        KName_
-    );
-
-    gradient() = q_/Kw;
+    gradient() = q_/K();
 
     fixedGradientFvPatchScalarField::updateCoeffs();
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
index ffbcf52ed5..0544b2bffa 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
@@ -31,9 +31,10 @@ Description
         myWallPatch
         {
             type            solidWallHeatFluxTemperature;
-            K               K;                 // Name of K field
-            q               uniform 1000;      // Heat flux / [W/m2]
-            value           300.0;             // Initial temperature / [K]
+            K               K;                  // Name of K field
+            q               uniform 1000;       // Heat flux / [W/m2]
+            value           uniform 300.0;      // Initial temperature / [K]
+            gradient        uniform 0.0;        // Initial gradient / [K/m]
         }
 
 
@@ -140,6 +141,11 @@ public:
 
     // Member functions
 
+        // Helper
+
+            //- Get K field on this patch
+            tmp<scalarField> K() const;
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field

From 1fe76fef6a9fa543ec63d346e6ec7674de522f6b Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 11:35:43 +0100
Subject: [PATCH 07/12] BUG: Load library containing the functionObjects before
 using them.

---
 .../incompressible/pimpleFoam/t-junction/system/controlDict   | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict b/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict
index f1c8b83674..14ca5f3091 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict
@@ -53,9 +53,11 @@ functions
 {
     probes
     {
-        type            probes;
         // Where to load it from
         functionObjectLibs ( "libsampling.so" );
+
+        type            probes;
+
         // Name of the directory for probe data
         name            probes;
 

From 5080bd32370d4e4438bce1b38e9a8b47c957c299 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 11:58:33 +0100
Subject: [PATCH 08/12] STYLE: convert to new time loop structure

---
 .../solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C        | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
index dd05bf88aa..30a5f6fb1a 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
+++ b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

From 21b283779cc9709a82f4ecd3cadc038673d15f7e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 11:59:31 +0100
Subject: [PATCH 09/12] STYLE: new time structure

---
 .../adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index 97992e8192..a18298f826 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

From d4ca02cf7e76e7abd477cd0616899a571520fadc Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 12:22:56 +0100
Subject: [PATCH 10/12] BUG: using Time::run() as a loop construct was not
 doing the functionObjects.

---
 src/OpenFOAM/db/Time/Time.C | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C
index 36bccb3911..35525f65f4 100644
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@@ -508,19 +508,11 @@ bool Foam::Time::run() const
         }
     }
 
-    return running;
-}
-
-
-bool Foam::Time::loop()
-{
-    bool running = run();
-
     if (running)
     {
         if (!subCycling_)
         {
-            readModifiedObjects();
+            const_cast<Time&>(*this).readModifiedObjects();
 
             if (timeIndex_ == startTimeIndex_)
             {
@@ -532,14 +524,22 @@ bool Foam::Time::loop()
             }
         }
 
-        // Check update the "running" status following the "++" operation
-        // to take into account possible side-effects from functionObjects
-        running = run();
+        // Update the "running" status following the
+        // possible side-effects from functionObjects
+        running = value() < (endTime_ - 0.5*deltaT_);
+    }
 
-        if (running)
-        {
-            operator++();
-        }
+    return running;
+}
+
+
+bool Foam::Time::loop()
+{
+    bool running = run();
+
+    if (running)
+    {
+        operator++();
     }
 
     return running;

From 01aa79be546301b1a92fa3804a51d968570736be Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 15:12:48 +0100
Subject: [PATCH 11/12] ENH: multiple faces between two cells gives problems
 for scotch

We filter them out now when constructing the adjacency structure.
---
 .../decompositionMethod/decompositionMethod.C | 249 ++++++++++++------
 .../decompositionMethod/decompositionMethod.H |   4 +
 2 files changed, 178 insertions(+), 75 deletions(-)

diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index f25f0ef448..844272e98b 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -215,6 +215,43 @@ void Foam::decompositionMethod::calcCellCells
 }
 
 
+// Return the minimum face between two cells. Only relevant for
+// cells with multiple faces inbetween.
+Foam::label Foam::decompositionMethod::masterFace
+(
+    const polyMesh& mesh,
+    const label own,
+    const label nei
+)
+{
+    label minFaceI = labelMax;
+
+    // Count multiple faces between own and nei only once
+    const cell& ownFaces = mesh.cells()[own];
+    forAll(ownFaces, i)
+    {
+        label otherFaceI = ownFaces[i];
+
+        if (mesh.isInternalFace(otherFaceI))
+        {
+            label nbrCellI =
+            (
+                mesh.faceNeighbour()[otherFaceI] != own
+              ? mesh.faceNeighbour()[otherFaceI]
+              : mesh.faceOwner()[otherFaceI]
+            );
+
+            if (nbrCellI == nei)
+            {
+                minFaceI = min(minFaceI, otherFaceI);
+            }
+        }
+    }
+
+    return minFaceI;
+}
+
+
 void Foam::decompositionMethod::calcCSR
 (
     const polyMesh& mesh,
@@ -222,61 +259,15 @@ void Foam::decompositionMethod::calcCSR
     List<int>& xadj
 )
 {
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    xadj.setSize(mesh.nCells()+1);
 
-    // Initialise the number of internal faces of the cells to twice the
-    // number of internal faces
-    label nInternalFaces = 2*mesh.nInternalFaces();
-
-    // Check the boundary for coupled patches and add to the number of
-    // internal faces
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            nInternalFaces += pbm[patchi].size();
-        }
-    }
-
-    // Create the adjncy array the size of the total number of internal and
-    // coupled faces
-    adjncy.setSize(nInternalFaces);
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        const labelList& cFaces = mesh.cells()[cellI];
-
-        forAll(cFaces, i)
-        {
-            label faceI = cFaces[i];
-
-            if
-            (
-                mesh.isInternalFace(faceI)
-             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
-            )
-            {
-                freeAdj++;
-            }
-        }
-    }
-    xadj[mesh.nCells()] = freeAdj;
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
+    // Count unique faces between cells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     labelList nFacesPerCell(mesh.nCells(), 0);
 
@@ -286,26 +277,95 @@ void Foam::decompositionMethod::calcCSR
         label own = mesh.faceOwner()[faceI];
         label nei = mesh.faceNeighbour()[faceI];
 
-        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+        if (faceI == masterFace(mesh, own, nei))
+        {
+            nFacesPerCell[own]++;
+            nFacesPerCell[nei]++;
+        }
     }
 
     // Coupled faces. Only cyclics done.
-    forAll(pbm, patchi)
+    HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
+
+    forAll(pbm, patchI)
     {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        if (isA<cyclicPolyPatch>(pbm[patchI]))
         {
-            const unallocLabelList& faceCells = pbm[patchi].faceCells();
+            const unallocLabelList& faceCells = pbm[patchI].faceCells();
 
             label sizeby2 = faceCells.size()/2;
 
-            for (label facei=0; facei<sizeby2; facei++)
+            for (label faceI=0; faceI<sizeby2; faceI++)
             {
-                label own = faceCells[facei];
-                label nei = faceCells[facei + sizeby2];
+                label own = faceCells[faceI];
+                label nei = faceCells[faceI + sizeby2];
 
-                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+                if (cellPair.insert(edge(own, nei)))
+                {
+                    nFacesPerCell[own]++;
+                    nFacesPerCell[nei]++;
+                }
+            }
+        }
+    }
+
+
+    // Size tables
+    // ~~~~~~~~~~~
+
+    // Sum nFacesPerCell
+    xadj.setSize(mesh.nCells()+1);
+
+    label nConnections = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = nConnections;
+        nConnections += nFacesPerCell[cellI];
+    }
+    xadj[mesh.nCells()] = nConnections;
+    adjncy.setSize(nConnections);
+
+
+
+    // Fill tables
+    // ~~~~~~~~~~~
+
+    nFacesPerCell = 0;
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = mesh.faceOwner()[faceI];
+        label nei = mesh.faceNeighbour()[faceI];
+
+        if (faceI == masterFace(mesh, own, nei))
+        {
+            adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+            adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+        }
+    }
+
+    // Coupled faces. Only cyclics done.
+    cellPair.clear();
+    forAll(pbm, patchI)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchI]))
+        {
+            const unallocLabelList& faceCells = pbm[patchI].faceCells();
+
+            label sizeby2 = faceCells.size()/2;
+
+            for (label faceI=0; faceI<sizeby2; faceI++)
+            {
+                label own = faceCells[faceI];
+                label nei = faceCells[faceI + sizeby2];
+
+                if (cellPair.insert(edge(own, nei)))
+                {
+                    adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+                    adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+                }
             }
         }
     }
@@ -320,12 +380,24 @@ void Foam::decompositionMethod::calcCSR
     List<int>& xadj
 )
 {
+    labelHashSet nbrCells;
+
     // Count number of internal faces
     label nConnections = 0;
 
     forAll(cellCells, coarseI)
     {
-        nConnections += cellCells[coarseI].size();
+        nbrCells.clear();
+
+        const labelList& cCells = cellCells[coarseI];
+
+        forAll(cCells, i)
+        {
+            if (nbrCells.insert(cCells[i]))
+            {
+                nConnections++;
+            }
+        }
     }
 
     // Create the adjncy array as twice the size of the total number of
@@ -343,11 +415,16 @@ void Foam::decompositionMethod::calcCSR
     {
         xadj[coarseI] = freeAdj;
 
+        nbrCells.clear();
+
         const labelList& cCells = cellCells[coarseI];
 
         forAll(cCells, i)
         {
-            adjncy[freeAdj++] = cCells[i];
+            if (nbrCells.insert(cCells[i]))
+            {
+                adjncy[freeAdj++] = cCells[i];
+            }
         }
     }
     xadj[cellCells.size()] = freeAdj;
@@ -413,15 +490,21 @@ void Foam::decompositionMethod::calcDistributedCSR
     // Number of faces per cell
     List<int> nFacesPerCell(mesh.nCells(), 0);
 
-    // Number of coupled faces
-    label nCoupledFaces = 0;
-
     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
     {
-        nFacesPerCell[faceOwner[faceI]]++;
-        nFacesPerCell[faceNeighbour[faceI]]++;
+        label own = faceOwner[faceI];
+        label nei = faceNeighbour[faceI];
+
+        if (faceI == masterFace(mesh, own, nei))
+        {
+            nFacesPerCell[own]++;
+            nFacesPerCell[nei]++;
+        }
     }
+
     // Handle coupled faces
+    HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
@@ -429,11 +512,18 @@ void Foam::decompositionMethod::calcDistributedCSR
         if (pp.coupled())
         {
             label faceI = pp.start();
+            label bFaceI = pp.start()-mesh.nInternalFaces();
 
             forAll(pp, i)
             {
-                nCoupledFaces++;
-                nFacesPerCell[faceOwner[faceI++]]++;
+                label own =  faceOwner[faceI];
+                label globalNei = globalNeighbour[bFaceI];
+                if (cellPair.insert(edge(own, globalNei)))
+                {
+                    nFacesPerCell[own]++;
+                }
+                faceI++;
+                bFaceI++;
             }
         }
     }
@@ -459,7 +549,7 @@ void Foam::decompositionMethod::calcDistributedCSR
     // Fill in adjncy
     // ~~~~~~~~~~~~~~
 
-    adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
+    adjncy.setSize(freeAdj);
 
     nFacesPerCell = 0;
 
@@ -469,10 +559,17 @@ void Foam::decompositionMethod::calcDistributedCSR
         label own = faceOwner[faceI];
         label nei = faceNeighbour[faceI];
 
-        adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
+        if (faceI == masterFace(mesh, own, nei))
+        {
+            adjncy[xadj[own] + nFacesPerCell[own]++] =
+                globalCells.toGlobal(nei);
+            adjncy[xadj[nei] + nFacesPerCell[nei]++] =
+                globalCells.toGlobal(own);
+        }
     }
+
     // For boundary faces is offsetted coupled neighbour
+    cellPair.clear();
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
@@ -485,9 +582,11 @@ void Foam::decompositionMethod::calcDistributedCSR
             forAll(pp, i)
             {
                 label own = faceOwner[faceI];
-                adjncy[xadj[own] + nFacesPerCell[own]++] =
-                    globalNeighbour[bFaceI];
-
+                label globalNei = globalNeighbour[bFaceI];
+                if (cellPair.insert(edge(own, globalNei)))
+                {
+                    adjncy[xadj[own] + nFacesPerCell[own]++] = globalNei;
+                }
                 faceI++;
                 bFaceI++;
             }
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
index d729fe9bd3..dabca144c8 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H
@@ -65,6 +65,10 @@ protected:
             labelListList& cellCells
         );
 
+        //- Calculate the minimum face between two neighbouring cells
+        //  (usually there is only one)
+        static label masterFace(const polyMesh&, const label, const label);
+
         // From mesh to compact row storage format
         // (like CompactListList)
         static void calcCSR

From 300ac4843cb1a25b4801d16c9f0443e7da792546 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 15 Apr 2010 18:15:09 +0100
Subject: [PATCH 12/12] BUG: fieldValue.C was defining duplicate IOField
 symbols

---
 .../functionObjects/field/Make/files          |  2 +-
 .../field/fieldValues/fieldValue/fieldValue.C |  5 --
 .../field/streamLine/streamLine.H             |  6 ---
 .../field/streamLine/streamLineParticle.C     | 18 ++-----
 .../field/streamLine/streamLineParticle.H     |  5 --
 .../field/streamLine/vectorIOFieldField.C     | 43 ++++++++++++++++
 .../field/streamLine/vectorIOFieldField.H     | 49 +++++++++++++++++++
 7 files changed, 98 insertions(+), 30 deletions(-)
 create mode 100644 src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.C
 create mode 100644 src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.H

diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files
index bfed281015..df93e62853 100644
--- a/src/postProcessing/functionObjects/field/Make/files
+++ b/src/postProcessing/functionObjects/field/Make/files
@@ -19,6 +19,6 @@ streamLine/streamLine.C
 streamLine/streamLineParticle.C
 streamLine/streamLineParticleCloud.C
 streamLine/streamLineFunctionObject.C
-
+streamLine/vectorIOFieldField.C
 
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
diff --git a/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.C b/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.C
index 1f0f13b184..1bf513e6ea 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.C
@@ -32,11 +32,6 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(fieldValue, 0);
-
-    defineTemplateTypeNameAndDebug(IOField<vector>, 0);
-    defineTemplateTypeNameAndDebug(IOField<sphericalTensor>, 0);
-    defineTemplateTypeNameAndDebug(IOField<symmTensor>, 0);
-    defineTemplateTypeNameAndDebug(IOField<tensor>, 0);
 }
 
 
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.H b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
index fe4f4d287f..0b8dd20ead 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
@@ -206,12 +206,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//#ifdef NoRepository
-//#   include "streamLineTemplates.C"
-//#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
index 9174e35781..ed702f489a 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
@@ -24,22 +24,14 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "streamLineParticle.H"
+#include "vectorIOFieldField.H"
 
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(streamLineParticle, 0);
-
     defineParticleTypeNameAndDebug(streamLineParticle, 0);
-
-    defineTemplateTypeNameAndDebugWithName
-    (
-        IOField<vectorField>,
-        "vectorFieldField",
-        0
-    );
 }
 
 
@@ -402,13 +394,13 @@ void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
     );
     c.checkFieldIOobject(c, lifeTime);
 
-    IOField<pointField> sampledPositions
+    vectorIOFieldField sampledPositions
     (
         c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
     );
     c.checkFieldIOobject(c, sampledPositions);
 
-//    IOField<pointField> sampleVelocity
+//    vectorIOFieldField sampleVelocity
 //    (
 //        c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
 //    );
@@ -436,12 +428,12 @@ void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
         c.fieldIOobject("lifeTime", IOobject::NO_READ),
         np
     );
-    IOField<pointField> sampledPositions
+    vectorIOFieldField sampledPositions
     (
         c.fieldIOobject("sampledPositions", IOobject::NO_READ),
         np
     );
-//    IOField<pointField> sampleVelocity
+//    vectorIOFieldField sampleVelocity
 //    (
 //        c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
 //        np
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
index e2de691282..2ed481c60c 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
@@ -134,11 +134,6 @@ private:
 
 public:
 
-    //- Run-time type information
-    TypeName("streamLineParticle");
-
-
-
     // Constructors
 
         //- Construct from components
diff --git a/src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.C b/src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.C
new file mode 100644
index 0000000000..7ecd4d86a2
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    vectorFieldField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vectorIOFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebugWithName
+    (
+        vectorIOFieldField,
+        "vectorFieldField",
+        0
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.H b/src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.H
new file mode 100644
index 0000000000..b57b7813fd
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/streamLine/vectorIOFieldField.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::vectorIOFieldField
+
+Description
+    vectorField with IO.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef vectorIOFieldField_H
+#define vectorIOFieldField_H
+
+#include "vectorField.H"
+#include "IOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOField<vectorField> vectorIOFieldField;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //