Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-dev-OpenCFD into develop

applications/solvers/basic/potentialFoam/createControls.H

@@ -1,9 +0,0 @@
-const dictionary& potentialFlow
-(
-    mesh.solutionDict().subDict("potentialFlow")
-);
-
-const int nNonOrthCorr
-(
-    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
-);

applications/solvers/basic/potentialFoam/createFields.H

@@ -12,6 +12,7 @@ volVectorField U
     mesh
 );
 
+// Initialise the velocity internal field to zero
 U = dimensionedVector("0", U.dimensions(), vector::zero);
 
 surfaceScalarField phi
@@ -34,7 +35,9 @@ if (args.optionFound("initialiseUBCs"))
 }
 
 
-// Default name for the pressure field
+// Construct a pressure field
+// If it is available read it otherwise construct from the velocity BCs
+// converting fixed-value BCs to zero-gradient and vice versa.
 word pName("p");
 
 // Update name of the pressure field from the command-line option
@@ -55,6 +58,9 @@ forAll(U.boundaryField(), patchi)
     }
 }
 
+// Note that registerObject is false for the pressure field. The pressure
+// field in this solver doesn't have a physical value during the solution.
+// It shouldn't be looked up and used by sub models or boundary conditions.
 Info<< "Constructing pressure field " << pName << nl << endl;
 volScalarField p
 (
@@ -64,13 +70,28 @@ volScalarField p
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
+        IOobject::NO_WRITE,
+        false
     ),
     mesh,
     dimensionedScalar(pName, sqr(dimVelocity), 0),
     pBCTypes
 );
 
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+if (args.optionFound("writep"))
+{
+    setRefCell
+    (
+        p,
+        potentialFlow.dict(),
+        pRefCell,
+        pRefValue
+    );
+}
+
+
 Info<< "Constructing velocity potential field Phi\n" << endl;
 volScalarField Phi
 (

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -24,13 +24,64 @@ License
 Application
     potentialFoam
 
-Description
+Group
-    Potential flow solver which solves for the velocity potential
+    grpBasicSolvers
-    from which the flux-field is obtained and velocity field by reconstructing
-    the flux.
 
-    This application is particularly useful to generate starting fields for
+Description
-    Navier-Stokes codes.
+    Potential flow solver.
+
+    \heading Solver details
+    The potential flow solution is typically employed to generate initial fields
+    for full Navier-Stokes codes.  The flow is evolved using the equation:
+
+    \f[
+        \laplacian \Phi = \div(\vec{U})
+    \f]
+
+    Where:
+    \vartable
+        \Phi      | Velocity potential [m2/s]
+        \vec{U}   | Velocity [m/s]
+    \endvartable
+
+    The corresponding pressure field could be calculated from the divergence
+    of the Euler equation:
+
+    \f[
+        \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+
+    but this generates excessive pressure variation in regions of large
+    velocity gradient normal to the flow direction.  A better option is to
+    calculate the pressure field corresponding to velocity variation along the
+    stream-lines:
+
+    \f[
+        \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+    where the flow direction tensor \f$\vec{F}\f$ is obtained from
+    \f[
+        \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
+    \f]
+
+    \heading Required fields
+    \plaintable
+        U         | Velocity [m/s]
+    \endplaintable
+
+    \heading Optional fields
+    \plaintable
+        p         | Kinematic pressure [m2/s2]
+        Phi       | Velocity potential [m2/s]
+                  | Generated from p (if present) or U if not present
+    \endplaintable
+
+    \heading Options
+    \plaintable
+        -writep   | write the Euler pressure
+        -writePhi | Write the final velocity potential
+        -initialiseUBCs | Update the velocity boundaries before solving for Phi
+    \endplaintable
 
 \*---------------------------------------------------------------------------*/
 
@@ -58,19 +109,13 @@ int main(int argc, char *argv[])
     argList::addBoolOption
     (
         "writePhi",
-        "Write the velocity potential field"
+        "Write the final velocity potential field"
     );
 
     argList::addBoolOption
     (
         "writep",
-        "Calculate and write the pressure field"
+        "Calculate and write the Euler pressure field"
-    );
-
-    argList::addBoolOption
-    (
-        "withFunctionObjects",
-        "execute functionObjects"
     );
 
     #include "setRootCase.H"
@@ -137,7 +182,7 @@ int main(int argc, char *argv[])
         Phi.write();
     }
 
-    // Calculate the pressure field
+    // Calculate the pressure field from the Euler equation
     if (args.optionFound("writep"))
     {
         Info<< nl << "Calculating approximate pressure field" << endl;

applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C

@@ -32,6 +32,7 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "dynamicFvMesh.H"
+#include "pimpleControl.H"
 #include "vtkSurfaceWriter.H"
 #include "cyclicAMIPolyPatch.H"
 
@@ -129,14 +130,24 @@ int main(int argc, char *argv[])
         Info<< "Writing VTK files with weights of AMI patches." << nl << endl;
     }
 
+    pimpleControl pimple(mesh);
+
+    bool moveMeshOuterCorrectors
+    (
+        pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
+    );
+
     while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << endl;
 
-        for (int i = 0; i<2; i++)
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
             {
                 mesh.update();
             }
+        }
 
         mesh.checkMesh(true);
 

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -31,18 +31,66 @@ setFormat raw;
 //      dx          : DX scalar or vector format
 //      vtk         : VTK ascii format
 //      raw         : x y z value format for use with e.g. gnuplot 'splot'.
+//      boundaryData: written in a form that can be used with timeVaryingMapped
+//                    boundary conditions (move to constant/boundaryData)
+//      starcd      : Star-CD format (geometry in .cel/.vrt, data in .usr)
+//      nastran     : Nastran (.nas) format
 //
 // Note:
-// other formats such as obj, stl, etc can also be written (by proxy)
+// - other formats such as obj, stl, etc can also be written (by proxy)
 //   but without any values!
+// - boundaryData format: typical use:
+//      - use a surfaces functionObject to sample a patch:
+//            type            surfaces;
+//            surfaceFormat   boundaryData;
+//            fields          ( p );
+//            surfaces
+//            (
+//                inlet
+//                {
+//                    type            patch;
+//                    patches         (inpletPatch);
+//                    interpolate     false;
+//                }
+//            );
+//      - the boundaryData writer will write postProcessing/surfaces/inlet/
+//      - move this to constant/boundaryData/inlet of the destination case
+//      - use a timeVaryingMappedFixedValue bc to read&interpolate and use
+//        as fixedValue.
+// - boundaryData: 'interpolate false' writes face centres, 'interpolate true'
+//   writes points. For 2D case the face centres might only be a single column
+//   so cannot be used in timeVaryingMapped with standard mapping (since
+//   uses triangulation to do interpolation).
+
 surfaceFormat vtk;
 
-// optionally define extra controls for the output formats
+// Optionally define extra controls for the output formats
 formatOptions
 {
     ensight
     {
+        // ascii/binary format
         format  ascii;
+        //collateTimes true;  // write single file containing multiple timesteps
+                              // (only for static surfaces)
+    }
+    nastran
+    {
+        // From OpenFOAM field name to Nastran field name
+        fields
+        (
+            (U      PLOAD4)
+            (p      PLOAD2)
+        );
+        // Optional scale
+        scale 1.0;
+        // Optional format
+        // short/long/free format
+        // short: scalar field width 8. default.
+        // long : scalar field width 16
+        // free : comma separated, field width according to writePrecision
+        //        setting
+        format  short;      // short/long/free
     }
 }
 
@@ -143,9 +191,19 @@ sets
         type        patchSeed;
         axis        xyz;
         patches     (".*Wall.*");
-        // Number of points to seed. Divided amongst all processors according
+
-        // to fraction of patches they hold.
+        // Subset patch faces by:
+
+        // 1. Number of points to seed. Divided amongst all processors
+        //   according to fraction of patches they hold.
         maxPoints   100;
+
+        // 2. Specified set of locations. This selects for every the specified
+        //    point the nearest patch face. (in addition the number of points
+        //    is also truncated by the maxPoints setting)
+        //    The difference with patchCloud is that this selects patch
+        //    face centres, not an arbitrary location on the face.
+        points   ((0.049 0.099 0.005)(0.051 0.054 0.005));
     }
 
 );
@@ -252,6 +310,8 @@ surfaces
         // zone            ABC;         // Optional: zone only
         // exposedPatchName fixedWalls; // Optional: zone only
 
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
+
         // regularise      false;    // Optional: do not simplify
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
                                      // to merge points (default=1e-6)
@@ -265,6 +325,9 @@ surfaces
         isoValue        0.5;
         interpolate     false;
         regularise      false;              // do not simplify
+
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
+
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
                                      // to merge points (default=1e-6)
     }
@@ -284,6 +347,8 @@ surfaces
         // zone             ABC;          // Optional: zone only
         // exposedPatchName fixedWalls;  // Optional: zone only
 
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
+
         // regularise      false;    // Optional: do not simplify
         // mergeTol        1e-10;    // Optional: fraction of mesh bounding box
                                      // to merge points (default=1e-6)
@@ -305,6 +370,8 @@ surfaces
                                 // of isoSurfaceCell
         interpolate     false;
         regularise      false;  // Optional: do not simplify
+        // bounds       (1 1 1)(2 2 2); // Optional: limit extent
+
         // mergeTol 1e-10;      // Optional: fraction of mesh bounding box
                                 // to merge points (default=1e-6)
     }

applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H

@@ -0,0 +1,10 @@
+//
+// addRegionOption.H
+// ~~~~~~~~~~~~~~~~~
+
+    Foam::argList::addOption
+    (
+        "regions",
+        "(name1 .. nameN)",
+        "specify alternative mesh regions"
+    );

applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C

@@ -38,6 +38,11 @@ Usage
     \param -region \<name\> \n
     Specify an alternative mesh region.
 
+    \param -regions (\<name1\> \<name2\> .. \<namen\>) \n
+    Specify alternative mesh regions. The region names will be sorted
+    alphabetically and a single composite name will be created
+        \<nameX\>_\<nameY\>.._\<nameZ\>
+
     On execution, the combined patch geometry (points and faces) are output
     to the communications directory.
 
@@ -59,6 +64,7 @@ SeeAlso
 int main(int argc, char *argv[])
 {
     #include "addRegionOption.H"
+    #include "addRegionsOption.H"
     argList::validArgs.append("patchGroup");
     argList::addOption
     (
@@ -68,16 +74,52 @@ int main(int argc, char *argv[])
     );
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createNamedMesh.H"
 
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    wordList regionNames(1, fvMesh::defaultRegion);
+    if (!args.optionReadIfPresent("region", regionNames[0]))
+    {
+        args.optionReadIfPresent("regions", regionNames);
+    }
 
-    const wordRe patchGroup(args[1]);
+    const wordRe patchGroup(args.argRead<wordRe>(1));
 
     fileName commsDir(runTime.path()/"comms");
     args.optionReadIfPresent("commsDir", commsDir);
 
-    externalCoupledFunctionObject::writeGeometry(mesh, commsDir, patchGroup);
+
+    // Make sure region names are in canonical order
+    stableSort(regionNames);
+
+
+    PtrList<const fvMesh> meshes(regionNames.size());
+    forAll(regionNames, i)
+    {
+        Info<< "Create mesh " << regionNames[i] << " for time = "
+            << runTime.timeName() << nl << endl;
+
+        meshes.set
+        (
+            i,
+            new fvMesh
+            (
+                Foam::IOobject
+                (
+                    regionNames[i],
+                    runTime.timeName(),
+                    runTime,
+                    Foam::IOobject::MUST_READ
+                )
+            )
+        );
+    }
+
+
+    externalCoupledFunctionObject::writeGeometry
+    (
+        UPtrList<const fvMesh>(meshes),
+        commsDir,
+        patchGroup
+    );
 
     Info<< "\nEnd\n" << endl;
 

etc/caseDicts/setConstraintTypes

@@ -39,11 +39,13 @@ nonuniformTransformCyclic
 processor
 {
     type  processor;
+    value $internalField;
 }
 
 processorCyclic
 {
     type  processorCyclic;
+    value $internalField;
 }
 
 symmetryPlane

src/OpenFOAM/db/error/IOerror.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -174,7 +174,7 @@ void Foam::IOerror::exit(const int)
         jobInfo.exit();
     }
 
-    if (abort_)
+    if (env("FOAM_ABORT"))
     {
         abort();
     }
@@ -215,7 +215,7 @@ void Foam::IOerror::abort()
         jobInfo.abort();
     }
 
-    if (abort_)
+    if (env("FOAM_ABORT"))
     {
         Perr<< endl << *this << endl
             << "\nFOAM aborting (FOAM_ABORT set)\n" << endl;

src/OpenFOAM/db/error/error.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,7 +40,6 @@ Foam::error::error(const string& title)
     functionName_("unknown"),
     sourceFileName_("unknown"),
     sourceFileLineNumber_(0),
-    abort_(env("FOAM_ABORT")),
     throwExceptions_(false),
     messageStreamPtr_(new OStringStream())
 {
@@ -61,7 +60,6 @@ Foam::error::error(const dictionary& errDict)
     functionName_(errDict.lookup("functionName")),
     sourceFileName_(errDict.lookup("sourceFileName")),
     sourceFileLineNumber_(readLabel(errDict.lookup("sourceFileLineNumber"))),
-    abort_(env("FOAM_ABORT")),
     throwExceptions_(false),
     messageStreamPtr_(new OStringStream())
 {
@@ -83,7 +81,6 @@ Foam::error::error(const error& err)
     functionName_(err.functionName_),
     sourceFileName_(err.sourceFileName_),
     sourceFileLineNumber_(err.sourceFileLineNumber_),
-    abort_(err.abort_),
     throwExceptions_(err.throwExceptions_),
     messageStreamPtr_(new OStringStream(*err.messageStreamPtr_))
 {
@@ -173,7 +170,7 @@ void Foam::error::exit(const int errNo)
         jobInfo.exit();
     }
 
-    if (abort_)
+    if (env("FOAM_ABORT"))
     {
         abort();
     }
@@ -214,7 +211,7 @@ void Foam::error::abort()
         jobInfo.abort();
     }
 
-    if (abort_)
+    if (env("FOAM_ABORT"))
     {
         Perr<< endl << *this << endl
             << "\nFOAM aborting (FOAM_ABORT set)\n" << endl;

src/OpenFOAM/db/error/error.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,11 +78,10 @@ protected:
         string sourceFileName_;
         label sourceFileLineNumber_;
 
-        bool abort_;
-
         bool throwExceptions_;
         OStringStream* messageStreamPtr_;
 
+
 public:
 
     // Constructors

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C

@@ -0,0 +1,242 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triangle.H"
+#include "triPoints.H"
+#include "plane.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+template<class AboveOp, class BelowOp>
+inline void Foam::triangle<Point, PointRef>::triSliceWithPlane
+(
+    const plane& pl,
+    const triPoints& tri,
+    AboveOp& aboveOp,
+    BelowOp& belowOp
+)
+{
+    // distance to cutting plane
+    FixedList<scalar, 3> d;
+
+    // determine how many of the points are above the cutting plane
+    label nPos = 0;
+    label posI = -1;
+    label negI = -1;
+    forAll(tri, i)
+    {
+        d[i] = ((tri[i] - pl.refPoint()) & pl.normal());
+
+        if (d[i] > 0)
+        {
+            nPos++;
+            posI = i;
+        }
+        else
+        {
+            negI = i;
+        }
+    }
+
+    if (nPos == 3)
+    {
+        aboveOp(tri);
+    }
+    else if (nPos == 2)
+    {
+        // point under the plane
+        label i0 = negI;
+
+        // indices of remaining points
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+
+        // determine the two intersection points
+        point p01 = planeIntersection(d, tri, i0, i1);
+        point p02 = planeIntersection(d, tri, i0, i2);
+
+        aboveOp(triPoints(tri[i1], tri[i2], p02));
+        aboveOp(triPoints(tri[i1], p02, p01));
+        belowOp(triPoints(tri[i0], p01, p02));
+    }
+    else if (nPos == 1)
+    {
+        // point above the plane
+        label i0 = posI;
+
+        // indices of remaining points
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+
+        // determine the two intersection points
+        point p01 = planeIntersection(d, tri, i0, i1);
+        point p02 = planeIntersection(d, tri, i0, i2);
+
+        belowOp(triPoints(tri[i1], tri[i2], p02));
+        belowOp(triPoints(tri[i1], p02, p01));
+        aboveOp(triPoints(tri[i0], p01, p02));
+    }
+    else
+    {
+        // All below
+        belowOp(tri);
+    }
+}
+
+
+template<class Point, class PointRef>
+template<class AboveOp, class BelowOp>
+inline void Foam::triangle<Point, PointRef>::sliceWithPlane
+(
+    const plane& pl,
+    AboveOp& aboveOp,
+    BelowOp& belowOp
+) const
+{
+    triSliceWithPlane(pl, triPoints(a_, b_, c_), aboveOp, belowOp);
+}
+
+
+template<class Point, class PointRef>
+template<class InsideOp, class OutsideOp>
+inline void Foam::triangle<Point, PointRef>::triangleOverlap
+(
+    const vector& n,
+    const triangle<Point, PointRef>& tgt,
+    InsideOp& insideOp,
+    OutsideOp& outsideOp
+) const
+{
+    // There are two possibilities with this algorithm - we either cut
+    // the outside triangles with all the edges or not (and keep them
+    // as disconnected triangles). In the first case
+    // we cannot do any evaluation short cut so we've chosen not to re-cut
+    // the outside triangles.
+
+
+    triIntersectionList insideTrisA;
+    label nInsideA = 0;
+    storeOp insideOpA(insideTrisA, nInsideA);
+
+    triIntersectionList outsideTrisA;
+    label nOutsideA = 0;
+    storeOp outsideOpA(outsideTrisA, nOutsideA);
+
+
+    const triPoints thisTri(a_, b_, c_);
+
+
+    // Cut original triangle with tgt edge 0.
+    // From *this to insideTrisA, outsideTrisA.
+    {
+        scalar s = Foam::mag(tgt.b() - tgt.a());
+        const plane pl0(tgt.a(), tgt.b(), tgt.b() + s*n);
+        triSliceWithPlane(pl0, thisTri, insideOpA, outsideOpA);
+    }
+
+    // Shortcut if nothing cut
+
+    if (insideOpA.nTris_ == 0)
+    {
+        outsideOp(thisTri);
+        return;
+    }
+
+    if (outsideOpA.nTris_ == 0)
+    {
+        insideOp(thisTri);
+        return;
+    }
+
+
+    // Cut all triangles with edge 1.
+    // From insideTrisA to insideTrisB, outsideTrisA
+
+    triIntersectionList insideTrisB;
+    label nInsideB = 0;
+    storeOp insideOpB(insideTrisB, nInsideB);
+
+    //triIntersectionList outsideTrisB;
+    //label nOutsideB = 0;
+    //storeOp outsideOpB(outsideTrisB, nOutsideB);
+
+    {
+        scalar s = Foam::mag(tgt.c() - tgt.b());
+        const plane pl0(tgt.b(), tgt.c(), tgt.c() + s*n);
+
+        for (label i = 0; i < insideOpA.nTris_; i++)
+        {
+            const triPoints& tri = insideOpA.tris_[i];
+            triSliceWithPlane(pl0, tri, insideOpB, outsideOpA);
+        }
+
+        //// Recut outside triangles (not necessary if only interested in
+        //// intersection properties)
+        //for (label i = 0; i < outsideOpA.nTris_; i++)
+        //{
+        //    const triPoints& tri = outsideOpA.tris_[i];
+        //    triSliceWithPlane(pl0, tri, outsideOpB, outsideOpB);
+        //}
+    }
+
+
+    // Cut all triangles with edge 2.
+    // From insideTrisB to insideTrisA, outsideTrisA
+    {
+        scalar s = Foam::mag(tgt.a() - tgt.c());
+        const plane pl0(tgt.c(), tgt.a(), tgt.a() + s*n);
+
+        insideOpA.nTris_ = 0;
+        //outsideOpA.nTris_ = 0;
+        for (label i = 0; i < insideOpB.nTris_; i++)
+        {
+            const triPoints& tri = insideOpB.tris_[i];
+            triSliceWithPlane(pl0, tri, insideOpA, outsideOpA);
+        }
+
+        //// Recut outside triangles (not necessary if only interested in
+        //// intersection properties)
+        //for (label i = 0; i < outsideOpB.nTris_; i++)
+        //{
+        //    const triPoints& tri = outsideOpB.tris_[i];
+        //    triSliceWithPlane(pl0, tri, outsideOpA, outsideOpA);
+        //}
+    }
+
+    // Transfer from A to argument
+    for (label i = 0; i < insideOpA.nTris_; i++)
+    {
+        insideOp(insideOpA.tris_[i]);
+    }
+
+    for (label i = 0; i < outsideOpA.nTris_; i++)
+    {
+        outsideOp(outsideOpA.tris_[i]);
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,6 +52,8 @@ namespace Foam
 
 class Istream;
 class Ostream;
+class triPoints;
+class plane;
 
 // Forward declaration of friend functions and operators
 
@@ -71,6 +73,8 @@ inline Ostream& operator<<
     const triangle<Point, PointRef>&
 );
 
+typedef triangle<point, const point&> triPointRef;
+
 
 /*---------------------------------------------------------------------------*\
                           Class triangle Declaration
@@ -79,13 +83,14 @@ inline Ostream& operator<<
 template<class Point, class PointRef>
 class triangle
 {
-    // Private data
-
-        PointRef a_, b_, c_;
-
-
 public:
 
+    // Public typedefs
+
+        //- Storage type for triangles originating from intersecting triangle
+        //  with another triangle
+        typedef FixedList<triPoints, 27> triIntersectionList;
+
         //- Return types for classify
         enum proxType
         {
@@ -95,11 +100,81 @@ public:
         };
 
 
+    // Public classes
+
+        // Classes for use in sliceWithPlane. What to do with decomposition
+        // of triangle.
+
+            //- Dummy
+            class dummyOp
+            {
+            public:
+                inline void operator()(const triPoints&);
+            };
+
+            //- Sum resulting areas
+            class sumAreaOp
+            {
+            public:
+                scalar area_;
+
+                inline sumAreaOp();
+
+                inline void operator()(const triPoints&);
+            };
+
+            //- Store resulting tris
+            class storeOp
+            {
+            public:
+                triIntersectionList& tris_;
+                label& nTris_;
+
+                inline storeOp(triIntersectionList&, label&);
+
+                inline void operator()(const triPoints&);
+            };
+
+
+private:
+
+    // Private data
+
+        PointRef a_, b_, c_;
+
+
+   // Private Member Functions
+
+        //- Helper: calculate intersection point
+        inline static point planeIntersection
+        (
+            const FixedList<scalar, 3>& d,
+            const triPoints& t,
+            const label negI,
+            const label posI
+        );
+
+        //- Helper: slice triangle with plane
+        template<class AboveOp, class BelowOp>
+        inline static void triSliceWithPlane
+        (
+            const plane& pl,
+            const triPoints& tri,
+            AboveOp& aboveOp,
+            BelowOp& belowOp
+        );
+
+
+public:
+
     // Constructors
 
         //- Construct from three points
         inline triangle(const Point& a, const Point& b, const Point& c);
 
+        //- Construct from three points
+        inline triangle(const FixedList<Point, 3>&);
+
         //- Construct from three points in the list of points
         //  The indices could be from triFace etc.
         inline triangle
@@ -237,6 +312,27 @@ public:
                 pointHit& edgePoint
             ) const;
 
+            //- Decompose triangle into triangles above and below plane
+            template<class AboveOp, class BelowOp>
+            inline void sliceWithPlane
+            (
+                const plane& pl,
+                AboveOp& aboveOp,
+                BelowOp& belowOp
+            ) const;
+
+            //- Decompose triangle into triangles inside and outside
+            //  (with respect to user provided normal) other
+            //  triangle.
+            template<class InsideOp, class OutsideOp>
+            inline void triangleOverlap
+            (
+                const vector& n,
+                const triangle<Point, PointRef>& tri,
+                InsideOp& insideOp,
+                OutsideOp& outsideOp
+            ) const;
+
 
     // IOstream operators
 
@@ -264,6 +360,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "triangle.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "IOstreams.H"
 #include "pointHit.H"
+#include "triPoints.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -43,6 +44,18 @@ inline Foam::triangle<Point, PointRef>::triangle
 {}
 
 
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::triangle
+(
+    const FixedList<Point, 3>& tri
+)
+:
+    a_(tri[0]),
+    b_(tri[1]),
+    c_(tri[2])
+{}
+
+
 template<class Point, class PointRef>
 inline Foam::triangle<Point, PointRef>::triangle
 (
@@ -804,6 +817,66 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
 }
 
 
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::dummyOp::operator()
+(
+    const triPoints&
+)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::sumAreaOp::sumAreaOp()
+:
+    area_(0.0)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::sumAreaOp::operator()
+(
+    const triPoints& tri
+)
+{
+    area_ += triangle<Point, const Point&>(tri).mag();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::storeOp::storeOp
+(
+    triIntersectionList& tris,
+    label& nTris
+)
+:
+    tris_(tris),
+    nTris_(nTris)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::storeOp::operator()
+(
+    const triPoints& tri
+)
+{
+    tris_[nTris_++] = tri;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::point Foam::triangle<Point, PointRef>::planeIntersection
+(
+    const FixedList<scalar, 3>& d,
+    const triPoints& t,
+    const label negI,
+    const label posI
+)
+{
+    return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI] + d[posI]);
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>

src/dynamicMesh/motionSmoother/motionSmootherAlgo.C

@@ -514,12 +514,31 @@ void Foam::motionSmootherAlgo::setDisplacement
     const labelList& cppMeshPoints =
         mesh.globalData().coupledPatch().meshPoints();
 
+    const labelList& ppMeshPoints = pp.meshPoints();
+
+    // Knock out displacement on points which are not on pp but are coupled
+    // to them since we want 'proper' values from displacement to take
+    // precedence.
+    {
+        PackedBoolList isPatchPoint(mesh.nPoints());
+        isPatchPoint.set(ppMeshPoints);
+        syncTools::syncPointList
+        (
+            mesh,
+            isPatchPoint,
+            maxEqOp<unsigned int>(),
+            0
+        );
         forAll(cppMeshPoints, i)
         {
-        displacement[cppMeshPoints[i]] = vector::zero;
+            label pointI = cppMeshPoints[i];
+            if (isPatchPoint[pointI])
+            {
+                displacement[pointI] = vector::zero;
+            }
+        }
     }
 
-    const labelList& ppMeshPoints = pp.meshPoints();
 
     // Set internal point data from displacement on combined patch points.
     forAll(ppMeshPoints, patchPointI)

src/dynamicMesh/motionSmoother/motionSmootherAlgo.H

@@ -56,7 +56,8 @@ Note
     and/or edges but no faces of pp). Hence we have to be careful when e.g.
     synchronising displacements that the value from the processor which has
     faces of pp get priority. This is currently handled in setDisplacement
-    by resetting the internal displacement to zero before doing anything
+    by resetting the internal displacement to zero on coupled points
+    that are coupled to patch points before doing anything
     else. The combine operator used will give preference to non-zero
     values.
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,12 +48,9 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
         return 1.0;
     }
 
-    //scalar totalOverFirst = totalThickness/firstLayerThickess;
+    const label maxIters = 20;
-
-    const label maxIters = 10;
     const scalar tol = 1e-8;
 
-
     if (mag(n-totalOverFirst) < tol)
     {
         return 1.0;
@@ -74,8 +71,6 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
         maxR = totalOverFirst/(n - 1);
     }
 
-    //Info<< "Solution bounds = (" << minR << ", " << maxR << ")" << nl << endl;
-
     // Starting guess
     scalar r = 0.5*(minR + maxR);
 
@@ -85,14 +80,9 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 
         const scalar fx = pow(r, n) - totalOverFirst*r - (1 - totalOverFirst);
         const scalar dfx = n*pow(r, n - 1) - totalOverFirst;
-
         r -= fx/dfx;
 
-        const scalar error = mag(r - prevr);
+        if (mag(r - prevr) < tol)
-
-        //Info<< i << " " << r << " Error = " << error << endl;
-
-        if (error < tol)
         {
             break;
         }
@@ -103,7 +93,6 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from dictionary
 Foam::layerParameters::layerParameters
 (
     const dictionary& dict,
@@ -417,8 +406,8 @@ Foam::scalar Foam::layerParameters::layerThickness
             }
             else
             {
-                return firstLayerThickess *
+                return firstLayerThickess
-                    (1.0 - pow(expansionRatio, nLayers))
+                   *(1.0 - pow(expansionRatio, nLayers))
                    /(1.0 - expansionRatio);
             }
         }
@@ -433,8 +422,8 @@ Foam::scalar Foam::layerParameters::layerThickness
             else
             {
                 scalar invExpansion = 1.0 / expansionRatio;
-                return finalLayerThickess *
+                return finalLayerThickess
-                    (1.0 - pow(invExpansion, nLayers))
+                   *(1.0 - pow(invExpansion, nLayers))
                    /(1.0 - invExpansion);
             }
         }

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -333,7 +333,6 @@ public:
                 const label nLayers,
                 const scalar expansionRatio
             ) const;
-
 };
 
 

src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H

@@ -33,8 +33,8 @@ SeeAlso
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef plane_H
+#ifndef planeExtrusion_H
-#define plane_H
+#define planeExtrusion_H
 
 #include "linearNormal.H"
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C

@@ -44,7 +44,7 @@ const Foam::scalar Foam::cyclicACMIPolyPatch::tolerance_ = 1e-6;
 
 void Foam::cyclicACMIPolyPatch::initPatchFaceAreas() const
 {
-    if (!empty() && faceAreas0_.empty())
+    if (!empty() && (faceAreas0_.empty() || boundaryMesh().mesh().moving()))
     {
         faceAreas0_ = faceAreas();
     }
@@ -52,9 +52,13 @@ void Foam::cyclicACMIPolyPatch::initPatchFaceAreas() const
     const cyclicACMIPolyPatch& nbrACMI =
         refCast<const cyclicACMIPolyPatch>(this->neighbPatch());
 
-    if (!nbrACMI.empty() && nbrACMI.faceAreas0().empty())
+    if
+    (
+        !nbrACMI.empty()
+     && (nbrACMI.faceAreas0().empty() || boundaryMesh().mesh().moving())
+    )
     {
-        nbrACMI.initPatchFaceAreas();
+        nbrACMI.faceAreas0_ = nbrACMI.faceAreas();
     }
 }
 
@@ -136,11 +140,13 @@ void Foam::cyclicACMIPolyPatch::setNeighbourFaceAreas() const
 
 void Foam::cyclicACMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
-    // Initialise the AMI so that base geometry (e.g. cell volumes) are
+    // Note: cyclicAMIPolyPatch clears AMI so do first
-    // correctly evaluated
-    resetAMI();
-
     cyclicAMIPolyPatch::initGeometry(pBufs);
+
+    // Initialise the AMI so that base geometry (e.g. cell volumes) are
+    // correctly evaluated before e.g. any of the processor patches gets
+    // hit (since uses cell volumes in its initGeometry)
+    resetAMI();
 }
 
 
@@ -156,7 +162,13 @@ void Foam::cyclicACMIPolyPatch::initMovePoints
     const pointField& p
 )
 {
+    // Note: cyclicAMIPolyPatch clears AMI so do first
     cyclicAMIPolyPatch::initMovePoints(pBufs, p);
+
+    // Initialise the AMI so that base geometry (e.g. cell volumes) are
+    // correctly evaluated before e.g. any of the processor patches gets
+    // hit (since uses cell volumes in its initGeometry)
+    resetAMI();
 }
 
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
 
     if (owner())
     {
-        const scalarField& w = srcMask_;
+        const scalarField& w = AMI().srcWeightsSum();
 
         tmp<Field<Type> > interpField(AMI().interpolateToSource(fldCouple));
 
@@ -45,7 +45,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIPolyPatch::interpolate
     }
     else
     {
-        const scalarField& w = neighbPatch().tgtMask();
+        const scalarField& w = neighbPatch().AMI().tgtWeightsSum();
 
         tmp<Field<Type> > interpField
         (
@@ -82,16 +82,18 @@ void Foam::cyclicACMIPolyPatch::interpolate
 
     if (owner())
     {
-        const scalarField& w = srcMask_;
+        const scalarField& w = AMI().srcWeightsSum();
 
         AMI().interpolateToSource(fldCouple, cop, result);
+
         result = result + (1.0 - w)*fldNonOverlap;
     }
     else
     {
-        const scalarField& w = neighbPatch().tgtMask();
+        const scalarField& w = neighbPatch().AMI().tgtWeightsSum();
 
         neighbPatch().AMI().interpolateToTarget(fldCouple, cop, result);
+
         result = result + (1.0 - w)*fldNonOverlap;
     }
 }

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -24,9 +24,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cyclicAMIPolyPatch.H"
+#include "transformField.H"
 #include "SubField.H"
+#include "polyMesh.H"
 #include "Time.H"
 #include "addToRunTimeSelectionTable.H"
+#include "faceAreaIntersect.H"
+#include "ops.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -172,18 +176,19 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
                 scalar errorPos = mag(transformedAreaPos + area0);
                 scalar errorNeg = mag(transformedAreaNeg + area0);
 
-                if (errorPos < errorNeg)
+                scalar scaledErrorPos = errorPos/(mag(area0) + ROOTVSMALL);
-                {
+                scalar scaledErrorNeg = errorNeg/(mag(area0) + ROOTVSMALL);
+
+                // One of the errors should be (close to) zero. If this is
+                // the reverse transformation flip the rotation angle.
                 revT = revTPos;
-                }
+                if (errorPos > errorNeg && scaledErrorNeg <= matchTolerance())
-                else
                 {
                     revT = revTNeg;
                     rotationAngle_ *= -1;
                 }
 
-                scalar areaError =
+                scalar areaError = min(scaledErrorPos, scaledErrorNeg);
-                    min(errorPos, errorNeg)/(mag(area0) + ROOTVSMALL);
 
                 if (areaError > matchTolerance())
                 {
@@ -406,6 +411,9 @@ void Foam::cyclicAMIPolyPatch::resetAMI
 
 void Foam::cyclicAMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
+    // The AMI is no longer valid. Leave it up to demand-driven calculation
+    AMIPtr_.clear();
+
     polyPatch::initGeometry(pBufs);
 }
 
@@ -431,6 +439,9 @@ void Foam::cyclicAMIPolyPatch::initMovePoints
     const pointField& p
 )
 {
+    // The AMI is no longer valid. Leave it up to demand-driven calculation
+    AMIPtr_.clear();
+
     polyPatch::initMovePoints(pBufs, p);
 
     // See below. Clear out any local geometry
@@ -447,19 +458,15 @@ void Foam::cyclicAMIPolyPatch::movePoints
     polyPatch::movePoints(pBufs, p);
 
     calcTransforms();
-
-    // Note: resetAMI is called whilst in geometry update. So the slave
-    // side might not have reached 'movePoints'. Is explicitly handled by
-    // - clearing geometry of neighbour inside initMovePoints
-    // - not using localPoints() inside resetAMI
-    resetAMI();
 }
 
 
 void Foam::cyclicAMIPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    polyPatch::initUpdateMesh(pBufs);
+    // The AMI is no longer valid. Leave it up to demand-driven calculation
     AMIPtr_.clear();
+
+    polyPatch::initUpdateMesh(pBufs);
 }
 
 

src/meshTools/regionSplit/regionSplit.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,8 @@ License
 #include "processorPolyPatch.H"
 #include "globalIndex.H"
 #include "syncTools.H"
+#include "FaceCellWave.H"
+#include "minData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,342 +42,95 @@ defineTypeNameAndDebug(regionSplit, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Handle (non-processor) coupled faces.
+void Foam::regionSplit::calcNonCompactRegionSplit
-void Foam::regionSplit::transferCoupledFaceRegion
-(
-    const label faceI,
-    const label otherFaceI,
-
-    labelList& faceRegion,
-    DynamicList<label>& newChangedFaces
-) const
-{
-    if (faceRegion[faceI] >= 0)
-    {
-        if (faceRegion[otherFaceI] == -1)
-        {
-            faceRegion[otherFaceI] = faceRegion[faceI];
-            newChangedFaces.append(otherFaceI);
-        }
-        else if (faceRegion[otherFaceI] == -2)
-        {
-            // otherFaceI blocked but faceI is not. Is illegal for coupled
-            // faces, not for explicit connections.
-        }
-        else if (faceRegion[otherFaceI] != faceRegion[faceI])
-        {
-            FatalErrorIn
-            (
-                  "regionSplit::transferCoupledFaceRegion"
-                  "(const label, const label, labelList&, labelList&) const"
-              )   << "Problem : coupled face " << faceI
-                  << " on patch " << mesh().boundaryMesh().whichPatch(faceI)
-                  << " has region " << faceRegion[faceI]
-                  << " but coupled face " << otherFaceI
-                  << " has region " << faceRegion[otherFaceI]
-                  << endl
-                  << "Is your blocked faces specification"
-                  << " synchronized across coupled boundaries?"
-                  << abort(FatalError);
-        }
-    }
-    else if (faceRegion[faceI] == -1)
-    {
-        if (faceRegion[otherFaceI] >= 0)
-        {
-            faceRegion[faceI] = faceRegion[otherFaceI];
-            newChangedFaces.append(faceI);
-        }
-        else if (faceRegion[otherFaceI] == -2)
-        {
-            // otherFaceI blocked but faceI is not. Is illegal for coupled
-            // faces, not for explicit connections.
-        }
-    }
-}
-
-
-void Foam::regionSplit::fillSeedMask
-(
-    const List<labelPair>& explicitConnections,
-    labelList& cellRegion,
-    labelList& faceRegion,
-    const label seedCellID,
-    const label markValue
-) const
-{
-    // Do seed cell
-    cellRegion[seedCellID] = markValue;
-
-
-    // Collect faces on seed cell
-    const cell& cFaces = mesh().cells()[seedCellID];
-
-    label nFaces = 0;
-
-    labelList changedFaces(cFaces.size());
-
-    forAll(cFaces, i)
-    {
-        label faceI = cFaces[i];
-
-        if (faceRegion[faceI] == -1)
-        {
-            faceRegion[faceI] = markValue;
-            changedFaces[nFaces++] = faceI;
-        }
-    }
-    changedFaces.setSize(nFaces);
-
-
-    // Loop over changed faces. FaceCellWave in small.
-
-    while (changedFaces.size())
-    {
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedFaces:"
-        //        << changedFaces.size() << endl;
-        //}
-
-        DynamicList<label> changedCells(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label faceI = changedFaces[i];
-
-            label own = mesh().faceOwner()[faceI];
-
-            if (cellRegion[own] == -1)
-            {
-                cellRegion[own] = markValue;
-                changedCells.append(own);
-            }
-
-            if (mesh().isInternalFace(faceI))
-            {
-                label nei = mesh().faceNeighbour()[faceI];
-
-                if (cellRegion[nei] == -1)
-                {
-                    cellRegion[nei] = markValue;
-                    changedCells.append(nei);
-                }
-            }
-        }
-
-
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedCells:"
-        //        << changedCells.size() << endl;
-        //}
-
-        // Loop over changedCells and collect faces
-        DynamicList<label> newChangedFaces(changedCells.size());
-
-        forAll(changedCells, i)
-        {
-            label cellI = changedCells[i];
-
-            const cell& cFaces = mesh().cells()[cellI];
-
-            forAll(cFaces, cFaceI)
-            {
-                label faceI = cFaces[cFaceI];
-
-                if (faceRegion[faceI] == -1)
-                {
-                    faceRegion[faceI] = markValue;
-                    newChangedFaces.append(faceI);
-                }
-            }
-        }
-
-
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedFaces before sync:"
-        //        << changedFaces.size() << endl;
-        //}
-
-
-        // Check for changes to any locally coupled face.
-        // Global connections are done later.
-
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if
-            (
-                isA<cyclicPolyPatch>(pp)
-             && refCast<const cyclicPolyPatch>(pp).owner()
-            )
-            {
-                // Transfer from neighbourPatch to here or vice versa.
-
-                const cyclicPolyPatch& cycPatch =
-                    refCast<const cyclicPolyPatch>(pp);
-
-                label faceI = cycPatch.start();
-
-                forAll(cycPatch, i)
-                {
-                    label otherFaceI = cycPatch.transformGlobalFace(faceI);
-
-                    transferCoupledFaceRegion
-                    (
-                        faceI,
-                        otherFaceI,
-                        faceRegion,
-                        newChangedFaces
-                    );
-
-                    faceI++;
-                }
-            }
-        }
-        forAll(explicitConnections, i)
-        {
-            transferCoupledFaceRegion
-            (
-                explicitConnections[i][0],
-                explicitConnections[i][1],
-                faceRegion,
-                newChangedFaces
-            );
-        }
-
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
-        //        << newChangedFaces.size() << endl;
-        //}
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-Foam::label Foam::regionSplit::calcLocalRegionSplit
 (
+    const globalIndex& globalFaces,
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
 
     labelList& cellRegion
 ) const
 {
-    if (debug)
+    // Field on cells and faces.
-    {
+    List<minData> cellData(mesh().nCells());
-        if (blockedFace.size())
+    List<minData> faceData(mesh().nFaces());
-        {
-            // Check that blockedFace is synced.
-            boolList syncBlockedFace(blockedFace);
-            syncTools::swapFaceList(mesh(), syncBlockedFace);
 
-            forAll(syncBlockedFace, faceI)
+    // Take over blockedFaces by seeding a negative number
+    // (so is always less than the decomposition)
+    label nUnblocked = 0;
+    forAll(faceData, faceI)
     {
-                if (syncBlockedFace[faceI] != blockedFace[faceI])
+        if (blockedFace.size() && blockedFace[faceI])
         {
-                    FatalErrorIn
+            faceData[faceI] = minData(-2);
-                    (
+        }
-                        "regionSplit::calcLocalRegionSplit(..)"
+        else
-                    )   << "Face " << faceI << " not synchronised. My value:"
+        {
-                        << blockedFace[faceI] << "  coupled value:"
+            nUnblocked++;
-                        << syncBlockedFace[faceI]
+        }
-                        << abort(FatalError);
+    }
-                }
+
-            }
+    // Seed unblocked faces
-        }
+    labelList seedFaces(nUnblocked);
-    }
+    List<minData> seedData(nUnblocked);
-
+    nUnblocked = 0;
-    // Region per face.
+
-    // -1 unassigned
+
-    // -2 blocked
+    forAll(faceData, faceI)
-    labelList faceRegion(mesh().nFaces(), -1);
+    {
-
+        if (blockedFace.empty() || !blockedFace[faceI])
-    if (blockedFace.size())
+        {
-    {
+            seedFaces[nUnblocked] = faceI;
-        forAll(blockedFace, faceI)
+            // Seed face with globally unique number
-        {
+            seedData[nUnblocked] = minData(globalFaces.toGlobal(faceI));
-            if (blockedFace[faceI])
+            nUnblocked++;
-            {
+        }
-                faceRegion[faceI] = -2;
+    }
-            }
+
-        }
+
-    }
+    // Propagate information inwards
-
+    FaceCellWave<minData> deltaCalc
-
-    // Assign local regions
-    // ~~~~~~~~~~~~~~~~~~~~
-
-    // Start with region 0
-    label nLocalRegions = 0;
-
-    label unsetCellI = 0;
-
-    do
-    {
-        // Find first unset cell
-
-        for (; unsetCellI < mesh().nCells(); unsetCellI++)
-        {
-            if (cellRegion[unsetCellI] == -1)
-            {
-                break;
-            }
-        }
-
-        if (unsetCellI >= mesh().nCells())
-        {
-            break;
-        }
-
-        fillSeedMask
     (
+        mesh(),
         explicitConnections,
-            cellRegion,
+        false,  // disable walking through cyclicAMI for backwards compatibility
-            faceRegion,
+        seedFaces,
-            unsetCellI,
+        seedData,
-            nLocalRegions
+        faceData,
+        cellData,
+        mesh().globalData().nTotalCells()+1
     );
 
-        // Current unsetCell has now been handled. Go to next region.
-        nLocalRegions++;
-        unsetCellI++;
-    }
-    while (true);
 
-
+    // And extract
-    if (debug)
+    cellRegion.setSize(mesh().nCells());
-    {
     forAll(cellRegion, cellI)
     {
-            if (cellRegion[cellI] < 0)
+        if (cellData[cellI].valid(deltaCalc.data()))
         {
-                FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
+            cellRegion[cellI] = cellData[cellI].data();
-                    << "cell:" << cellI << " region:" << cellRegion[cellI]
-                    << abort(FatalError);
-            }
         }
+        else
+        {
+            // Unvisited cell -> only possible if surrounded by blocked faces.
+            // If so make up region from any of the faces
+            const cell& cFaces = mesh().cells()[cellI];
+            label faceI = cFaces[0];
 
-        forAll(faceRegion, faceI)
+            if (blockedFace.size() && !blockedFace[faceI])
             {
-            if (faceRegion[faceI] == -1)
+                FatalErrorIn("regionSplit::calcNonCompactRegionSplit(..)")
-            {
+                    << "Problem: unblocked face " << faceI
-                FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
+                    << " at " << mesh().faceCentres()[faceI]
-                    << "face:" << faceI << " region:" << faceRegion[faceI]
+                    << " on unassigned cell " << cellI
-                    << abort(FatalError);
+                    << mesh().cellCentres()[faceI]
+                    << exit(FatalError);
+            }
+            cellRegion[cellI] = globalFaces.toGlobal(faceI);
         }
     }
 }
 
-    return nLocalRegions;
-}
-
 
 Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 (
@@ -388,176 +143,142 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 {
     // See header in regionSplit.H
 
-    // 1. Do local analysis
+
-    label nLocalRegions = calcLocalRegionSplit
+    if (!doGlobalRegions)
+    {
+        // Block all parallel faces to avoid comms across
+        boolList coupledOrBlockedFace(blockedFace);
+        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+        if (coupledOrBlockedFace.size())
+        {
+            forAll(pbm, patchI)
+            {
+                const polyPatch& pp = pbm[patchI];
+                if (isA<processorPolyPatch>(pp))
+                {
+                    label faceI = pp.start();
+                    forAll(pp, i)
+                    {
+                        coupledOrBlockedFace[faceI++] = true;
+                    }
+                }
+            }
+        }
+
+        // Create dummy (local only) globalIndex
+        labelList offsets(Pstream::nProcs()+1, 0);
+        for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++)
+        {
+            offsets[i] = mesh().nFaces();
+        }
+        const globalIndex globalRegions(offsets.xfer());
+
+        // Minimise regions across connected cells
+        // Note: still uses global decisions so all processors are running
+        //       in lock-step, i.e. slowest determines overall time.
+        //       To avoid this we could switch off Pstream::parRun.
+        calcNonCompactRegionSplit
         (
+            globalRegions,
+            coupledOrBlockedFace,
+            explicitConnections,
+            cellRegion
+        );
+
+        // Compact
+        Map<label> globalToCompact(mesh().nCells()/8);
+        forAll(cellRegion, cellI)
+        {
+            label region = cellRegion[cellI];
+
+            label globalRegion;
+
+            Map<label>::const_iterator fnd = globalToCompact.find(region);
+            if (fnd == globalToCompact.end())
+            {
+                globalRegion = globalRegions.toGlobal(globalToCompact.size());
+                globalToCompact.insert(region, globalRegion);
+            }
+            else
+            {
+                globalRegion = fnd();
+            }
+            cellRegion[cellI] = globalRegion;
+        }
+
+
+        // Return globalIndex with size = localSize and all regions local
+        labelList compactOffsets(Pstream::nProcs()+1, 0);
+        for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++)
+        {
+            compactOffsets[i] = globalToCompact.size();
+        }
+
+        return autoPtr<globalIndex>(new globalIndex(compactOffsets.xfer()));
+    }
+
+
+
+    // Initial global region numbers
+    const globalIndex globalRegions(mesh().nFaces());
+
+    // Minimise regions across connected cells (including parallel)
+    calcNonCompactRegionSplit
+    (
+        globalRegions,
         blockedFace,
         explicitConnections,
         cellRegion
     );
 
-    if (!doGlobalRegions)
+
+    // Now our cellRegion will have
+    // - non-local regions (i.e. originating from other processors)
+    // - non-compact locally originating regions
+    // so we'll need to compact
+
+    // 4a: count per originating processor the number of regions
+    labelList nOriginating(Pstream::nProcs(), 0);
     {
-        return autoPtr<globalIndex>(new globalIndex(nLocalRegions));
+        labelHashSet haveRegion(mesh().nCells()/8);
-    }
 
-
-    // 2. Assign global regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-    // Offset local regions to create unique global regions.
-
-    globalIndex globalRegions(nLocalRegions);
-
-
-    // Convert regions to global ones
         forAll(cellRegion, cellI)
         {
-        cellRegion[cellI] = globalRegions.toGlobal(cellRegion[cellI]);
+            label region = cellRegion[cellI];
-    }
 
-
+            // Count originating processor. Use isLocal as efficiency since
-    // 3. Merge global regions
+            // most cells are locally originating.
-    // ~~~~~~~~~~~~~~~~~~~~~~~
+            if (globalRegions.isLocal(region))
-    // Regions across non-blocked proc patches get merged.
-    // This will set merged global regions to be the min of both.
-    // (this will create gaps in the global region list so they will get
-    // merged later on)
-
-    while (true)
             {
-        if (debug)
+                if (haveRegion.insert(region))
                 {
-            Pout<< nl << "-- Starting Iteration --" << endl;
+                    nOriginating[Pstream::myProcNo()]++;
-        }
-
-
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
-
-        labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces(), -1);
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (pp.coupled())
-            {
-                const labelList& patchCells = pp.faceCells();
-                SubList<label> patchNbrRegion
-                (
-                    nbrRegion,
-                    pp.size(),
-                    pp.start()-mesh().nInternalFaces()
-                );
-
-                forAll(patchCells, i)
-                {
-                    label faceI = pp.start()+i;
-                    if (!blockedFace.size() || !blockedFace[faceI])
-                    {
-                        patchNbrRegion[i] = cellRegion[patchCells[i]];
                 }
             }
-            }
+            else
-        }
-        syncTools::swapBoundaryFaceList(mesh(), nbrRegion);
-
-        Map<label> globalToMerged(mesh().nFaces()-mesh().nInternalFaces());
-
-        forAll(patches, patchI)
             {
-            const polyPatch& pp = patches[patchI];
+                label procI = globalRegions.whichProcID(region);
-
+                if (haveRegion.insert(region))
-            if (pp.coupled())
                 {
-                const labelList& patchCells = pp.faceCells();
+                    nOriginating[procI]++;
-                SubList<label> patchNbrRegion
-                (
-                    nbrRegion,
-                    pp.size(),
-                    pp.start()-mesh().nInternalFaces()
-                );
-
-                forAll(patchCells, i)
-                {
-                    label faceI = pp.start()+i;
-
-                    if (!blockedFace.size() || !blockedFace[faceI])
-                    {
-                        if (patchNbrRegion[i] < cellRegion[patchCells[i]])
-                        {
-                            //Pout<< "on patch:" << pp.name()
-                            //    << " cell:" << patchCells[i]
-                            //    << " at:"
-                            // << mesh().cellCentres()[patchCells[i]]
-                            //    << " was:" << cellRegion[patchCells[i]]
-                            //    << " nbr:" << patchNbrRegion[i]
-                            //    << endl;
-
-                            globalToMerged.insert
-                            (
-                                cellRegion[patchCells[i]],
-                                patchNbrRegion[i]
-                            );
                 }
             }
         }
     }
-        }
-
-
-        label nMerged = returnReduce(globalToMerged.size(), sumOp<label>());
-
-        if (debug)
-        {
-            Pout<< "nMerged:" << nMerged << endl;
-        }
-
-        if (nMerged == 0)
-        {
-            break;
-        }
-
-        // Renumber the regions according to the globalToMerged
-        forAll(cellRegion, cellI)
-        {
-            label regionI = cellRegion[cellI];
-            Map<label>::const_iterator iter = globalToMerged.find(regionI);
-            if (iter != globalToMerged.end())
-            {
-                 cellRegion[cellI] = iter();
-            }
-        }
-    }
-
-
-    // Now our cellRegion will have non-local elements in it. So compact
-    // it.
-
-    // 4a: count. Use a labelHashSet to count regions only once.
-    label nCompact = 0;
-    {
-        labelHashSet localRegion(mesh().nFaces()-mesh().nInternalFaces());
-        forAll(cellRegion, cellI)
-        {
-            if
-            (
-                globalRegions.isLocal(cellRegion[cellI])
-             && localRegion.insert(cellRegion[cellI])
-            )
-            {
-                nCompact++;
-            }
-        }
-    }
 
     if (debug)
     {
-        Pout<< "Compacted from " << nLocalRegions
+        Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
-            << " down to " << nCompact << " local regions." << endl;
+            << " local regions." << endl;
     }
 
 
     // Global numbering for compacted local regions
-    autoPtr<globalIndex> globalCompactPtr(new globalIndex(nCompact));
+    autoPtr<globalIndex> globalCompactPtr
+    (
+        new globalIndex(nOriginating[Pstream::myProcNo()])
+    );
     const globalIndex& globalCompact = globalCompactPtr();
 
 
@@ -568,12 +289,15 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
     // labelList.
 
     // Local compaction map
-    Map<label> globalToCompact(2*nCompact);
+    Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
     // Remote regions we want the compact number for
     List<labelHashSet> nonLocal(Pstream::nProcs());
     forAll(nonLocal, procI)
     {
-        nonLocal[procI].resize((nLocalRegions-nCompact)/Pstream::nProcs());
+        if (procI != Pstream::myProcNo())
+        {
+            nonLocal[procI].resize(2*nOriginating[procI]);
+        }
     }
 
     forAll(cellRegion, cellI)
@@ -581,15 +305,12 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
         label region = cellRegion[cellI];
         if (globalRegions.isLocal(region))
         {
-            Map<label>::const_iterator iter = globalToCompact.find(region);
+            // Insert new compact region (if not yet present)
-            if (iter == globalToCompact.end())
+            globalToCompact.insert
-            {
-                label compactRegion = globalCompact.toGlobal
             (
-                    globalToCompact.size()
+                region,
+                globalCompact.toGlobal(globalToCompact.size())
             );
-                globalToCompact.insert(region, compactRegion);
-            }
         }
         else
         {
@@ -603,22 +324,17 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
     // Convert the nonLocal (labelHashSets) to labelLists.
 
     labelListList sendNonLocal(Pstream::nProcs());
-    labelList nNonLocal(Pstream::nProcs(), 0);
     forAll(sendNonLocal, procI)
     {
-        sendNonLocal[procI].setSize(nonLocal[procI].size());
+        sendNonLocal[procI] = nonLocal[procI].toc();
-        forAllConstIter(labelHashSet, nonLocal[procI], iter)
-        {
-            sendNonLocal[procI][nNonLocal[procI]++] = iter.key();
-        }
     }
 
     if (debug)
     {
-        forAll(nNonLocal, procI)
+        forAll(sendNonLocal, procI)
         {
             Pout<< "    from processor " << procI
-                << " want " << nNonLocal[procI]
+                << " want " << sendNonLocal[procI].size()
                 << " region numbers."
                 << endl;
         }
@@ -627,7 +343,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 
 
     // Get the wanted region labels into recvNonLocal
-    labelListList recvNonLocal;
+    labelListList recvNonLocal(Pstream::nProcs());
     labelListList sizes;
     Pstream::exchange<labelList, label>
     (
@@ -655,6 +371,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 
     // Send back (into recvNonLocal)
     recvNonLocal.clear();
+    recvNonLocal.setSize(sendWantedLocal.size());
     sizes.clear();
     Pstream::exchange<labelList, label>
     (

src/meshTools/regionSplit/regionSplit.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,9 @@ Description
 
     Can optionally keep all regions local to the processor.
 
+    Note: does not walk across cyclicAMI/cyclicACMI - since not 'coupled()'
+          at the patch level.
+
 
 SourceFiles
     regionSplit.C
@@ -126,30 +129,10 @@ class regionSplit
 
     // Private Member Functions
 
-        //- Transfer faceRegion data from one face to the other (or vice versa)
+        //- Calculate region split in non-compact (global) numbering.
-        void transferCoupledFaceRegion
+        void calcNonCompactRegionSplit
-        (
-            const label faceI,
-            const label otherFaceI,
-
-            labelList& faceRegion,
-            DynamicList<label>& newChangedFaces
-        ) const;
-
-        //- Given a seed cell label, fill cellRegion/faceRegion with markValue
-        //  for contiguous region around it
-        void fillSeedMask
-        (
-            const List<labelPair>& explicitConnections,
-            labelList& cellRegion,
-            labelList& faceRegion,
-            const label seedCellID,
-            const label markValue
-        ) const;
-
-        //- Calculate local region split. Return number of regions.
-        label calcLocalRegionSplit
         (
+            const globalIndex& globalFaces,
             const boolList& blockedFace,
             const List<labelPair>& explicitConnections,
             labelList& cellRegion

src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C

@@ -31,6 +31,7 @@ License
 #include "volFields.H"
 #include "globalIndex.H"
 #include "fvMesh.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -65,11 +66,16 @@ Foam::fileName Foam::externalCoupledFunctionObject::baseDir() const
 Foam::fileName Foam::externalCoupledFunctionObject::groupDir
 (
     const fileName& commsDir,
-    const word& regionName,
+    const word& regionGroupName,
     const wordRe& groupName
 )
 {
-    fileName result(commsDir/regionName/string::validate<fileName>(groupName));
+    fileName result
+    (
+        commsDir
+       /regionGroupName
+       /string::validate<fileName>(groupName)
+    );
     result.clean();
 
     return result;
@@ -126,11 +132,11 @@ void Foam::externalCoupledFunctionObject::removeReadFiles() const
 
     if (log_) Info<< type() << ": removing all read files" << endl;
 
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
+        const word& compName = regionGroupNames_[regionI];
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+
-        const labelList& groups = regionToGroups_[regionName];
+        const labelList& groups = regionToGroups_[compName];
         forAll(groups, i)
         {
             label groupI = groups[i];
@@ -141,7 +147,7 @@ void Foam::externalCoupledFunctionObject::removeReadFiles() const
                 const word& fieldName = groupReadFields_[groupI][fieldI];
                 rm
                 (
-                    groupDir(commsDir_, mesh.dbDir(), groupName)
+                    groupDir(commsDir_, compName, groupName)
                   / fieldName + ".in"
                 );
             }
@@ -159,22 +165,22 @@ void Foam::externalCoupledFunctionObject::removeWriteFiles() const
 
     if (log_) Info<< type() << ": removing all write files" << endl;
 
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
+        const word& compName = regionGroupNames_[regionI];
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+
-        const labelList& groups = regionToGroups_[regionName];
+        const labelList& groups = regionToGroups_[compName];
         forAll(groups, i)
         {
             label groupI = groups[i];
             const wordRe& groupName = groupNames_[groupI];
 
-            forAll(groupWriteFields_[groupI], fieldI)
+            forAll(groupReadFields_[groupI], fieldI)
             {
-                const word& fieldName = groupWriteFields_[groupI][fieldI];
+                const word& fieldName = groupReadFields_[groupI][fieldI];
                 rm
                 (
-                    groupDir(commsDir_, mesh.dbDir(), groupName)
+                    groupDir(commsDir_, compName, groupName)
                   / fieldName + ".out"
                 );
             }
@@ -376,12 +382,21 @@ void Foam::externalCoupledFunctionObject::readLines
 
 void Foam::externalCoupledFunctionObject::writeGeometry
 (
-    const fvMesh& mesh,
+    const UPtrList<const fvMesh>& meshes,
     const fileName& commsDir,
     const wordRe& groupName
 )
 {
-    fileName dir(groupDir(commsDir, mesh.dbDir(), groupName));
+    wordList regionNames(meshes.size());
+    forAll(meshes, i)
+    {
+        regionNames[i] = meshes[i].dbDir();
+    }
+
+    // Make sure meshes are provided in sorted order
+    checkOrder(regionNames);
+
+    fileName dir(groupDir(commsDir, compositeName(regionNames), groupName));
 
     //if (log_)
     {
@@ -397,6 +412,14 @@ void Foam::externalCoupledFunctionObject::writeGeometry
         osFacesPtr.reset(new OFstream(dir/"patchFaces"));
     }
 
+
+    DynamicList<face> allMeshesFaces;
+    DynamicField<point> allMeshesPoints;
+
+    forAll(meshes, meshI)
+    {
+        const fvMesh& mesh = meshes[meshI];
+
         const labelList patchIDs
         (
             mesh.boundaryMesh().patchSet
@@ -405,91 +428,194 @@ void Foam::externalCoupledFunctionObject::writeGeometry
             ).sortedToc()
         );
 
+        // Count faces
+        label nFaces = 0;
         forAll(patchIDs, i)
         {
-        label patchI = patchIDs[i];
+            nFaces += mesh.boundaryMesh()[patchIDs[i]].size();
+        }
 
-        const polyPatch& p = mesh.boundaryMesh()[patchI];
+        // Collect faces
+        DynamicList<label> allFaceIDs(nFaces);
+        forAll(patchIDs, i)
+        {
+            const polyPatch& p = mesh.boundaryMesh()[patchIDs[i]];
+
+            forAll(p, pI)
+            {
+                allFaceIDs.append(p.start()+pI);
+            }
+        }
+
+        // Construct overall patch
+        indirectPrimitivePatch allPatch
+        (
+            IndirectList<face>(mesh.faces(), allFaceIDs),
+            mesh.points()
+        );
 
         labelList pointToGlobal;
         labelList uniquePointIDs;
         mesh.globalData().mergePoints
         (
-            p.meshPoints(),
+            allPatch.meshPoints(),
-            p.meshPointMap(),
+            allPatch.meshPointMap(),
             pointToGlobal,
             uniquePointIDs
         );
 
         label procI = Pstream::myProcNo();
 
-        List<pointField> allPoints(Pstream::nProcs());
+        List<pointField> collectedPoints(Pstream::nProcs());
-        allPoints[procI] = pointField(mesh.points(), uniquePointIDs);
+        collectedPoints[procI] = pointField(mesh.points(), uniquePointIDs);
-        Pstream::gatherList(allPoints);
+        Pstream::gatherList(collectedPoints);
 
-        List<faceList> allFaces(Pstream::nProcs());
+        List<faceList> collectedFaces(Pstream::nProcs());
-        faceList& patchFaces = allFaces[procI];
+        faceList& patchFaces = collectedFaces[procI];
-        patchFaces = p.localFaces();
+        patchFaces = allPatch.localFaces();
         forAll(patchFaces, faceI)
         {
             inplaceRenumber(pointToGlobal, patchFaces[faceI]);
         }
-        Pstream::gatherList(allFaces);
+        Pstream::gatherList(collectedFaces);
 
         if (Pstream::master())
         {
-            pointField pts
+            // Append and renumber
-            (
+            label nPoints = allMeshesPoints.size();
-                ListListOps::combine<pointField>
+
-                (
+            forAll(collectedPoints, procI)
-                    allPoints,
+            {
-                    accessOp<pointField>()
+                allMeshesPoints.append(collectedPoints[procI]);
-                )
+
-            );
+            }
+            face newFace;
+            forAll(collectedFaces, procI)
+            {
+                const faceList& procFaces = collectedFaces[procI];
+
+                forAll(procFaces, faceI)
+                {
+                    const face& f = procFaces[faceI];
+
+                    newFace.setSize(f.size());
+                    forAll(f, fp)
+                    {
+                        newFace[fp] = f[fp]+nPoints;
+                    }
+                    allMeshesFaces.append(newFace);
+                }
+
+                nPoints += collectedPoints[procI].size();
+            }
+        }
 
         //if (log_)
         {
-                Info<< typeName << ": for patch " << p.name()
+            Info<< typeName << ": for mesh " << mesh.name()
-                    << " writing " << pts.size() << " points to "
+                << " writing " << allMeshesPoints.size() << " points to "
                 << osPointsPtr().name() << endl;
+            Info<< typeName << ": for mesh " << mesh.name()
+                << " writing " << allMeshesFaces.size() << " faces to "
+                << osFacesPtr().name() << endl;
+        }
     }
 
     // Write points
-            osPointsPtr() << patchKey.c_str() << p.name() << pts << endl;
+    if (osPointsPtr.valid())
-
-            faceList fcs
-            (
-                ListListOps::combine<faceList>(allFaces, accessOp<faceList>())
-            );
-
-            //if (log_)
     {
-                Info<< typeName << ": for patch " << p.name()
+        osPointsPtr() << allMeshesPoints << endl;
-                    << " writing " << fcs.size() << " faces to "
-                    << osFacesPtr().name() << endl;
     }
 
     // Write faces
-            osFacesPtr() << patchKey.c_str() << p.name() << fcs << endl;
+    if (osFacesPtr.valid())
+    {
+        osFacesPtr() << allMeshesFaces << endl;
     }
 }
+
+
+Foam::word Foam::externalCoupledFunctionObject::compositeName
+(
+    const wordList& regionNames
+)
+{
+    if (regionNames.size() == 0)
+    {
+        FatalErrorIn
+        (
+            "externalCoupledFunctionObject::compositeName(const wordList&)"
+        )   << "Empty regionNames" << abort(FatalError);
+        return word::null;
+    }
+    else if (regionNames.size() == 1)
+    {
+        if (regionNames[0] == polyMesh::defaultRegion)
+        {
+            // For compatibility with single region cases suppress single
+            // region name
+            return word("");
+        }
+        else
+        {
+            return regionNames[0];
+        }
+    }
+    else
+    {
+        // Enforce lexical ordering
+        checkOrder(regionNames);
+
+        word composite(regionNames[0]);
+        for (label i = 1; i < regionNames.size(); i++)
+        {
+            composite += "_" + regionNames[i];
+        }
+
+        return composite;
+    }
+}
+
+
+void Foam::externalCoupledFunctionObject::checkOrder
+(
+    const wordList& regionNames
+)
+{
+    labelList order;
+    sortedOrder(regionNames, order);
+    if (order != identity(regionNames.size()))
+    {
+        FatalErrorIn
+        (
+            "externalCoupledFunctionObject::checkOrder(const wordList&)"
+        )   << "regionNames " << regionNames << " not in alphabetical order :"
+            << order << exit(FatalError);
+    }
 }
 
 
 void Foam::externalCoupledFunctionObject::readData()
 {
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
+        const word& compName = regionGroupNames_[regionI];
-        const labelList& groups = regionToGroups_[regionName];
+        const wordList& regionNames = regionGroupRegions_[regionI];
 
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+        // Get the meshes for the region-group
+        UPtrList<const fvMesh> meshes(regionNames.size());
+        forAll(regionNames, j)
+        {
+            const word& regionName = regionNames[j];
+            meshes.set(j, &time_.lookupObject<fvMesh>(regionName));
+        }
+
+        const labelList& groups = regionToGroups_[compName];
 
         forAll(groups, i)
         {
             label groupI = groups[i];
             const wordRe& groupName = groupNames_[groupI];
-            const labelList& patchIDs = groupPatchIDs_[groupI];
             const wordList& fieldNames = groupReadFields_[groupI];
 
             forAll(fieldNames, fieldI)
@@ -498,37 +624,32 @@ void Foam::externalCoupledFunctionObject::readData()
 
                 bool ok = readData<scalar>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<vector>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<sphericalTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<symmTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<tensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
 
@@ -538,7 +659,7 @@ void Foam::externalCoupledFunctionObject::readData()
                     (
                         "void Foam::externalCoupledFunctionObject::readData()"
                     )
-                        << "Field " << fieldName << " in region " << mesh.name()
+                        << "Field " << fieldName << " in regions " << compName
                         << " was not found." << endl;
                 }
             }
@@ -549,56 +670,59 @@ void Foam::externalCoupledFunctionObject::readData()
 
 void Foam::externalCoupledFunctionObject::writeData() const
 {
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
+        const word& compName = regionGroupNames_[regionI];
-        const labelList& groups = regionToGroups_[regionName];
+        const wordList& regionNames = regionGroupRegions_[regionI];
 
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+        // Get the meshes for the region-group
+        UPtrList<const fvMesh> meshes(regionNames.size());
+        forAll(regionNames, j)
+        {
+            const word& regionName = regionNames[j];
+            meshes.set(j, &time_.lookupObject<fvMesh>(regionName));
+        }
+
+        const labelList& groups = regionToGroups_[compName];
 
         forAll(groups, i)
         {
             label groupI = groups[i];
             const wordRe& groupName = groupNames_[groupI];
-            const labelList& patchIDs = groupPatchIDs_[groupI];
             const wordList& fieldNames = groupWriteFields_[groupI];
 
             forAll(fieldNames, fieldI)
             {
                 const word& fieldName = fieldNames[fieldI];
+
                 bool ok = writeData<scalar>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<vector>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<sphericalTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<symmTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<tensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
 
@@ -608,7 +732,7 @@ void Foam::externalCoupledFunctionObject::writeData() const
                     (
                         "void Foam::externalCoupledFunctionObject::writeData()"
                     )
-                        << "Field " << fieldName << " in region " << mesh.name()
+                        << "Field " << fieldName << " in regions " << compName
                         << " was not found." << endl;
                 }
             }
@@ -625,12 +749,20 @@ void Foam::externalCoupledFunctionObject::initialise()
     }
 
     // Write the geometry if not already there
-    forAll(regionNames_, regionI)
+    forAll(regionGroupRegions_, i)
     {
-        const word& regionName = regionNames_[regionI];
+        const word& compName = regionGroupNames_[i];
-        const labelList& groups = regionToGroups_[regionName];
+        const wordList& regionNames = regionGroupRegions_[i];
 
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+        // Get the meshes for the region-group
+        UPtrList<const fvMesh> meshes(regionNames.size());
+        forAll(regionNames, j)
+        {
+            const word& regionName = regionNames[j];
+            meshes.set(j, &time_.lookupObject<fvMesh>(regionName));
+        }
+
+        const labelList& groups = regionToGroups_[compName];
 
         forAll(groups, i)
         {
@@ -640,14 +772,16 @@ void Foam::externalCoupledFunctionObject::initialise()
             bool exists = false;
             if (Pstream::master())
             {
-                fileName dir(groupDir(commsDir_, mesh.dbDir(), groupName));
+                fileName dir(groupDir(commsDir_, compName, groupName));
 
-                exists = isFile(dir/"patchPoints") || isFile(dir/"patchFaces");
+                exists =
+                    isFile(dir/"patchPoints")
+                 || isFile(dir/"patchFaces");
             }
 
             if (!returnReduce(exists, orOp<bool>()))
             {
-                writeGeometry(mesh, commsDir_, groupName);
+                writeGeometry(meshes, commsDir_, groupName);
             }
         }
     }
@@ -802,6 +936,12 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
     initByExternal_ = readBool(dict.lookup("initByExternal"));
     log_ = dict.lookupOrDefault("log", false);
 
+
+    // Get names of all fvMeshes (and derived types)
+    wordList allRegionNames(time_.lookupClass<fvMesh>().sortedToc());
+
+
+
     const dictionary& allRegionsDict = dict.subDict("regions");
 
     forAllConstIter(dictionary, allRegionsDict, iter)
@@ -818,9 +958,16 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
                 << exit(FatalIOError);
         }
 
-        const word& regionName = iter().keyword();
+        const wordRe regionGroupName(iter().keyword());
         const dictionary& regionDict = iter().dict();
-        regionNames_.append(regionName);
+
+        labelList regionIDs = findStrings(regionGroupName, allRegionNames);
+
+        const wordList regionNames(allRegionNames, regionIDs);
+
+        regionGroupNames_.append(compositeName(regionNames));
+        regionGroupRegions_.append(regionNames);
+
 
         forAllConstIter(dictionary, regionDict, regionIter)
         {
@@ -844,7 +991,7 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
 
             HashTable<labelList>::iterator fnd = regionToGroups_.find
             (
-                regionName
+                regionGroupNames_.last()
             );
             if (fnd != regionToGroups_.end())
             {
@@ -852,21 +999,15 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
             }
             else
             {
-                regionToGroups_.insert(regionName, labelList(1, nGroups));
+                regionToGroups_.insert
+                (
+                    regionGroupNames_.last(),
+                    labelList(1, nGroups)
+                );
             }
             groupNames_.append(groupName);
             groupReadFields_.append(readFields);
             groupWriteFields_.append(writeFields);
-
-            // Pre-calculate the patchIDs
-            const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
-            groupPatchIDs_.append
-            (
-                mesh.boundaryMesh().patchSet
-                (
-                    List<wordRe>(1, groupName)
-                ).sortedToc()
-            );
         }
     }
 
@@ -875,25 +1016,26 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
     if (log_)
     {
         Info<< type() << ": Communicating with regions:" << endl;
-        forAll(regionNames_, regionI)
+        forAll(regionGroupNames_, rgI)
         {
-            const word& regionName = regionNames_[regionI];
+            //const wordList& regionNames = regionGroupRegions_[rgI];
-            const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+            const word& compName = regionGroupNames_[rgI];
 
-            Info<< "Region: " << mesh.name() << endl << incrIndent;
+            Info<< "Region: " << compName << endl << incrIndent;
-            const labelList& groups = regionToGroups_[regionName];
+            const labelList& groups = regionToGroups_[compName];
             forAll(groups, i)
             {
                 label groupI = groups[i];
                 const wordRe& groupName = groupNames_[groupI];
-                const labelList& patchIDs = groupPatchIDs_[groupI];
 
-                Info<< indent << "Group: " << groupName << "\t"
+                Info<< indent << "patchGroup: " << groupName << "\t"
-                    << " patches: " << patchIDs << endl
-                    << incrIndent
-                    << indent << "Reading fields: " << groupReadFields_[groupI]
                     << endl
-                    << indent << "Writing fields: " << groupWriteFields_[groupI]
+                    << incrIndent
+                    << indent << "Reading fields: "
+                    << groupReadFields_[groupI]
+                    << endl
+                    << indent << "Writing fields: "
+                    << groupWriteFields_[groupI]
                     << endl
                     << decrIndent;
             }
@@ -907,17 +1049,17 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
     //       should already be written - but just make sure
     if (Pstream::master())
     {
-        forAll(regionNames_, regionI)
+        forAll(regionGroupNames_, rgI)
         {
-            const word& regionName = regionNames_[regionI];
+            const word& compName = regionGroupNames_[rgI];
-            const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+
-            const labelList& groups = regionToGroups_[regionName];
+            const labelList& groups = regionToGroups_[compName];
             forAll(groups, i)
             {
                 label groupI = groups[i];
                 const wordRe& groupName = groupNames_[groupI];
 
-                fileName dir(groupDir(commsDir_, mesh.dbDir(), groupName));
+                fileName dir(groupDir(commsDir_, compName, groupName));
                 if (!isDir(dir))
                 {
                     if (log_)
@@ -925,7 +1067,6 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
                         Info<< type() << ": creating communications directory "
                             << dir << endl;
                     }
-
                     mkDir(dir);
                 }
             }

src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.H

@@ -48,7 +48,10 @@ Description
     which gets read/written on the master processor only. In the
     communications directory the structure will be
 
-        <regionName>/<patchGroup>/<fieldName>.[in|out]
+        <regionsName>/<patchGroup>/<fieldName>.[in|out]
+
+    (where regionsName is either the name of a single region or a composite
+     of multiple region names)
 
     At start-up, the boundary creates a lock file, i.e..
 
@@ -58,13 +61,13 @@ Description
     execution the boundary values are written to files (one per region,
     per patch(group), per field), e.g.
 
-        <regionName>/<patchGroup>/<fieldName>.out
+        <regionsName>/<patchGroup>/<fieldName>.out
 
     The lock file is then removed, instructing the external source to take
     control of the program execution. When ready, the external program
     should create the return values, e.g. to files
 
-        <regionName>/<patchGroup>/<fieldName>.in
+        <regionsName>/<patchGroup>/<fieldName>.in
 
     ... and then re-instate the lock file. The functionObject will then
     read these values, apply them to the boundary conditions and pass
@@ -82,7 +85,7 @@ Description
 
         regions
         {
-            region0
+            "(region1|region0)"         // Name of region(s)
             {
                 TPatchGroup             // Name of patch(group)
                 {
@@ -95,7 +98,7 @@ Description
     \endverbatim
 
     This reads/writes (on the master processor) the directory:
-        comms/region0/TPatchGroup/
+        comms/region0_region1/TPatchGroup/
     with contents:
         patchPoints     (collected points)
         patchFaces      (collected faces)
@@ -120,6 +123,7 @@ SourceFiles
 #include "wordReList.H"
 #include "scalarField.H"
 #include "Switch.H"
+#include "UPtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -167,18 +171,18 @@ class externalCoupledFunctionObject
         //- Log flag
         bool log_;
 
-        //- Names of regions
+        //- Names of (composite) regions
-        DynamicList<word> regionNames_;
+        DynamicList<word> regionGroupNames_;
 
-        // Per region the indices of the group information
+        // Per (composite) region the names of the regions
+        DynamicList<wordList> regionGroupRegions_;
+
+        // Per (composite) region the indices of the group information
         HashTable<labelList> regionToGroups_;
 
         // Per group the names of the patches/patchGroups
         DynamicList<wordRe> groupNames_;
 
-        // Per group the indices of the patches
-        DynamicList<labelList> groupPatchIDs_;
-
         // Per group the names of the fields to read
         DynamicList<wordList> groupReadFields_;
 
@@ -195,7 +199,7 @@ class externalCoupledFunctionObject
         static fileName groupDir
         (
             const fileName& commsDir,
-            const word& regionName,
+            const word& regionsName,
             const wordRe& groupName
         );
 
@@ -225,9 +229,8 @@ class externalCoupledFunctionObject
         template<class Type>
         bool readData
         (
-            const fvMesh& mesh,
+            const UPtrList<const fvMesh>& meshes,
             const wordRe& groupName,
-            const labelList& patchIDs,
             const word& fieldName
         );
         //- Read data for all regions, all fields
@@ -237,9 +240,8 @@ class externalCoupledFunctionObject
         template<class Type>
         bool writeData
         (
-            const fvMesh& mesh,
+            const UPtrList<const fvMesh>& meshes,
             const wordRe& groupName,
-            const labelList& patchIDs,
             const word& fieldName
         ) const;
 
@@ -275,6 +277,7 @@ class externalCoupledFunctionObject
         template<class Type>
         static tmp<Field<Type> > gatherAndCombine(const Field<Type>& fld);
 
+        static void checkOrder(const wordList&);
 
         //- Disallow default bitwise copy construc
         externalCoupledFunctionObject(const externalCoupledFunctionObject&);
@@ -356,10 +359,14 @@ public:
 
         // Other
 
+            //- Create single name by appending words (in sorted order),
+            //  separated by '_'
+            static word compositeName(const wordList&);
+
             //- Write geometry for the group/patch
             static void writeGeometry
             (
-                const fvMesh& mesh,
+                const UPtrList<const fvMesh>& meshes,
                 const fileName& commsDir,
                 const wordRe& groupName
             );

src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObjectTemplates.C

@@ -40,25 +40,20 @@ License
 template<class Type>
 bool Foam::externalCoupledFunctionObject::readData
 (
-    const fvMesh& mesh,
+    const UPtrList<const fvMesh>& meshes,
     const wordRe& groupName,
-    const labelList& patchIDs,
     const word& fieldName
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
     typedef externalCoupledMixedFvPatchField<Type> patchFieldType;
 
-    if (!mesh.foundObject<volFieldType>(fieldName))
+    wordList regionNames(meshes.size());
+    forAll(meshes, i)
     {
-        return false;
+        regionNames[i] = meshes[i].dbDir();
     }
 
-    const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
-    const typename volFieldType::GeometricBoundaryField& bf =
-        cvf.boundaryField();
-
-
     // File only opened on master; contains data for all processors, for all
     // patchIDs.
     autoPtr<IFstream> masterFilePtr;
@@ -66,7 +61,7 @@ bool Foam::externalCoupledFunctionObject::readData
     {
         const fileName transferFile
         (
-            groupDir(commsDir_, mesh.dbDir(), groupName)
+            groupDir(commsDir_, compositeName(regionNames), groupName)
           / fieldName + ".in"
         );
 
@@ -82,18 +77,46 @@ bool Foam::externalCoupledFunctionObject::readData
             (
                 "void externalCoupledFunctionObject::readData"
                 "("
-                    "const fvMesh&, "
+                    "const UPtrList<const fvMesh>&, "
                     "const wordRe&, "
-                    "const labelList&, "
                     "const word&"
                ")",
                 masterFilePtr()
-            )   << "Cannot open file for region " << mesh.name()
+            )   << "Cannot open file for region " << compositeName(regionNames)
-                << ", field " << fieldName << ", patches " << patchIDs
+                << ", field " << fieldName
                 << exit(FatalIOError);
         }
     }
 
+
+    label nFound = 0;
+
+
+    forAll(meshes, i)
+    {
+        const fvMesh& mesh = meshes[i];
+
+        if (!mesh.foundObject<volFieldType>(fieldName))
+        {
+            continue;
+        }
+
+        nFound++;
+
+        const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
+        const typename volFieldType::GeometricBoundaryField& bf =
+            cvf.boundaryField();
+
+
+        // Get the patches
+        const labelList patchIDs
+        (
+            mesh.boundaryMesh().patchSet
+            (
+                List<wordRe>(1, groupName)
+            ).sortedToc()
+        );
+
         // Handle column-wise reading of patch data. Supports most easy types
         forAll(patchIDs, i)
         {
@@ -101,8 +124,8 @@ bool Foam::externalCoupledFunctionObject::readData
 
             if (isA<patchFieldType>(bf[patchI]))
             {
-            // Explicit handling of externalCoupledObjectMixed bcs - they have
+                // Explicit handling of externalCoupledMixed bcs - they
-            // specialised reading routines.
+                // have specialised reading routines.
 
                 patchFieldType& pf = const_cast<patchFieldType&>
                 (
@@ -141,7 +164,8 @@ bool Foam::externalCoupledFunctionObject::readData
                     data
                 );
 
-            mixedFvPatchField<Type>& pf = const_cast<mixedFvPatchField<Type>&>
+                mixedFvPatchField<Type>& pf =
+                const_cast<mixedFvPatchField<Type>&>
                 (
                     refCast<const mixedFvPatchField<Type> >
                     (
@@ -154,12 +178,22 @@ bool Foam::externalCoupledFunctionObject::readData
                 direction columnI = 2*pTraits<Type>::nComponents;
 
                 Field<Type>& refValue = pf.refValue();
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
                 {
                     refValue.replace(cmpt, data[columnI++]);
                 }
                 Field<Type>& refGrad = pf.refGrad();
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
                 {
                     refGrad.replace(cmpt, data[columnI++]);
                 }
@@ -192,9 +226,18 @@ bool Foam::externalCoupledFunctionObject::readData
 
                 // Transfer gradient to bc
                 Field<Type>& gradient = pf.gradient();
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
                 {
-                gradient.replace(cmpt, data[pTraits<Type>::nComponents+cmpt]);
+                    gradient.replace
+                    (
+                        cmpt,
+                        data[pTraits<Type>::nComponents+cmpt]
+                    );
                 }
 
                 // Update the value from the read coefficicient. Bypass any
@@ -215,7 +258,12 @@ bool Foam::externalCoupledFunctionObject::readData
 
                 // Transfer read value to bc
                 Field<Type> value(bf[patchI].size());
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
                 {
                     value.replace(cmpt, data[cmpt]);
                 }
@@ -241,9 +289,8 @@ bool Foam::externalCoupledFunctionObject::readData
                 (
                     "void externalCoupledFunctionObject::readData"
                     "("
-                    "const fvMesh&, "
+                        "const UPtrList<const fvMesh>&, "
                         "const wordRe&, "
-                    "const labelList&, "
                         "const word&"
                    ")"
                 )
@@ -255,8 +302,9 @@ bool Foam::externalCoupledFunctionObject::readData
 
             initialised_ = true;
         }
+    }
 
-    return true;
+    return nFound > 0;
 }
 
 
@@ -308,25 +356,20 @@ Foam::externalCoupledFunctionObject::gatherAndCombine
 template<class Type>
 bool Foam::externalCoupledFunctionObject::writeData
 (
-    const fvMesh& mesh,
+    const UPtrList<const fvMesh>& meshes,
     const wordRe& groupName,
-    const labelList& patchIDs,
     const word& fieldName
 ) const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
     typedef externalCoupledMixedFvPatchField<Type> patchFieldType;
 
-    if (!mesh.foundObject<volFieldType>(fieldName))
+    wordList regionNames(meshes.size());
+    forAll(meshes, i)
     {
-        return false;
+        regionNames[i] = meshes[i].dbDir();
     }
 
-    const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
-    const typename volFieldType::GeometricBoundaryField& bf =
-        cvf.boundaryField();
-
-
     // File only opened on master; contains data for all processors, for all
     // patchIDs
     autoPtr<OFstream> masterFilePtr;
@@ -334,7 +377,7 @@ bool Foam::externalCoupledFunctionObject::writeData
     {
         const fileName transferFile
         (
-            groupDir(commsDir_, mesh.dbDir(), groupName)
+            groupDir(commsDir_, compositeName(regionNames), groupName)
           / fieldName + ".out"
         );
 
@@ -350,14 +393,13 @@ bool Foam::externalCoupledFunctionObject::writeData
             (
                 "externalCoupledFunctionObject::writeData"
                 "("
-                    "const fvMesh&, "
+                    "const UPtrList<const fvMesh>&, "
                     "const wordRe&, "
-                    "const labelList&, "
                     "const word&"
                 ") const",
                 masterFilePtr()
-            )   << "Cannot open file for region " << mesh.name()
+            )   << "Cannot open file for region " << compositeName(regionNames)
-                << ", field " << fieldName << ", patches " << patchIDs
+                << ", field " << fieldName
                 << exit(FatalIOError);
         }
     }
@@ -365,6 +407,33 @@ bool Foam::externalCoupledFunctionObject::writeData
 
     bool headerDone = false;
 
+    label nFound = 0;
+
+    forAll(meshes, i)
+    {
+        const fvMesh& mesh = meshes[i];
+
+        if (!mesh.foundObject<volFieldType>(fieldName))
+        {
+            continue;
+        }
+
+        nFound++;
+
+        const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
+        const typename volFieldType::GeometricBoundaryField& bf =
+            cvf.boundaryField();
+
+
+        // Get the patches
+        const labelList patchIDs
+        (
+            mesh.boundaryMesh().patchSet
+            (
+                List<wordRe>(1, groupName)
+            ).sortedToc()
+        );
+
         // Handle column-wise writing of patch data. Supports most easy types
         forAll(patchIDs, i)
         {
@@ -374,8 +443,8 @@ bool Foam::externalCoupledFunctionObject::writeData
 
             if (isA<patchFieldType>(bf[patchI]))
             {
-            // Explicit handling of externalCoupledObjectMixed bcs - they have
+                // Explicit handling of externalCoupledMixed bcs - they
-            // specialised writing routines
+                // have specialised writing routines
 
                 const patchFieldType& pf = refCast<const patchFieldType>
                 (
@@ -451,8 +520,9 @@ bool Foam::externalCoupledFunctionObject::writeData
                 }
             }
         }
+    }
 
-    return true;
+    return nFound > 0;
 }
 
 

src/sampling/sampledSet/face/faceOnlySet.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -169,6 +169,7 @@ void Foam::faceOnlySet::calcSamples
         // Set points and cell/face labels to empty lists
         //Info<< "calcSamples : Both start_ and end_ outside domain"
         //    << endl;
+        const_cast<polyMesh&>(mesh()).moving(oldMoving);
 
         return;
     }

src/sampling/sampledSet/patchSeed/patchSeedSet.C

@@ -30,9 +30,8 @@ License
 #include "treeDataFace.H"
 #include "Time.H"
 #include "meshTools.H"
-//#include "Random.H"
-// For 'facePoint' helper function only
 #include "mappedPatchBase.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -85,19 +84,145 @@ void Foam::patchSeedSet::calcSamples
     }
 
 
-    label totalSize = returnReduce(sz, sumOp<label>());
+    if (!rndGenPtr_.valid())
+    {
+        rndGenPtr_.reset(new Random(0));
+    }
+    Random& rndGen = rndGenPtr_();
+
+
+    if (selectedLocations_.size())
+    {
+        DynamicList<label> newPatchFaces(patchFaces.size());
+
+        // Find the nearest patch face
+        {
+            // 1. All processors find nearest local patch face for all
+            //    selectedLocations
+
+            // All the info for nearest. Construct to miss
+            List<mappedPatchBase::nearInfo> nearest(selectedLocations_.size());
+
+            const indirectPrimitivePatch pp
+            (
+                IndirectList<face>(mesh().faces(), patchFaces),
+                mesh().points()
+            );
+
+            treeBoundBox patchBb
+            (
+                treeBoundBox(pp.points(), pp.meshPoints()).extend
+                (
+                    rndGen,
+                    1e-4
+                )
+            );
+            patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+            indexedOctree<treeDataFace> boundaryTree
+            (
+                treeDataFace    // all information needed to search faces
+                (
+                    false,      // do not cache bb
+                    mesh(),
+                    patchFaces  // boundary faces only
+                ),
+                patchBb,        // overall search domain
+                8,              // maxLevel
+                10,             // leafsize
+                3.0             // duplicity
+            );
+
+            // Get some global dimension so all points are equally likely
+            // to be found
+            const scalar globalDistSqr
+            (
+                //magSqr
+                //(
+                //    boundBox
+                //    (
+                //        pp.points(),
+                //        pp.meshPoints(),
+                //        true
+                //    ).span()
+                //)
+                GREAT
+            );
+
+            forAll(selectedLocations_, sampleI)
+            {
+                const point& sample = selectedLocations_[sampleI];
+
+                pointIndexHit& nearInfo = nearest[sampleI].first();
+                nearInfo = boundaryTree.findNearest
+                (
+                    sample,
+                    globalDistSqr
+                );
+
+                if (!nearInfo.hit())
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() =
+                        Pstream::myProcNo();
+                }
+                else
+                {
+                    point fc(pp[nearInfo.index()].centre(pp.points()));
+                    nearInfo.setPoint(fc);
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() =
+                        Pstream::myProcNo();
+                }
+            }
+
+
+            // 2. Reduce on master. Select nearest processor.
+
+            // Find nearest. Combine on master.
+            Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
+            Pstream::listCombineScatter(nearest);
+
+
+            // 3. Pick up my local faces that have won
+
+            forAll(nearest, sampleI)
+            {
+                if (nearest[sampleI].first().hit())
+                {
+                    label procI = nearest[sampleI].second().second();
+                    label index = nearest[sampleI].first().index();
+
+                    if (procI == Pstream::myProcNo())
+                    {
+                        newPatchFaces.append(pp.addressing()[index]);
+                    }
+                }
+            }
+        }
+
+        if (debug)
+        {
+            Pout<< "Found " << newPatchFaces.size()
+                << " out of " << selectedLocations_.size()
+                << " on local processor" << endl;
+        }
+
+        patchFaces.transfer(newPatchFaces);
+    }
 
 
     // Shuffle and truncate if in random mode
+    label totalSize = returnReduce(patchFaces.size(), sumOp<label>());
+
     if (maxPoints_ < totalSize)
     {
         // Check what fraction of maxPoints_ I need to generate locally.
-        label myMaxPoints = label(scalar(sz)/totalSize*maxPoints_);
+        label myMaxPoints =
+            label(scalar(patchFaces.size())/totalSize*maxPoints_);
 
-        rndGenPtr_.reset(new Random(123456));
+        labelList subset = identity(patchFaces.size());
-        Random& rndGen = rndGenPtr_();
-
-        labelList subset = identity(sz);
         for (label iter = 0; iter < 4; iter++)
         {
             forAll(subset, i)
@@ -115,7 +240,7 @@ void Foam::patchSeedSet::calcSamples
         if (debug)
         {
             Pout<< "In random mode : selected " << patchFaces.size()
-                << " faces out of " << sz << endl;
+                << " faces out of " << patchFaces.size() << endl;
         }
     }
 
@@ -135,6 +260,9 @@ void Foam::patchSeedSet::calcSamples
     forAll(patchFaces, i)
     {
         label faceI = patchFaces[i];
+
+        // Slightly shift point in since on warped face face-diagonal
+        // decomposition might be outside cell for face-centre decomposition!
         pointIndexHit info = mappedPatchBase::facePoint
         (
             mesh(),
@@ -217,9 +345,15 @@ Foam::patchSeedSet::patchSeedSet
             wordReList(dict.lookup("patches"))
         )
     ),
-    //searchDist_(readScalar(dict.lookup("maxDistance"))),
+    maxPoints_(readLabel(dict.lookup("maxPoints"))),
-    //offsetDist_(readScalar(dict.lookup("offsetDist"))),
+    selectedLocations_
-    maxPoints_(readLabel(dict.lookup("maxPoints")))
+    (
+        dict.lookupOrDefault<pointField>
+        (
+            "points",
+            pointField(0)
+        )
+    )
 {
     genSamples();
 

src/sampling/sampledSet/patchSeed/patchSeedSet.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,18 +58,15 @@ class patchSeedSet
         //- Patches to sample
         const labelHashSet patchSet_;
 
-//        //- Maximum distance to look for nearest
+        //- Maximum number of patch faces to seed (if in random subset mode)
-//        const scalar searchDist_;
-//
-//        //- Offset distance
-//        const scalar offsetDist_;
-
-        //- Maximum number of patch faces to seed
         const label maxPoints_;
 
         //- Random number generator (if maxPoints < num patch faces)
         autoPtr<Random> rndGenPtr_;
 
+        //- Patch faces to seed selected based on nearness to supplied points
+        const pointField selectedLocations_;
+
 
     // Private Member Functions
 

src/sampling/sampledSet/uniform/uniformSet.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -276,6 +276,8 @@ void Foam::uniformSet::calcSamples
         // (or is along edge)
         // Set points and cell/face labels to empty lists
 
+        const_cast<polyMesh&>(mesh()).moving(oldMoving);
+
         return;
     }
 

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -278,7 +278,8 @@ void Foam::distanceSurface::createGeometry()
                 cellDistance,
                 pointDistance_,
                 distance_,
-                regularise_
+                regularise_,
+                bounds_
             )
         );
     }
@@ -291,7 +292,8 @@ void Foam::distanceSurface::createGeometry()
                 cellDistance,
                 pointDistance_,
                 distance_,
-                regularise_
+                regularise_,
+                bounds_
             )
         );
     }
@@ -336,6 +338,7 @@ Foam::distanceSurface::distanceSurface
     cell_(dict.lookupOrDefault("cell", true)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     zoneKey_(keyType::null),
     needsUpdate_(true),
     isoSurfCellPtr_(NULL),
@@ -364,7 +367,8 @@ Foam::distanceSurface::distanceSurface
     const bool signedDistance,
     const bool cell,
     const Switch regularise,
-    const Switch average
+    const Switch average,
+    const boundBox& bounds
 )
 :
     sampledSurface(name, mesh, interpolate),
@@ -390,6 +394,7 @@ Foam::distanceSurface::distanceSurface
     cell_(cell),
     regularise_(regularise),
     average_(average),
+    bounds_(bounds),
     zoneKey_(keyType::null),
     needsUpdate_(true),
     isoSurfCellPtr_(NULL),

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -75,6 +75,9 @@ class distanceSurface
         //- Whether to recalculate cell values as average of point values
         const Switch average_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- If restricted to zones, name of this zone or a regular expression
         keyType zoneKey_;
 
@@ -144,7 +147,8 @@ public:
             const bool signedDistance,
             const bool cell,
             const Switch regularise,
-            const Switch average
+            const Switch average,
+            const boundBox& bounds = boundBox::greatBox
         );
 
 

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,12 +32,36 @@ License
 #include "surfaceFields.H"
 #include "OFstream.H"
 #include "meshTools.H"
+#include "triSurfaceSearch.H"
+#include "surfaceIntersection.H"
+#include "intersectedSurface.H"
+#include "searchableBox.H"
+#include "triSurfaceMesh.H"
+#include "triPoints.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
 defineTypeNameAndDebug(isoSurface, 0);
+
+// Helper class for slicing triangles
+class storeOp
+{
+public:
+    DynamicList<triPoints>& tris_;
+
+    inline storeOp(DynamicList<triPoints>& tris)
+    :
+        tris_(tris)
+    {}
+
+    inline void operator()(const triPoints& tri)
+    {
+        tris_.append(tri);
+    }
+};
+
 }
 
 
@@ -473,43 +497,6 @@ void Foam::isoSurface::calcCutTypes
 }
 
 
-// Return the two common points between two triangles
-Foam::labelPair Foam::isoSurface::findCommonPoints
-(
-    const labelledTri& tri0,
-    const labelledTri& tri1
-)
-{
-    labelPair common(-1, -1);
-
-    label fp0 = 0;
-    label fp1 = findIndex(tri1, tri0[fp0]);
-
-    if (fp1 == -1)
-    {
-        fp0 = 1;
-        fp1 = findIndex(tri1, tri0[fp0]);
-    }
-
-    if (fp1 != -1)
-    {
-        // So tri0[fp0] is tri1[fp1]
-
-        // Find next common point
-        label fp0p1 = tri0.fcIndex(fp0);
-        label fp1p1 = tri1.fcIndex(fp1);
-        label fp1m1 = tri1.rcIndex(fp1);
-
-        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
-        {
-            common[0] = tri0[fp0];
-            common[1] = tri0[fp0p1];
-        }
-    }
-    return common;
-}
-
-
 // Caculate centre of surface.
 Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
 {
@@ -523,73 +510,6 @@ Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
 }
 
 
-// Replace surface (localPoints, localTris) with single point. Returns
-// point. Destructs arguments.
-Foam::pointIndexHit Foam::isoSurface::collapseSurface
-(
-    pointField& localPoints,
-    DynamicList<labelledTri, 64>& localTris
-)
-{
-    pointIndexHit info(false, vector::zero, localTris.size());
-
-    if (localTris.size() == 1)
-    {
-        const labelledTri& tri = localTris[0];
-        info.setPoint(tri.centre(localPoints));
-        info.setHit();
-    }
-    else if (localTris.size() == 2)
-    {
-        // Check if the two triangles share an edge.
-        const labelledTri& tri0 = localTris[0];
-        const labelledTri& tri1 = localTris[0];
-
-        labelPair shared = findCommonPoints(tri0, tri1);
-
-        if (shared[0] != -1)
-        {
-            info.setPoint
-            (
-                0.5
-              * (
-                    tri0.centre(localPoints)
-                  + tri1.centre(localPoints)
-                )
-            );
-            info.setHit();
-        }
-    }
-    else if (localTris.size())
-    {
-        // Check if single region. Rare situation.
-        triSurface surf
-        (
-            localTris,
-            geometricSurfacePatchList(0),
-            localPoints,
-            true
-        );
-        localTris.clearStorage();
-
-        labelList faceZone;
-        label nZones = surf.markZones
-        (
-            boolList(surf.nEdges(), false),
-            faceZone
-        );
-
-        if (nZones == 1)
-        {
-            info.setPoint(calcCentre(surf));
-            info.setHit();
-        }
-    }
-
-    return info;
-}
-
-
 // Determine per cell centre whether all the intersections get collapsed
 // to a single point
 void Foam::isoSurface::calcSnappedCc
@@ -1133,6 +1053,220 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
 }
 
 
+void Foam::isoSurface::trimToPlanes
+(
+    const PtrList<plane>& planes,
+    const triPointRef& tri,
+    DynamicList<point>& newTriPoints
+)
+{
+    // Buffer for generated triangles
+    DynamicList<triPoints> insideTrisA;
+    storeOp insideOpA(insideTrisA);
+
+    // Buffer for generated triangles
+    DynamicList<triPoints> insideTrisB;
+    storeOp insideOpB(insideTrisB);
+
+    triPointRef::dummyOp dop;
+
+    // Store starting triangle in insideTrisA
+    insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
+
+
+    bool useA = true;
+
+    forAll(planes, faceI)
+    {
+        const plane& pl = planes[faceI];
+
+        if (useA)
+        {
+            insideOpB.tris_.clear();
+            forAll(insideOpA.tris_, i)
+            {
+                const triPoints& tri = insideOpA.tris_[i];
+                triPointRef(tri).sliceWithPlane(pl, insideOpB, dop);
+            }
+        }
+        else
+        {
+            insideOpA.tris_.clear();
+            forAll(insideOpB.tris_, i)
+            {
+                const triPoints& tri = insideOpB.tris_[i];
+                triPointRef(tri).sliceWithPlane(pl, insideOpA, dop);
+            }
+        }
+        useA = !useA;
+    }
+
+
+    // Transfer
+    if (useA)
+    {
+        forAll(insideOpA.tris_, i)
+        {
+            const triPoints& tri = insideOpA.tris_[i];
+            newTriPoints.append(tri[0]);
+            newTriPoints.append(tri[1]);
+            newTriPoints.append(tri[2]);
+        }
+    }
+    else
+    {
+        forAll(insideOpB.tris_, i)
+        {
+            const triPoints& tri = insideOpB.tris_[i];
+            newTriPoints.append(tri[0]);
+            newTriPoints.append(tri[1]);
+            newTriPoints.append(tri[2]);
+        }
+    }
+}
+
+
+void Foam::isoSurface::trimToBox
+(
+    const treeBoundBox& bb,
+    DynamicList<point>& triPoints,
+    DynamicList<label>& triMap      // trimmed to original tris
+)
+{
+    if (debug)
+    {
+        Pout<< "isoSurface : trimming to " << bb << endl;
+    }
+
+    // Generate inwards pointing planes
+    PtrList<plane> planes(6);
+    const pointField pts(bb.treeBoundBox::points());
+    forAll(treeBoundBox::faces, faceI)
+    {
+        const face& f = treeBoundBox::faces[faceI];
+        const vector& n = treeBoundBox::faceNormals[faceI];
+        planes.set(faceI, new plane(pts[f[0]], -n));
+    }
+
+    label nTris = triPoints.size()/3;
+
+    DynamicList<point> newTriPoints(triPoints.size()/16);
+    triMap.setCapacity(nTris/16);
+
+    label vertI = 0;
+    for (label triI = 0; triI < nTris; triI++)
+    {
+        const point& p0 = triPoints[vertI++];
+        const point& p1 = triPoints[vertI++];
+        const point& p2 = triPoints[vertI++];
+
+        label oldNPoints = newTriPoints.size();
+        trimToPlanes
+        (
+            planes,
+            triPointRef(p0, p1, p2),
+            newTriPoints
+        );
+
+        label nCells = (newTriPoints.size()-oldNPoints)/3;
+        for (label i = 0; i < nCells; i++)
+        {
+            triMap.append(triI);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurface : trimmed from " << nTris
+            << " down to " << triMap.size()
+            << " triangles." << endl;
+    }
+
+    triPoints.transfer(newTriPoints);
+}
+
+
+void Foam::isoSurface::trimToBox
+(
+    const treeBoundBox& bb,
+    DynamicList<point>& triPoints,  // new points
+    DynamicList<label>& triMap,     // map from (new) triangle to original
+    labelList& triPointMap,         // map from (new) point to original
+    labelList& interpolatedPoints,  // labels of newly introduced points
+    List<FixedList<label, 3> >& interpolatedOldPoints,// and their interpolation
+    List<FixedList<scalar, 3> >& interpolationWeights
+)
+{
+    const List<point> oldTriPoints(triPoints);
+
+    // Trim triPoints, return map
+    trimToBox(bb, triPoints, triMap);
+
+
+    // Find point correspondence:
+    // - one-to-one for preserved original points (triPointMap)
+    // - interpolation for newly introduced points
+    //   (interpolatedOldPoints)
+    label sz = oldTriPoints.size()/100;
+    DynamicList<label> dynInterpolatedPoints(sz);
+    DynamicList<FixedList<label, 3> > dynInterpolatedOldPoints(sz);
+    DynamicList<FixedList<scalar, 3> > dynInterpolationWeights(sz);
+
+
+    triPointMap.setSize(triPoints.size());
+    forAll(triMap, triI)
+    {
+        label oldTriI = triMap[triI];
+
+        // Find point correspondence. Assumes coordinate is bit-copy.
+        for (label i = 0; i < 3; i++)
+        {
+            label pointI = 3*triI+i;
+            const point& pt = triPoints[pointI];
+
+            // Compare to old-triangle's points
+            label matchPointI = -1;
+            for (label j = 0; j < 3; j++)
+            {
+                label oldPointI = 3*oldTriI+j;
+                if (pt == oldTriPoints[oldPointI])
+                {
+                    matchPointI = oldPointI;
+                    break;
+                }
+            }
+
+            triPointMap[pointI] = matchPointI;
+
+            // If new point: calculate and store interpolation
+            if (matchPointI == -1)
+            {
+                dynInterpolatedPoints.append(pointI);
+
+                FixedList<label, 3> oldPoints;
+                oldPoints[0] = 3*oldTriI;
+                oldPoints[1] = 3*oldTriI+1;
+                oldPoints[2] = 3*oldTriI+2;
+                dynInterpolatedOldPoints.append(oldPoints);
+
+                triPointRef tri(oldTriPoints, oldPoints);
+                scalarList bary;
+                tri.barycentric(pt, bary);
+                FixedList<scalar, 3> weights;
+                weights[0] = bary[0];
+                weights[1] = bary[1];
+                weights[2] = bary[2];
+                dynInterpolationWeights.append(weights);
+            }
+        }
+    }
+
+    interpolatedPoints.transfer(dynInterpolatedPoints);
+    interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
+    interpolationWeights.transfer(dynInterpolationWeights);
+}
+
+
 // Does face use valid vertices?
 bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
 {
@@ -1204,437 +1338,6 @@ bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
 }
 
 
-void Foam::isoSurface::calcAddressing
-(
-    const triSurface& surf,
-    List<FixedList<label, 3> >& faceEdges,
-    labelList& edgeFace0,
-    labelList& edgeFace1,
-    Map<labelList>& edgeFacesRest
-) const
-{
-    const pointField& points = surf.points();
-
-    pointField edgeCentres(3*surf.size());
-    label edgeI = 0;
-    forAll(surf, triI)
-    {
-        const labelledTri& tri = surf[triI];
-        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
-        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
-    }
-
-    pointField mergedCentres;
-    labelList oldToMerged;
-    bool hasMerged = mergePoints
-    (
-        edgeCentres,
-        mergeDistance_,
-        false,
-        oldToMerged,
-        mergedCentres
-    );
-
-    if (debug)
-    {
-        Pout<< "isoSurface : detected "
-            << mergedCentres.size()
-            << " geometric edges on " << surf.size() << " triangles." << endl;
-    }
-
-    if (!hasMerged)
-    {
-        return;
-    }
-
-
-    // Determine faceEdges
-    faceEdges.setSize(surf.size());
-    edgeI = 0;
-    forAll(surf, triI)
-    {
-        faceEdges[triI][0] = oldToMerged[edgeI++];
-        faceEdges[triI][1] = oldToMerged[edgeI++];
-        faceEdges[triI][2] = oldToMerged[edgeI++];
-    }
-
-
-    // Determine edgeFaces
-    edgeFace0.setSize(mergedCentres.size());
-    edgeFace0 = -1;
-    edgeFace1.setSize(mergedCentres.size());
-    edgeFace1 = -1;
-    edgeFacesRest.clear();
-
-    // Overflow edge faces for geometric shared edges that turned
-    // out to be different anyway.
-    EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
-
-    forAll(oldToMerged, oldEdgeI)
-    {
-        label triI = oldEdgeI / 3;
-        label edgeI = oldToMerged[oldEdgeI];
-
-        if (edgeFace0[edgeI] == -1)
-        {
-            // First triangle for edge
-            edgeFace0[edgeI] = triI;
-        }
-        else
-        {
-            //- Check that the two triangles actually topologically
-            //  share an edge
-            const labelledTri& prevTri = surf[edgeFace0[edgeI]];
-            const labelledTri& tri = surf[triI];
-
-            label fp = oldEdgeI % 3;
-
-            edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
-            label prevTriIndex = -1;
-
-            forAll(prevTri, i)
-            {
-                if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
-                {
-                    prevTriIndex = i;
-                    break;
-                }
-            }
-
-            if (prevTriIndex == -1)
-            {
-                // Different edge. Store for later.
-                EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
-
-                if (iter != extraEdgeFaces.end())
-                {
-                    labelList& eFaces = iter();
-                    label sz = eFaces.size();
-                    eFaces.setSize(sz+1);
-                    eFaces[sz] = triI;
-                }
-                else
-                {
-                    extraEdgeFaces.insert(e, labelList(1, triI));
-                }
-            }
-            else if (edgeFace1[edgeI] == -1)
-            {
-                edgeFace1[edgeI] = triI;
-            }
-            else
-            {
-                //WarningIn("orientSurface(triSurface&)")
-                //    << "Edge " << edgeI << " with centre "
-                //    << mergedCentres[edgeI]
-                //    << " used by more than two triangles: "
-                //    << edgeFace0[edgeI] << ", "
-                //    << edgeFace1[edgeI] << " and " << triI << endl;
-                Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
-
-                if (iter != edgeFacesRest.end())
-                {
-                    labelList& eFaces = iter();
-                    label sz = eFaces.size();
-                    eFaces.setSize(sz+1);
-                    eFaces[sz] = triI;
-                }
-                else
-                {
-                    edgeFacesRest.insert(edgeI, labelList(1, triI));
-                }
-            }
-        }
-    }
-
-    // Add extraEdgeFaces
-    edgeI = edgeFace0.size();
-
-    edgeFace0.setSize(edgeI + extraEdgeFaces.size());
-    edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
-
-    forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
-    {
-        const labelList& eFaces = iter();
-
-        // The current edge will become edgeI. Replace all occurrences in
-        // faceEdges
-        forAll(eFaces, i)
-        {
-            label triI = eFaces[i];
-            const labelledTri& tri = surf[triI];
-
-            FixedList<label, 3>& fEdges = faceEdges[triI];
-            forAll(tri, fp)
-            {
-                edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
-                if (e == iter.key())
-                {
-                    fEdges[fp] = edgeI;
-                    break;
-                }
-            }
-        }
-
-
-        // Add face to edgeFaces
-
-        edgeFace0[edgeI] = eFaces[0];
-
-        if (eFaces.size() >= 2)
-        {
-            edgeFace1[edgeI] = eFaces[1];
-
-            if (eFaces.size() > 2)
-            {
-                edgeFacesRest.insert
-                (
-                    edgeI,
-                    SubList<label>(eFaces, eFaces.size()-2, 2)
-                );
-            }
-        }
-
-        edgeI++;
-    }
-}
-
-
-void Foam::isoSurface::walkOrientation
-(
-    const triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const label seedTriI,
-    labelList& flipState
-)
-{
-    // Do walk for consistent orientation.
-    DynamicList<label> changedFaces(surf.size());
-
-    changedFaces.append(seedTriI);
-
-    while (changedFaces.size())
-    {
-        DynamicList<label> newChangedFaces(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label triI = changedFaces[i];
-            const labelledTri& tri = surf[triI];
-            const FixedList<label, 3>& fEdges = faceEdges[triI];
-
-            forAll(fEdges, fp)
-            {
-                label edgeI = fEdges[fp];
-
-                // my points:
-                label p0 = tri[fp];
-                label p1 = tri[tri.fcIndex(fp)];
-
-                label nbrI =
-                (
-                    edgeFace0[edgeI] != triI
-                  ? edgeFace0[edgeI]
-                  : edgeFace1[edgeI]
-                );
-
-                if (nbrI != -1 && flipState[nbrI] == -1)
-                {
-                    const labelledTri& nbrTri = surf[nbrI];
-
-                    // nbr points
-                    label nbrFp = findIndex(nbrTri, p0);
-
-                    if (nbrFp == -1)
-                    {
-                        FatalErrorIn("isoSurface::walkOrientation(..)")
-                            << "triI:" << triI
-                            << " tri:" << tri
-                            << " p0:" << p0
-                            << " p1:" << p1
-                            << " fEdges:" << fEdges
-                            << " edgeI:" << edgeI
-                            << " edgeFace0:" << edgeFace0[edgeI]
-                            << " edgeFace1:" << edgeFace1[edgeI]
-                            << " nbrI:" << nbrI
-                            << " nbrTri:" << nbrTri
-                            << abort(FatalError);
-                    }
-
-
-                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-
-                    bool sameOrientation = (p1 == nbrP1);
-
-                    if (flipState[triI] == 0)
-                    {
-                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-                    }
-                    else
-                    {
-                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-                    }
-                    newChangedFaces.append(nbrI);
-                }
-            }
-        }
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-void Foam::isoSurface::orientSurface
-(
-    triSurface& surf,
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest
-)
-{
-    // -1 : unvisited
-    //  0 : leave as is
-    //  1 : flip
-    labelList flipState(surf.size(), -1);
-
-    label seedTriI = 0;
-
-    while (true)
-    {
-        // Find first unvisited triangle
-        for
-        (
-            ;
-            seedTriI < surf.size() && flipState[seedTriI] != -1;
-            seedTriI++
-        )
-        {}
-
-        if (seedTriI == surf.size())
-        {
-            break;
-        }
-
-        // Note: Determine orientation of seedTriI?
-        // for now assume it is ok
-        flipState[seedTriI] = 0;
-
-        walkOrientation
-        (
-            surf,
-            faceEdges,
-            edgeFace0,
-            edgeFace1,
-            seedTriI,
-            flipState
-        );
-    }
-
-    // Do actual flipping
-    surf.clearOut();
-    forAll(surf, triI)
-    {
-        if (flipState[triI] == 1)
-        {
-            labelledTri tri(surf[triI]);
-
-            surf[triI][0] = tri[0];
-            surf[triI][1] = tri[2];
-            surf[triI][2] = tri[1];
-        }
-        else if (flipState[triI] == -1)
-        {
-            FatalErrorIn
-            (
-                "isoSurface::orientSurface(triSurface&, const label)"
-            )   << "problem" << abort(FatalError);
-        }
-    }
-}
-
-
-// Checks if triangle is connected through edgeI only.
-bool Foam::isoSurface::danglingTriangle
-(
-    const FixedList<label, 3>& fEdges,
-    const labelList& edgeFace1
-)
-{
-    label nOpen = 0;
-    forAll(fEdges, i)
-    {
-        if (edgeFace1[fEdges[i]] == -1)
-        {
-            nOpen++;
-        }
-    }
-
-    if (nOpen == 1 || nOpen == 2 || nOpen == 3)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Mark triangles to keep. Returns number of dangling triangles.
-Foam::label Foam::isoSurface::markDanglingTriangles
-(
-    const List<FixedList<label, 3> >& faceEdges,
-    const labelList& edgeFace0,
-    const labelList& edgeFace1,
-    const Map<labelList>& edgeFacesRest,
-    boolList& keepTriangles
-)
-{
-    keepTriangles.setSize(faceEdges.size());
-    keepTriangles = true;
-
-    label nDangling = 0;
-
-    // Remove any dangling triangles
-    forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
-    {
-        // These are all the non-manifold edges. Filter out all triangles
-        // with only one connected edge (= this edge)
-
-        label edgeI = iter.key();
-        const labelList& otherEdgeFaces = iter();
-
-        // Remove all dangling triangles
-        if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace0[edgeI]] = false;
-            nDangling++;
-        }
-        if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
-        {
-            keepTriangles[edgeFace1[edgeI]] = false;
-            nDangling++;
-        }
-        forAll(otherEdgeFaces, i)
-        {
-            label triI = otherEdgeFaces[i];
-            if (danglingTriangle(faceEdges[triI], edgeFace1))
-            {
-                keepTriangles[triI] = false;
-                nDangling++;
-            }
-        }
-    }
-    return nDangling;
-}
-
-
 Foam::triSurface Foam::isoSurface::subsetMesh
 (
     const triSurface& s,
@@ -1715,6 +1418,7 @@ Foam::isoSurface::isoSurface
     const scalarField& pVals,
     const scalar iso,
     const bool regularise,
+    const boundBox& bounds,
     const scalar mergeTol
 )
 :
@@ -1722,6 +1426,7 @@ Foam::isoSurface::isoSurface
     pVals_(pVals),
     iso_(iso),
     regularise_(regularise),
+    bounds_(bounds),
     mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
     if (debug)
@@ -1957,8 +1662,8 @@ Foam::isoSurface::isoSurface
     }
 
 
-
+    {
-    DynamicList<point> triPoints(nCutCells_);
+        DynamicList<point> triPoints(3*nCutCells_);
         DynamicList<label> triMeshCells(nCutCells_);
 
         generateTriPoints
@@ -1973,8 +1678,8 @@ Foam::isoSurface::isoSurface
             snappedCc,
             snappedPoint,
 
-        triPoints,
+            triPoints,      // 3 points of the triangle
-        triMeshCells
+            triMeshCells    // per triangle the originating cell
         );
 
         if (debug)
@@ -1984,6 +1689,27 @@ Foam::isoSurface::isoSurface
                 << " unmerged points." << endl;
         }
 
+        label nOldPoints = triPoints.size();
+
+        // Trimmed to original triangle
+        DynamicList<label> trimTriMap;
+        // Trimmed to original point
+        labelList trimTriPointMap;
+        if (bounds_ != boundBox::greatBox)
+        {
+            trimToBox
+            (
+                treeBoundBox(bounds_),
+                triPoints,              // new points
+                trimTriMap,             // map from (new) triangle to original
+                trimTriPointMap,        // map from (new) point to original
+                interpolatedPoints_,    // labels of newly introduced points
+                interpolatedOldPoints_, // and their interpolation
+                interpolationWeights_
+            );
+            triMeshCells = labelField(triMeshCells, trimTriMap);
+        }
+
 
         // Merge points and compact out non-valid triangles
         labelList triMap;           // merged to unmerged triangle
@@ -2004,11 +1730,35 @@ Foam::isoSurface::isoSurface
                 << " merged triangles." << endl;
         }
 
+
+        if (bounds_ != boundBox::greatBox)
+        {
+            // Adjust interpolatedPoints_
+            inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
+
+            // Adjust triPointMergeMap_
+            labelList newTriPointMergeMap(nOldPoints, -1);
+            forAll(trimTriPointMap, trimPointI)
+            {
+                label oldPointI = trimTriPointMap[trimPointI];
+                if (oldPointI >= 0)
+                {
+                    label pointI = triPointMergeMap_[trimPointI];
+                    if (pointI >= 0)
+                    {
+                        newTriPointMergeMap[oldPointI] = pointI;
+                    }
+                }
+            }
+            triPointMergeMap_.transfer(newTriPointMergeMap);
+        }
+
         meshCells_.setSize(triMap.size());
         forAll(triMap, i)
         {
             meshCells_[i] = triMeshCells[triMap[i]];
         }
+    }
 
     if (debug)
     {
@@ -2020,75 +1770,7 @@ Foam::isoSurface::isoSurface
             // Copied from surfaceCheck
             validTri(*this, triI);
         }
-    }
 
-
-    if (false)
-    {
-        List<FixedList<label, 3> > faceEdges;
-        labelList edgeFace0, edgeFace1;
-        Map<labelList> edgeFacesRest;
-
-
-        while (true)
-        {
-            // Calculate addressing
-            calcAddressing
-            (
-                *this,
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest
-            );
-
-            // See if any dangling triangles
-            boolList keepTriangles;
-            label nDangling = markDanglingTriangles
-            (
-                faceEdges,
-                edgeFace0,
-                edgeFace1,
-                edgeFacesRest,
-                keepTriangles
-            );
-
-            if (debug)
-            {
-                Pout<< "isoSurface : detected " << nDangling
-                    << " dangling triangles." << endl;
-            }
-
-            if (nDangling == 0)
-            {
-                break;
-            }
-
-            // Create face map (new to old)
-            labelList subsetTriMap(findIndices(keepTriangles, true));
-
-            labelList subsetPointMap;
-            labelList reversePointMap;
-            triSurface::operator=
-            (
-                subsetMesh
-                (
-                    *this,
-                    subsetTriMap,
-                    reversePointMap,
-                    subsetPointMap
-                )
-            );
-            meshCells_ = labelField(meshCells_, subsetTriMap);
-            inplaceRenumber(reversePointMap, triPointMergeMap_);
-        }
-
-        orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-    }
-
-
-    if (debug)
-    {
         fileName stlFile = mesh_.time().path() + ".stl";
         Pout<< "Dumping surface to " << stlFile << endl;
         triSurface::write(stlFile);

src/sampling/sampledSurface/isoSurface/isoSurface.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,8 @@ namespace Foam
 {
 
 class fvMesh;
+class plane;
+class treeBoundBox;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurface Declaration
@@ -116,10 +118,12 @@ class isoSurface
         //- Regularise?
         const Switch regularise_;
 
+        //- Optional bounds
+        const boundBox bounds_;
+
         //- When to merge points
         const scalar mergeDistance_;
 
-
         //- Whether face might be cut
         List<cellCutType> faceCutType_;
 
@@ -135,6 +139,15 @@ class isoSurface
         //- For every unmerged triangle point the point in the triSurface
         labelList triPointMergeMap_;
 
+        //- triSurface points that have weighted interpolation
+        DynamicList<label> interpolatedPoints_;
+
+        //- corresponding original, unmerged points
+        DynamicList<FixedList<label, 3> > interpolatedOldPoints_;
+
+        //- corresponding weights
+        DynamicList<FixedList<scalar, 3> > interpolationWeights_;
+
 
     // Private Member Functions
 
@@ -193,20 +206,8 @@ class isoSurface
             const scalarField& pVals
         );
 
-        static labelPair findCommonPoints
-        (
-            const labelledTri&,
-            const labelledTri&
-        );
-
         static point calcCentre(const triSurface&);
 
-        static pointIndexHit collapseSurface
-        (
-            pointField& localPoints,
-            DynamicList<labelledTri, 64>& localTris
-        );
-
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
@@ -323,6 +324,17 @@ class isoSurface
             DynamicList<label>& triMeshCells
         ) const;
 
+        template<class Type>
+        static tmp<Field<Type> > interpolate
+        (
+            const label nPoints,
+            const labelList& triPointMergeMap,
+            const labelList& interpolatedPoints,
+            const List<FixedList<label, 3> >& interpolatedOldPoints,
+            const List<FixedList<scalar, 3> >& interpolationWeights,
+            const DynamicList<Type>& unmergedValues
+        );
+
         triSurface stitchTriPoints
         (
             const bool checkDuplicates,
@@ -331,57 +343,38 @@ class isoSurface
             labelList& triMap               // merged to unmerged triangle
         ) const;
 
+        //- Trim triangle to planes
+        static void trimToPlanes
+        (
+            const PtrList<plane>& planes,
+            const triPointRef& tri,
+            DynamicList<point>& newTriPoints
+        );
+
+        //- Trim all triangles to box
+        static void trimToBox
+        (
+            const treeBoundBox& bb,
+            DynamicList<point>& triPoints,
+            DynamicList<label>& triMeshCells
+        );
+
+        //- Trim all triangles to box. Determine interpolation
+        //  for existing and new points
+        static void trimToBox
+        (
+            const treeBoundBox& bb,
+            DynamicList<point>& triPoints,
+            DynamicList<label>& triMap,
+            labelList& triPointMap,
+            labelList& interpolatedPoints,
+            List<FixedList<label, 3> >& interpolatedOldPoints,
+            List<FixedList<scalar, 3> >& interpolationWeights
+        );
+
         //- Check single triangle for (topological) validity
         static bool validTri(const triSurface&, const label);
 
-        //- Determine edge-face addressing
-        void calcAddressing
-        (
-            const triSurface& surf,
-            List<FixedList<label, 3> >& faceEdges,
-            labelList& edgeFace0,
-            labelList& edgeFace1,
-            Map<labelList>& edgeFacesRest
-        ) const;
-
-        //- Determine orientation
-        static void walkOrientation
-        (
-            const triSurface& surf,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const label seedTriI,
-            labelList& flipState
-        );
-
-        //- Orient surface
-        static void orientSurface
-        (
-            triSurface&,
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest
-        );
-
-        //- Is triangle (given by 3 edges) not fully connected?
-        static bool danglingTriangle
-        (
-            const FixedList<label, 3>& fEdges,
-            const labelList& edgeFace1
-        );
-
-        //- Mark all non-fully connected triangles
-        static label markDanglingTriangles
-        (
-            const List<FixedList<label, 3> >& faceEdges,
-            const labelList& edgeFace0,
-            const labelList& edgeFace1,
-            const Map<labelList>& edgeFacesRest,
-            boolList& keepTriangles
-        );
-
         static triSurface subsetMesh
         (
             const triSurface& s,
@@ -392,6 +385,10 @@ class isoSurface
 
 public:
 
+    //- Declare friendship with isoSurfaceCell to share some functionality
+    friend class isoSurfaceCell;
+
+
     //- Runtime type information
     TypeName("isoSurface");
 
@@ -407,24 +404,19 @@ public:
             const scalarField& pointIsoVals,
             const scalar iso,
             const bool regularise,
+            const boundBox& bounds = boundBox::greatBox,
             const scalar mergeTol = 1e-6    // fraction of bounding box
         );
 
 
     // Member Functions
 
-        //- For every face original cell in mesh
+        //- For every triangle the original cell in mesh
         const labelList& meshCells() const
         {
             return meshCells_;
         }
 
-        //- For every unmerged triangle point the point in the triSurface
-        const labelList& triPointMergeMap() const
-        {
-            return triPointMergeMap_;
-        }
-
         //- Interpolates cCoords,pCoords. Uses the references to the original
         //  fields used to create the iso surface.
         template<class Type>

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,9 @@ License
 #include "tetMatcher.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "triPoints.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1222,144 +1225,6 @@ void Foam::isoSurfaceCell::calcAddressing
 }
 
 
-//void Foam::isoSurfaceCell::walkOrientation
-//(
-//    const triSurface& surf,
-//    const List<FixedList<label, 3> >& faceEdges,
-//    const labelList& edgeFace0,
-//    const labelList& edgeFace1,
-//    const label seedTriI,
-//    labelList& flipState
-//)
-//{
-//    // Do walk for consistent orientation.
-//    DynamicList<label> changedFaces(surf.size());
-//
-//    changedFaces.append(seedTriI);
-//
-//    while (changedFaces.size())
-//    {
-//        DynamicList<label> newChangedFaces(changedFaces.size());
-//
-//        forAll(changedFaces, i)
-//        {
-//            label triI = changedFaces[i];
-//            const labelledTri& tri = surf[triI];
-//            const FixedList<label, 3>& fEdges = faceEdges[triI];
-//
-//            forAll(fEdges, fp)
-//            {
-//                label edgeI = fEdges[fp];
-//
-//                // my points:
-//                label p0 = tri[fp];
-//                label p1 = tri[tri.fcIndex(fp)];
-//
-//                label nbrI =
-//                (
-//                    edgeFace0[edgeI] != triI
-//                  ? edgeFace0[edgeI]
-//                  : edgeFace1[edgeI]
-//                );
-//
-//                if (nbrI != -1 && flipState[nbrI] == -1)
-//                {
-//                    const labelledTri& nbrTri = surf[nbrI];
-//
-//                    // nbr points
-//                    label nbrFp = findIndex(nbrTri, p0);
-//                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-//
-//                    bool sameOrientation = (p1 == nbrP1);
-//
-//                    if (flipState[triI] == 0)
-//                    {
-//                        flipState[nbrI] = (sameOrientation ? 0 : 1);
-//                    }
-//                    else
-//                    {
-//                        flipState[nbrI] = (sameOrientation ? 1 : 0);
-//                    }
-//                    newChangedFaces.append(nbrI);
-//                }
-//            }
-//        }
-//
-//        changedFaces.transfer(newChangedFaces);
-//    }
-//}
-//
-//
-//void Foam::isoSurfaceCell::orientSurface
-//(
-//    triSurface& surf,
-//    const List<FixedList<label, 3> >& faceEdges,
-//    const labelList& edgeFace0,
-//    const labelList& edgeFace1,
-//    const Map<labelList>& edgeFacesRest
-//)
-//{
-//    // -1 : unvisited
-//    //  0 : leave as is
-//    //  1 : flip
-//    labelList flipState(surf.size(), -1);
-//
-//    label seedTriI = 0;
-//
-//    while (true)
-//    {
-//        // Find first unvisited triangle
-//        for
-//        (
-//            ;
-//            seedTriI < surf.size() && flipState[seedTriI] != -1;
-//            seedTriI++
-//        )
-//        {}
-//
-//        if (seedTriI == surf.size())
-//        {
-//            break;
-//        }
-//
-//        // Note: Determine orientation of seedTriI?
-//        // for now assume it is ok
-//        flipState[seedTriI] = 0;
-//
-//        walkOrientation
-//        (
-//            surf,
-//            faceEdges,
-//            edgeFace0,
-//            edgeFace1,
-//            seedTriI,
-//            flipState
-//        );
-//    }
-//
-//    // Do actual flipping
-//    surf.clearOut();
-//    forAll(surf, triI)
-//    {
-//        if (flipState[triI] == 1)
-//        {
-//            labelledTri tri(surf[triI]);
-//
-//            surf[triI][0] = tri[0];
-//            surf[triI][1] = tri[2];
-//            surf[triI][2] = tri[1];
-//        }
-//        else if (flipState[triI] == -1)
-//        {
-//            FatalErrorIn
-//            (
-//                "isoSurfaceCell::orientSurface(triSurface&, const label)"
-//            )   << "problem" << abort(FatalError);
-//        }
-//    }
-//}
-
-
 // Checks if triangle is connected through edgeI only.
 bool Foam::isoSurfaceCell::danglingTriangle
 (
@@ -1517,6 +1382,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     const scalarField& pVals,
     const scalar iso,
     const bool regularise,
+    const boundBox& bounds,
     const scalar mergeTol
 )
 :
@@ -1524,6 +1390,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     cVals_(cVals),
     pVals_(pVals),
     iso_(iso),
+    bounds_(bounds),
     mergeDistance_(mergeTol*mesh.bounds().mag())
 {
     if (debug)
@@ -1607,7 +1474,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
-
+    {
         DynamicList<point> triPoints(nCutCells_);
         DynamicList<label> triMeshCells(nCutCells_);
 
@@ -1633,8 +1500,32 @@ Foam::isoSurfaceCell::isoSurfaceCell
                 << " unmerged triangles." << endl;
         }
 
+
+        label nOldPoints = triPoints.size();
+
+        // Trimmed to original triangle
+        DynamicList<label> trimTriMap;
+        // Trimmed to original point
+        labelList trimTriPointMap;
+        if (bounds_ != boundBox::greatBox)
+        {
+            isoSurface::trimToBox
+            (
+                treeBoundBox(bounds_),
+                triPoints,              // new points
+                trimTriMap,             // map from (new) triangle to original
+                trimTriPointMap,        // map from (new) point to original
+                interpolatedPoints_,    // labels of newly introduced points
+                interpolatedOldPoints_, // and their interpolation
+                interpolationWeights_
+            );
+            triMeshCells = labelField(triMeshCells, trimTriMap);
+        }
+
+
+
         // Merge points and compact out non-valid triangles
-    labelList triMap;           // merged to unmerged triangle
+        labelList triMap;
         triSurface::operator=
         (
             stitchTriPoints
@@ -1642,7 +1533,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
                 regularise,         // check for duplicate tris
                 triPoints,
                 triPointMergeMap_,  // unmerged to merged point
-            triMap
+                triMap              // merged to unmerged triangle
             )
         );
 
@@ -1652,11 +1543,35 @@ Foam::isoSurfaceCell::isoSurfaceCell
                 << " merged triangles." << endl;
         }
 
+        if (bounds_ != boundBox::greatBox)
+        {
+            // Adjust interpolatedPoints_
+            inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
+
+            // Adjust triPointMergeMap_
+            labelList newTriPointMergeMap(nOldPoints, -1);
+            forAll(trimTriPointMap, trimPointI)
+            {
+                label oldPointI = trimTriPointMap[trimPointI];
+                if (oldPointI >= 0)
+                {
+                    label pointI = triPointMergeMap_[trimPointI];
+                    if (pointI >= 0)
+                    {
+                        newTriPointMergeMap[oldPointI] = pointI;
+                    }
+                }
+            }
+            triPointMergeMap_.transfer(newTriPointMergeMap);
+        }
+
         meshCells_.setSize(triMap.size());
         forAll(triMap, i)
         {
             meshCells_[i] = triMeshCells[triMap[i]];
         }
+    }
+
 
     if (debug)
     {
@@ -1730,8 +1645,6 @@ Foam::isoSurfaceCell::isoSurfaceCell
             meshCells_ = labelField(meshCells_, subsetTriMap);
             inplaceRenumber(reversePointMap, triPointMergeMap_);
         }
-
-        //orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
     }
 }
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,7 @@ SourceFiles
 #include "labelPair.H"
 #include "pointIndexHit.H"
 #include "PackedBoolList.H"
+#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,6 +56,7 @@ namespace Foam
 {
 
 class polyMesh;
+class plane;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurfaceCell Declaration
@@ -91,6 +93,9 @@ class isoSurfaceCell
         //- isoSurfaceCell value
         const scalar iso_;
 
+        //- Optional bounds
+        const boundBox bounds_;
+
         //- When to merge points
         const scalar mergeDistance_;
 
@@ -106,6 +111,15 @@ class isoSurfaceCell
         //- For every unmerged triangle point the point in the triSurface
         labelList triPointMergeMap_;
 
+        //- triSurface points that have weighted interpolation
+        DynamicList<label> interpolatedPoints_;
+
+        //- corresponding original, unmerged points
+        DynamicList<FixedList<label, 3> > interpolatedOldPoints_;
+
+        //- corresponding weights
+        DynamicList<FixedList<scalar, 3> > interpolationWeights_;
+
 
     // Private Member Functions
 
@@ -214,19 +228,15 @@ class isoSurfaceCell
         void generateTriPoints
         (
             const DynamicList<Type>& snapped,
-            const scalar isoVal0,
             const scalar s0,
             const Type& p0,
             const label p0Index,
-            const scalar isoVal1,
             const scalar s1,
             const Type& p1,
             const label p1Index,
-            const scalar isoVal2,
             const scalar s2,
             const Type& p2,
             const label p2Index,
-            const scalar isoVal3,
             const scalar s3,
             const Type& p3,
             const label p3Index,
@@ -271,27 +281,6 @@ class isoSurfaceCell
             Map<labelList>& edgeFacesRest
         ) const;
 
-        ////- Determine orientation
-        //static void walkOrientation
-        //(
-        //    const triSurface& surf,
-        //    const List<FixedList<label, 3> >& faceEdges,
-        //    const labelList& edgeFace0,
-        //    const labelList& edgeFace1,
-        //    const label seedTriI,
-        //    labelList& flipState
-        //);
-
-        ////- Orient surface
-        //static void orientSurface
-        //(
-        //    triSurface&,
-        //    const List<FixedList<label, 3> >& faceEdges,
-        //    const labelList& edgeFace0,
-        //    const labelList& edgeFace1,
-        //    const Map<labelList>& edgeFacesRest
-        //);
-
         //- Is triangle (given by 3 edges) not fully connected?
         static bool danglingTriangle
         (
@@ -317,8 +306,6 @@ class isoSurfaceCell
             labelList& newToOldPoints
         );
 
-        //- Combine all triangles inside a cell into a minimal triangulation
-        void combineCellTriangles();
 
 public:
 
@@ -336,6 +323,7 @@ public:
             const scalarField& pointValues,
             const scalar iso,
             const bool regularise,
+            const boundBox& bounds = boundBox::greatBox,
             const scalar mergeTol = 1e-6    // fraction of bounding box
         );
 
@@ -348,24 +336,6 @@ public:
             return meshCells_;
         }
 
-        //- For every unmerged triangle point the point in the triSurface
-        const labelList triPointMergeMap() const
-        {
-            return triPointMergeMap_;
-        }
-
-
-        //- Interpolates cCoords,pCoords. Takes the original fields
-        //  used to create the iso surface.
-        template<class Type>
-        tmp<Field<Type> > interpolate
-        (
-            const scalarField& cVals,
-            const scalarField& pVals,
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
-
         //- Interpolates cCoords,pCoords.
         template<class Type>
         tmp<Field<Type> > interpolate

src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C

@@ -26,6 +26,7 @@ License
 #include "isoSurfaceCell.H"
 #include "polyMesh.H"
 #include "tetMatcher.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -76,22 +77,18 @@ void Foam::isoSurfaceCell::generateTriPoints
 (
     const DynamicList<Type>& snapped,
 
-    const scalar isoVal0,
     const scalar s0,
     const Type& p0,
     const label p0Index,
 
-    const scalar isoVal1,
     const scalar s1,
     const Type& p1,
     const label p1Index,
 
-    const scalar isoVal2,
     const scalar s2,
     const Type& p2,
     const label p2Index,
 
-    const scalar isoVal3,
     const scalar s3,
     const Type& p3,
     const label p3Index,
@@ -124,160 +121,196 @@ void Foam::isoSurfaceCell::generateTriPoints
         case 0x0F:
         break;
 
-        case 0x0E:
         case 0x01:
+        case 0x0E:
         {
-            // 0 is common point. Orient such that normal points in positive
+            pts.append
-            // gradient direction
+            (
-            if (isoVal0 >= isoVal1)
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+            if (triIndex == 0x0E)
             {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+                // Flip normals
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+                label sz = pts.size();
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                Swap(pts[sz-2], pts[sz-1]);
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
             }
         }
         break;
 
-        case 0x0D:
         case 0x02:
+        case 0x0D:
         {
-            // 1 is common point
+            pts.append
-            if (isoVal1 >= isoVal0)
+            (
+                generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+
+            if (triIndex == 0x0D)
             {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+                // Flip normals
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                label sz = pts.size();
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+                Swap(pts[sz-2], pts[sz-1]);
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
             }
         }
         break;
 
-        case 0x0C:
         case 0x03:
+        case 0x0C:
         {
-            Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type p0p2 =
-            Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type p1p3 =
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
 
-            if (isoVal0 >= isoVal3)
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+            pts.append(p1p3);
+            pts.append(p0p2);
+
+            pts.append(p1p3);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+            pts.append(p0p2);
+
+            if (triIndex == 0x0C)
             {
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                // Flip normals
-                pts.append(s02);
+                label sz = pts.size();
-                pts.append(s13);
+                Swap(pts[sz-5], pts[sz-4]);
-                pts.append(s13);
+                Swap(pts[sz-2], pts[sz-1]);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s02);
-            }
-            else
-            {
-                pts.append(s02);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(s13);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s13);
-                pts.append(s02);
             }
         }
         break;
 
-        case 0x0B:
         case 0x04:
+        case 0x0B:
         {
-            // 2 is common point
+            pts.append
-            if (isoVal2 >= isoVal0)
+            (
+                generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)
+            );
+
+            if (triIndex == 0x0B)
             {
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+                // Flip normals
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+                label sz = pts.size();
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+                Swap(pts[sz-2], pts[sz-1]);
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
             }
         }
         break;
 
-        case 0x0A:
         case 0x05:
+        case 0x0A:
         {
-            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p0p1 =
-            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p2p3 =
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            if (isoVal3 >= isoVal0)
+            pts.append(p0p1);
+            pts.append(p2p3);
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+            );
+
+            pts.append(p0p1);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+            );
+            pts.append(p2p3);
+
+            if (triIndex == 0x0A)
             {
-                pts.append(s01);
+                // Flip normals
-                pts.append(s23);
+                label sz = pts.size();
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+                Swap(pts[sz-5], pts[sz-4]);
-                pts.append(s01);
+                Swap(pts[sz-2], pts[sz-1]);
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s23);
-            }
-            else
-            {
-                pts.append(s23);
-                pts.append(s01);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-                pts.append(s01);
-                pts.append(s23);
             }
         }
         break;
 
-        case 0x09:
         case 0x06:
+        case 0x09:
         {
-            Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p0p1 =
-            Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+                generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type p2p3 =
+                generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
 
-            if (isoVal3 >= isoVal1)
+            pts.append(p0p1);
+            pts.append
+            (
+                generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+            );
+            pts.append(p2p3);
+
+            pts.append(p0p1);
+            pts.append(p2p3);
+            pts.append
+            (
+                generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+            );
+
+            if (triIndex == 0x09)
             {
-                pts.append(s01);
+                // Flip normals
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+                label sz = pts.size();
-                pts.append(s23);
+                Swap(pts[sz-5], pts[sz-4]);
-                pts.append(s01);
+                Swap(pts[sz-2], pts[sz-1]);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(s23);
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-                pts.append(s01);
-                pts.append(s23);
-                pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-                pts.append(s01);
-                pts.append(s23);
             }
         }
         break;
 
-        case 0x07:
         case 0x08:
+        case 0x07:
         {
-            // 3 is common point
+            pts.append
-            if (isoVal3 >= isoVal0)
+            (
+                generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)
+            );
+            pts.append
+            (
+                generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)
+            );
+            if (triIndex == 0x07)
             {
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+                // Flip normals
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+                label sz = pts.size();
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+                Swap(pts[sz-2], pts[sz-1]);
-            }
-            else
-            {
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-                pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
             }
         }
         break;
@@ -341,22 +374,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                     (
                         snappedPoints,
 
-                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
-                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
-                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
-                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -370,22 +399,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                     (
                         snappedPoints,
 
-                        pVals_[f0[0]],
                         pVals[f0[0]],
                         pCoords[f0[0]],
                         snappedPoint[f0[0]],
 
-                        pVals_[f0[1]],
                         pVals[f0[1]],
                         pCoords[f0[1]],
                         snappedPoint[f0[1]],
 
-                        pVals_[f0[2]],
                         pVals[f0[2]],
                         pCoords[f0[2]],
                         snappedPoint[f0[2]],
 
-                        pVals_[oppositeI],
                         pVals[oppositeI],
                         pCoords[oppositeI],
                         snappedPoint[oppositeI],
@@ -411,7 +436,6 @@ void Foam::isoSurfaceCell::generateTriPoints
                     }
 
                     label fp = f.fcIndex(fp0);
-
                     for (label i = 2; i < f.size(); i++)
                     {
                         label nextFp = f.fcIndex(fp);
@@ -425,22 +449,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
-                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
-                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
-                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
-                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -454,22 +474,18 @@ void Foam::isoSurfaceCell::generateTriPoints
                             (
                                 snappedPoints,
 
-                                pVals_[tri[0]],
                                 pVals[tri[0]],
                                 pCoords[tri[0]],
                                 snappedPoint[tri[0]],
 
-                                pVals_[tri[1]],
                                 pVals[tri[1]],
                                 pCoords[tri[1]],
                                 snappedPoint[tri[1]],
 
-                                pVals_[tri[2]],
                                 pVals[tri[2]],
                                 pCoords[tri[2]],
                                 snappedPoint[tri[2]],
 
-                                cVals_[cellI],
                                 cVals[cellI],
                                 cCoords[cellI],
                                 snappedCc[cellI],
@@ -495,7 +511,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 
     if (countNotFoundTets > 0)
     {
-        WarningIn("Foam::isoSurfaceCell::generateTriPoints")
+        WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
             << "Could not find " << countNotFoundTets
             << " tet base points, which may lead to inverted triangles."
             << endl;
@@ -510,13 +526,11 @@ template<class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::isoSurfaceCell::interpolate
 (
-    const scalarField& cVals,
-    const scalarField& pVals,
     const Field<Type>& cCoords,
     const Field<Type>& pCoords
 ) const
 {
-    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<Type> triPoints(3*nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     // Dummy snap data
@@ -524,59 +538,6 @@ Foam::isoSurfaceCell::interpolate
     labelList snappedCc(mesh_.nCells(), -1);
     labelList snappedPoint(mesh_.nPoints(), -1);
 
-
-    generateTriPoints
-    (
-        cVals,
-        pVals,
-
-        cCoords,
-        pCoords,
-
-        snappedPoints,
-        snappedCc,
-        snappedPoint,
-
-        triPoints,
-        triMeshCells
-    );
-
-
-    // One value per point
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] = triPoints[i];
-        }
-    }
-
-    return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type> >
-Foam::isoSurfaceCell::interpolate
-(
-    const Field<Type>& cCoords,
-    const Field<Type>& pCoords
-) const
-{
-    DynamicList<Type> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells_);
-
-    // Dummy snap data
-    DynamicList<Type> snappedPoints;
-    labelList snappedCc(mesh_.nCells(), -1);
-    labelList snappedPoint(mesh_.nPoints(), -1);
-
-
     generateTriPoints
     (
         cVals_,
@@ -593,22 +554,15 @@ Foam::isoSurfaceCell::interpolate
         triMeshCells
     );
 
-
+    return isoSurface::interpolate
-    // One value per point
+    (
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+        points().size(),
-    Field<Type>& values = tvalues();
+        triPointMergeMap_,
-
+        interpolatedPoints_,
-    forAll(triPoints, i)
+        interpolatedOldPoints_,
-    {
+        interpolationWeights_,
-        label mergedPointI = triPointMergeMap_[i];
+        triPoints
-
+    );
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] = triPoints[i];
-        }
-    }
-
-    return tvalues;
 }
 
 

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -120,7 +120,7 @@ Foam::isoSurface::adaptPatchFields
         {
             fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
             (
-                fld.boundaryField()[patchI]
+                sliceFld.boundaryField()[patchI]
             );
 
             const scalarField& w = mesh.weights().boundaryField()[patchI];
@@ -189,6 +189,8 @@ Type Foam::isoSurface::generatePoint
 }
 
 
+// Note: cannot use simpler isoSurfaceCell::generateTriPoints since
+//       the need here to sometimes pass in remote 'snappedPoints'
 template<class Type>
 void Foam::isoSurface::generateTriPoints
 (
@@ -240,8 +242,8 @@ void Foam::isoSurface::generateTriPoints
         case 0x0F:
         break;
 
-        case 0x0E:
         case 0x01:
+        case 0x0E:
             points.append
             (
                 generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
@@ -254,10 +256,16 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
             );
+            if (triIndex == 0x0E)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0D:
         case 0x02:
+        case 0x0D:
             points.append
             (
                 generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
@@ -270,33 +278,47 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
             );
+            if (triIndex == 0x0D)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0C:
         case 0x03:
+        case 0x0C:
             {
-            Type tp1 =
+                Type p0p2 =
                     generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
-            Type tp2 =
+                Type p1p3 =
                     generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
 
                 points.append
                 (
                     generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
                 );
-            points.append(tp1);
+                points.append(p1p3);
-            points.append(tp2);
+                points.append(p0p2);
-            points.append(tp2);
+
+                points.append(p1p3);
                 points.append
                 (
                     generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
                 );
-            points.append(tp1);
+                points.append(p0p2);
+            }
+            if (triIndex == 0x0C)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
             }
         break;
 
-        case 0x0B:
         case 0x04:
+        case 0x0B:
             {
                 points.append
                 (
@@ -311,56 +333,78 @@ void Foam::isoSurface::generateTriPoints
                     generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
                 );
             }
+            if (triIndex == 0x0B)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
 
-        case 0x0A:
         case 0x05:
+        case 0x0A:
             {
-            Type tp0 =
+                Type p0p1 =
                     generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
-            Type tp1 =
+                Type p2p3 =
                     generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
-            points.append(tp0);
+                points.append(p0p1);
-            points.append(tp1);
+                points.append(p2p3);
                 points.append
                 (
                     generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
                 );
-            points.append(tp0);
+
+                points.append(p0p1);
                 points.append
                 (
                     generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
                 );
-            points.append(tp1);
+                points.append(p2p3);
+            }
+            if (triIndex == 0x0A)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
             }
         break;
 
-        case 0x09:
         case 0x06:
+        case 0x09:
             {
-            Type tp0 =
+                Type p0p1 =
                     generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
-            Type tp1 =
+                Type p2p3 =
                     generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
-            points.append(tp0);
+                points.append(p0p1);
                 points.append
                 (
                     generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
                 );
-            points.append(tp1);
+                points.append(p2p3);
-            points.append(tp0);
+
+                points.append(p0p1);
+                points.append(p2p3);
                 points.append
                 (
                     generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
                 );
-            points.append(tp1);
+            }
+            if (triIndex == 0x09)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-5], points[sz-4]);
+                Swap(points[sz-2], points[sz-1]);
             }
         break;
 
-        case 0x07:
         case 0x08:
+        case 0x07:
             points.append
             (
                 generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
@@ -373,6 +417,12 @@ void Foam::isoSurface::generateTriPoints
             (
                 generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
             );
+            if (triIndex == 0x07)
+            {
+                // Flip normals
+                label sz = points.size();
+                Swap(points[sz-2], points[sz-1]);
+            }
         break;
     }
 }
@@ -689,6 +739,69 @@ void Foam::isoSurface::generateTriPoints
 //}
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurface::interpolate
+(
+    const label nPoints,
+    const labelList& triPointMergeMap,
+    const labelList& interpolatedPoints,
+    const List<FixedList<label, 3> >& interpolatedOldPoints,
+    const List<FixedList<scalar, 3> >& interpolationWeights,
+    const DynamicList<Type>& unmergedValues
+)
+{
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(nPoints, pTraits<Type>::zero));
+    Field<Type>& values = tvalues();
+
+
+    // Pass1: unweighted average of merged point values
+    {
+        labelList nValues(values.size(), 0);
+
+        forAll(unmergedValues, i)
+        {
+            label mergedPointI = triPointMergeMap[i];
+
+            if (mergedPointI >= 0)
+            {
+                values[mergedPointI] += unmergedValues[i];
+                nValues[mergedPointI]++;
+            }
+        }
+
+        forAll(values, i)
+        {
+            if (nValues[i] > 0)
+            {
+                values[i] /= scalar(nValues[i]);
+            }
+        }
+    }
+
+
+    // Pass2: weighted average for remaining values (from clipped triangles)
+
+    forAll(interpolatedPoints, i)
+    {
+        label pointI = interpolatedPoints[i];
+        const FixedList<label, 3>& oldPoints = interpolatedOldPoints[i];
+        const FixedList<scalar, 3>& w = interpolationWeights[i];
+
+        // Note: zeroing should not be necessary if interpolation only done
+        //       for newly introduced points (i.e. not in triPointMergeMap)
+        values[pointI] = pTraits<Type>::zero;
+        forAll(oldPoints, j)
+        {
+            values[pointI] = w[j]*unmergedValues[oldPoints[j]];
+        }
+    }
+
+    return tvalues;
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::isoSurface::interpolate
@@ -707,7 +820,7 @@ Foam::isoSurface::interpolate
     > > c2(adaptPatchFields(cCoords));
 
 
-    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<Type> triPoints(3*nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
     // Dummy snap data
@@ -731,52 +844,15 @@ Foam::isoSurface::interpolate
         triMeshCells
     );
 
-
+    return interpolate
-    // One value per point
-    tmp<Field<Type> > tvalues
     (
-        new Field<Type>(points().size(), pTraits<Type>::zero)
+        points().size(),
+        triPointMergeMap_,
+        interpolatedPoints_,
+        interpolatedOldPoints_,
+        interpolationWeights_,
+        triPoints
     );
-    Field<Type>& values = tvalues();
-    labelList nValues(values.size(), 0);
-
-    forAll(triPoints, i)
-    {
-        label mergedPointI = triPointMergeMap_[i];
-
-        if (mergedPointI >= 0)
-        {
-            values[mergedPointI] += triPoints[i];
-            nValues[mergedPointI]++;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "nValues:" << values.size() << endl;
-        label nMult = 0;
-        forAll(nValues, i)
-        {
-            if (nValues[i] == 0)
-            {
-                FatalErrorIn("isoSurface::interpolate(..)")
-                    << "point:" << i << " nValues:" << nValues[i]
-                    << abort(FatalError);
-            }
-            else if (nValues[i] > 1)
-            {
-                nMult++;
-            }
-        }
-        Pout<< "Of which mult:" << nMult << endl;
-    }
-
-    forAll(values, i)
-    {
-        values[i] /= scalar(nValues[i]);
-    }
-
-    return tvalues;
 }
 
 

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -122,18 +122,52 @@ void Foam::sampledIsoSurface::getIsoFields() const
     // Get pointField
     // ~~~~~~~~~~~~~~
 
+    // In case of multiple iso values we don't want to calculate multiple e.g.
+    // "volPointInterpolate(p)" so register it and re-use it. This is the
+    // same as the 'cache' functionality from volPointInterpolate but
+    // unfortunately that one does not guarantee that the field pointer
+    // remain: e.g. some other
+    // functionObject might delete the cached version.
+    // (volPointInterpolation::interpolate with cache=false deletes any
+    //  registered one or if mesh.changing())
+
     if (!subMeshPtr_.valid())
     {
-        word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+        const word pointFldName =
+            "volPointInterpolate_"
+          + type()
+          + "("
+          + isoField_
+          + ')';
 
         if (fvm.foundObject<pointScalarField>(pointFldName))
         {
             if (debug)
             {
-                Info<< "sampledIsoSurface::getIsoFields() : lookup pointField "
+                Info<< "sampledIsoSurface::getIsoFields() :"
-                    << pointFldName << endl;
+                    << " lookup pointField " << pointFldName << endl;
             }
-            pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+            const pointScalarField& pfld = fvm.lookupObject<pointScalarField>
+            (
+                pointFldName
+            );
+
+            if (!pfld.upToDate(*volFieldPtr_))
+            {
+                if (debug)
+                {
+                    Info<< "sampledIsoSurface::getIsoFields() :"
+                        << " updating pointField " << pointFldName << endl;
+                }
+                // Update the interpolated value
+                volPointInterpolation::New(fvm).interpolate
+                (
+                    *volFieldPtr_,
+                    const_cast<pointScalarField&>(pfld)
+                );
+            }
+
+            pointFieldPtr_ = &pfld;
         }
         else
         {
@@ -142,32 +176,22 @@ void Foam::sampledIsoSurface::getIsoFields() const
             if (debug)
             {
                 Info<< "sampledIsoSurface::getIsoFields() :"
-                    << "checking pointField " << pointFldName
+                    << " creating and storing pointField "
-                    << " for same time " << fvm.time().timeName()
+                    << pointFldName << " for time "
-                    << endl;
+                    << fvm.time().timeName() << endl;
             }
 
-            if
+            tmp<pointScalarField> tpfld
             (
-                storedPointFieldPtr_.empty()
+                volPointInterpolation::New(fvm).interpolate
-             || (fvm.time().timeName() != storedPointFieldPtr_().instance())
+                (
+                    *volFieldPtr_,
+                    pointFldName,
+                    false
                 )
-            {
-                if (debug)
-                {
-                    Info<< "sampledIsoSurface::getIsoFields() :"
-                        << " interpolating volField " << volFieldPtr_->name()
-                        << " to get pointField " << pointFldName << endl;
-                }
-
-                storedPointFieldPtr_.reset
-                (
-                    volPointInterpolation::New(fvm)
-                    .interpolate(*volFieldPtr_).ptr()
             );
-                storedPointFieldPtr_->checkOut();
+            pointFieldPtr_ = tpfld.ptr();
-                pointFieldPtr_ = storedPointFieldPtr_.operator->();
+            const_cast<pointScalarField*>(pointFieldPtr_)->store();
-            }
         }
 
 
@@ -233,8 +257,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
 
         // Pointfield on submesh
 
-        word pointFldName =
+        const word pointFldName =
-            "volPointInterpolate("
+            "volPointInterpolate_"
+          + type()
+          + "("
           + volSubFieldPtr_->name()
           + ')';
 
@@ -245,11 +271,28 @@ void Foam::sampledIsoSurface::getIsoFields() const
                 Info<< "sampledIsoSurface::getIsoFields() :"
                     << " submesh lookup pointField " << pointFldName << endl;
             }
-            storedPointSubFieldPtr_.clear();
+            const pointScalarField& pfld = subFvm.lookupObject<pointScalarField>
-            pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
             (
                 pointFldName
             );
+
+            if (!pfld.upToDate(*volSubFieldPtr_))
+            {
+                if (debug)
+                {
+                    Info<< "sampledIsoSurface::getIsoFields() :"
+                        << " updating submesh pointField "
+                        << pointFldName << endl;
+                }
+                // Update the interpolated value
+                volPointInterpolation::New(subFvm).interpolate
+                (
+                    *volSubFieldPtr_,
+                    const_cast<pointScalarField&>(pfld)
+                );
+            }
+
+            pointFieldPtr_ = &pfld;
         }
         else
         {
@@ -260,15 +303,15 @@ void Foam::sampledIsoSurface::getIsoFields() const
                     << volSubFieldPtr_->name()
                     << " to get submesh pointField " << pointFldName << endl;
             }
-            storedPointSubFieldPtr_.reset
+            tmp<pointScalarField> tpfld
             (
                 volPointInterpolation::New
                 (
                     subFvm
-                ).interpolate(*volSubFieldPtr_).ptr()
+                ).interpolate(*volSubFieldPtr_)
             );
-            storedPointSubFieldPtr_->checkOut();
+            pointSubFieldPtr_ = tpfld.ptr();
-            pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
+            const_cast<pointScalarField*>(pointSubFieldPtr_)->store();
         }
 
 
@@ -349,28 +392,34 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
     if (subMeshPtr_.valid())
     {
+        const volScalarField& vfld = *volSubFieldPtr_;
+
         surfPtr_.reset
         (
             new isoSurface
             (
-                *volSubFieldPtr_,
+                vfld,
                 *pointSubFieldPtr_,
                 isoVal_,
                 regularise_,
+                bounds_,
                 mergeTol_
             )
         );
     }
     else
     {
+        const volScalarField& vfld = *volFieldPtr_;
+
         surfPtr_.reset
         (
             new isoSurface
             (
-                *volFieldPtr_,
+                vfld,
                 *pointFieldPtr_,
                 isoVal_,
                 regularise_,
+                bounds_,
                 mergeTol_
             )
         );
@@ -412,6 +461,7 @@ Foam::sampledIsoSurface::sampledIsoSurface
     sampledSurface(name, mesh, dict),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
@@ -422,7 +472,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
     prevTimeIndex_(-1),
     storedVolFieldPtr_(NULL),
     volFieldPtr_(NULL),
-    storedPointFieldPtr_(NULL),
     pointFieldPtr_(NULL)
 {
     if (!sampledSurface::interpolate())

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H

@@ -63,6 +63,9 @@ class sampledIsoSurface
         //- Iso value
         const scalar isoVal_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- Merge tolerance
         const scalar mergeTol_;
 
@@ -94,7 +97,6 @@ class sampledIsoSurface
             mutable const volScalarField* volFieldPtr_;
 
             //- Cached pointfield
-            mutable autoPtr<pointScalarField> storedPointFieldPtr_;
             mutable const pointScalarField* pointFieldPtr_;
 
             // And on subsetted mesh
@@ -107,7 +109,6 @@ class sampledIsoSurface
                 mutable const volScalarField* volSubFieldPtr_;
 
                 //- Cached pointfield
-                mutable autoPtr<pointScalarField> storedPointSubFieldPtr_;
                 mutable const pointScalarField* pointSubFieldPtr_;
 
 

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,7 +148,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             cellAvg,
             pointFld().internalField(),
             isoVal_,
-            regularise_
+            regularise_,
+            bounds_
         );
 
         const_cast<sampledIsoSurfaceCell&>
@@ -166,7 +167,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             cellFld.internalField(),
             pointFld().internalField(),
             isoVal_,
-            regularise_
+            regularise_,
+            bounds_
         );
 
         const_cast<sampledIsoSurfaceCell&>
@@ -185,6 +187,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             << "    average        : " << average_ << nl
             << "    isoField       : " << isoField_ << nl
             << "    isoValue       : " << isoVal_ << nl
+            << "    bounds         : " << bounds_ << nl
             << "    points         : " << points().size() << nl
             << "    tris           : " << triSurface::size() << nl
             << "    cut cells      : " << meshCells_.size() << endl;
@@ -206,6 +209,7 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
     sampledSurface(name, mesh, dict),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", true)),
     zoneKey_(keyType::null),

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H

@@ -62,6 +62,9 @@ class sampledIsoSurfaceCell
         //- Iso value
         const scalar isoVal_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- Whether to coarse
         const Switch regularise_;
 

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -230,6 +230,7 @@ void Foam::sampledCuttingPlane::createGeometry()
             pointDistance_,
             0.0,
             regularise_,
+            bounds_,
             mergeTol_
         )
         //new isoSurfaceCell
@@ -262,6 +263,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
 :
     sampledSurface(name, mesh, dict),
     plane_(dict),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -60,6 +60,9 @@ class sampledCuttingPlane
         //- Plane
         const plane plane_;
 
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
+
         //- Merge tolerance
         const scalar mergeTol_;
 

src/sampling/sampledSurface/sampledPatch/sampledPatch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,11 +84,7 @@ const Foam::labelList& Foam::sampledPatch::patchIDs() const
 {
     if (patchIDs_.empty())
     {
-        patchIDs_ = mesh().boundaryMesh().patchSet
+        patchIDs_ = mesh().boundaryMesh().patchSet(patchNames_).sortedToc();
-        (
-            patchNames_,
-            false
-        ).sortedToc();
     }
     return patchIDs_;
 }

src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -162,19 +162,7 @@ void Foam::sampledSurfaces::write()
 
         const label nFields = classifyFields();
 
-        if (Pstream::master())
+        // write geometry first if required,
-        {
-            if (debug)
-            {
-                Pout<< "Creating directory "
-                    << outputPath_/obr_.time().timeName() << nl << endl;
-
-            }
-
-            mkDir(outputPath_/obr_.time().timeName());
-        }
-
-        // Write geometry first if required,
         // or when no fields would otherwise be written
         if (nFields == 0 || formatter_->separateGeometry())
         {

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -644,8 +644,25 @@ bool Foam::sampledTriSurfaceMesh::update()
         return false;
     }
 
+    // Calculate surface and mesh overlap bounding box
+    treeBoundBox bb
+    (
+        surface_.triSurface::points(),
+        surface_.triSurface::meshPoints()
+    );
+    bb.min() = max(bb.min(), mesh().bounds().min());
+    bb.max() = min(bb.max(), mesh().bounds().max());
+
+    // Extend a bit
+    const vector span(bb.span());
+
+    bb.min() -= 0.5*span;
+    bb.max() += 0.5*span;
+
+    bb.inflate(1e-6);
+
     // Mesh search engine, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), polyMesh::FACE_PLANES);
+    meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);
 
     return update(meshSearcher);
 }

src/triSurface/triSurface/interfaces/NAS/readNAS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ Description
 
     - Uses the Ansa "$ANSA_NAME" or the Hypermesh "$HMNAME COMP" extensions
       to obtain patch names.
-    - Handles Nastran short and long formats, but not free format.
+    - Handles Nastran short, long, and comma-separated free formats.
     - Properly handles the Nastran compact floating point notation: \n
     \verbatim
         GRID          28        10.20269-.030265-2.358-8
@@ -69,6 +69,32 @@ static scalar parseNASCoord(const string& s)
 }
 
 
+// Read a column of a given width from either a fixed-format NAS file, or a
+// comma-separated free-format NAS file
+static std::string readNASToken
+(
+    const string& line,
+    const size_t& width,
+    size_t& index
+)
+{
+    size_t indexStart, indexEnd;
+
+    indexStart = index;
+
+    indexEnd = line.find(',', indexStart);
+    index = indexEnd + 1;
+
+    if (indexEnd == std::string::npos)
+    {
+        indexEnd = indexStart + width;
+        index = indexEnd;
+    }
+
+    return line.substr(indexStart, indexEnd - indexStart);
+}
+
+
 bool triSurface::readNAS(const fileName& fName)
 {
     IFstream is(fName);
@@ -104,6 +130,7 @@ bool triSurface::readNAS(const fileName& fName)
 
     while (is.good())
     {
+        size_t linei = 0;
         string line;
         is.getLine(line);
 
@@ -197,17 +224,16 @@ bool triSurface::readNAS(const fileName& fName)
         }
 
         // Read first word
-        IStringStream lineStream(line);
+        word cmd(IStringStream(readNASToken(line, 8, linei))());
-        word cmd;
-        lineStream >> cmd;
 
         if (cmd == "CTRIA3")
         {
-            label groupId = readLabel(IStringStream(line.substr(16,8))());
+            readNASToken(line, 8, linei);
-            label a = readLabel(IStringStream(line.substr(24,8))());
+            label groupId =
-            label b = readLabel(IStringStream(line.substr(32,8))());
+                readLabel(IStringStream(readNASToken(line, 8, linei))());
-            label c = readLabel(IStringStream(line.substr(40,8))());
+            label a = readLabel(IStringStream(readNASToken(line, 8, linei))());
-
+            label b = readLabel(IStringStream(readNASToken(line, 8, linei))());
+            label c = readLabel(IStringStream(readNASToken(line, 8, linei))());
 
             // Convert group into patch
             Map<label>::const_iterator iter = groupToPatch.find(groupId);
@@ -228,11 +254,13 @@ bool triSurface::readNAS(const fileName& fName)
         }
         else if (cmd == "CQUAD4")
         {
-            label groupId = readLabel(IStringStream(line.substr(16,8))());
+            readNASToken(line, 8, linei);
-            label a = readLabel(IStringStream(line.substr(24,8))());
+            label groupId =
-            label b = readLabel(IStringStream(line.substr(32,8))());
+                readLabel(IStringStream(readNASToken(line, 8, linei))());
-            label c = readLabel(IStringStream(line.substr(40,8))());
+            label a = readLabel(IStringStream(readNASToken(line, 8, linei))());
-            label d = readLabel(IStringStream(line.substr(48,8))());
+            label b = readLabel(IStringStream(readNASToken(line, 8, linei))());
+            label c = readLabel(IStringStream(readNASToken(line, 8, linei))());
+            label d = readLabel(IStringStream(readNASToken(line, 8, linei))());
 
             // Convert group into patch
             Map<label>::const_iterator iter = groupToPatch.find(groupId);
@@ -255,7 +283,8 @@ bool triSurface::readNAS(const fileName& fName)
         else if (cmd == "PSHELL")
         {
             // Read shell type since group gives patchnames
-            label groupId = readLabel(IStringStream(line.substr(8,8))());
+            label groupId =
+                readLabel(IStringStream(readNASToken(line, 8, linei))());
             if (groupId == ansaId && ansaType == "PSHELL")
             {
                 groupToName.insert(groupId, string::validate<word>(ansaName));
@@ -264,10 +293,12 @@ bool triSurface::readNAS(const fileName& fName)
         }
         else if (cmd == "GRID")
         {
-            label index = readLabel(IStringStream(line.substr(8,8))());
+            label index =
-            scalar x = parseNASCoord(line.substr(24, 8));
+                readLabel(IStringStream(readNASToken(line, 8, linei))());
-            scalar y = parseNASCoord(line.substr(32, 8));
+            readNASToken(line, 8, linei);
-            scalar z = parseNASCoord(line.substr(40, 8));
+            scalar x = parseNASCoord(readNASToken(line, 8, linei));
+            scalar y = parseNASCoord(readNASToken(line, 8, linei));
+            scalar z = parseNASCoord(readNASToken(line, 8, linei));
 
             indices.append(index);
             points.append(point(x, y, z));
@@ -279,11 +310,13 @@ bool triSurface::readNAS(const fileName& fName)
             // Typical line (spaces compacted)
             // GRID*      126   0 -5.55999875E+02 -5.68730474E+02
             // *         2.14897901E+02
+            label index =
+                readLabel(IStringStream(readNASToken(line, 8, linei))());
+            readNASToken(line, 8, linei);
+            scalar x = parseNASCoord(readNASToken(line, 16, linei));
+            scalar y = parseNASCoord(readNASToken(line, 16, linei));
 
-            label index = readLabel(IStringStream(line.substr(8,16))());
+            linei = 0;
-            scalar x = parseNASCoord(line.substr(40, 16));
-            scalar y = parseNASCoord(line.substr(56, 16));
-
             is.getLine(line);
             if (line[0] != '*')
             {
@@ -295,7 +328,8 @@ bool triSurface::readNAS(const fileName& fName)
                     << " line:" << is.lineNumber()
                     << exit(FatalError);
             }
-            scalar z = parseNASCoord(line.substr(8, 16));
+            readNASToken(line, 8, linei);
+            scalar z = parseNASCoord(readNASToken(line, 16, linei));
 
             indices.append(index);
             points.append(point(x, y, z));

src/triSurface/triSurface/interfaces/OBJ/readOBJ.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,7 +53,9 @@ bool Foam::triSurface::readOBJ(const fileName& OBJfileName)
 
         label sz = line.size();
 
-        if (sz && line[sz-1] == '\\')
+        if (sz)
+        {
+            if (line[sz-1] == '\\')
             {
                 line.substr(0, sz-1);
                 line += getLineNoComment(OBJfile);
@@ -120,7 +122,11 @@ bool Foam::triSurface::readOBJ(const fileName& OBJfileName)
                     }
                     else
                     {
-                    vertexSpec = line.substr(startNum, line.size() - startNum);
+                        vertexSpec = line.substr
+                        (
+                            startNum,
+                            line.size() - startNum
+                        );
                     }
 
                     string::size_type slashPos = vertexSpec.find('/');
@@ -155,6 +161,7 @@ bool Foam::triSurface::readOBJ(const fileName& OBJfileName)
                 }
             }
         }
+    }
 
     points.shrink();
     faces.shrink();

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/T

@@ -1,56 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      T;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 1 0 0 0];
-
-internalField   uniform 293;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            zeroGradient;
-    }
-
-    topAndBottom
-    {
-        type            zeroGradient;
-    }
-
-    hot
-    {
-        type            externalCoupledTemperature;
-        commsDir        "${FOAM_CASE}/comms";
-        fileName        "data";
-        initByExternal  yes;
-        log             true;
-        value           uniform 307.75; // 34.6 degC
-    }
-
-    cold
-    {
-        type            externalCoupledTemperature;
-        commsDir        "${FOAM_CASE}/comms";
-        fileName        "data";
-        initByExternal  yes;
-        log             true;
-        value           uniform 288.15; // 15 degC
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/alphat

@@ -1,51 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      alphat;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    topAndBottom
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    hot
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    cold
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/epsilon

@@ -1,47 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      epsilon;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -3 0 0 0 0];
-
-internalField   uniform 4e-06;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
-    }
-    topAndBottom
-    {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
-    }
-    hot
-    {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
-    }
-    cold
-    {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/omega

@@ -1,47 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      omega;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 -1 0 0 0 0];
-
-internalField   uniform 0.12;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-    topAndBottom
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-    hot
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-    cold
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allclean

@@ -1,14 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/CleanFunctions
-
-cleanCase
-
-rm -rf comms
-
-killall externalSolver > /dev/null 2>&1
-
-
-# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun

@@ -1,14 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-./Allrun.pre
-
-runApplication $(getApplication) &
-
-./externalSolver
-
-
-# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun-parallel

@@ -1,16 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-./Allrun.pre
-
-runApplication decomposePar
-
-runParallel $(getApplication) 4 &
-
-./externalSolver
-
-
-# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun.pre

@@ -1,11 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-runApplication blockMesh
-runApplication createExternalCoupledPatchGeometry T
-
-
-# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/README

@@ -1,5 +0,0 @@
-Example of an explicit coupling between OpenFOAM and an external application
-using the externalCoupled boundary conditions.
-
-The case is based on the buoyantCavity tutorial case, whereby on each iteration
-the 'hot' and 'cold' patch temperatures are incremented by 1K.

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/T

@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+        value           uniform 300;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/U

@@ -10,41 +10,21 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0 0 0);
+internalField   uniform (0.01 0 0);
 
 boundaryField
 {
-    frontAndBack
+    ".*"
     {
-        type            fixedValue;
+        type            calculated;
-        value           uniform (0 0 0);
+        value           uniform (0.01 0 0);
-    }
-
-    topAndBottom
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    hot
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    cold
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
     }
 }
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/epsilon

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/k

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0.1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,30 +20,11 @@ internalField   uniform 1e5;
 
 boundaryField
 {
-    frontAndBack
+    ".*"
     {
         type            calculated;
-        value           $internalField;
+        value           uniform 1e5;
-    }
-
-    topAndBottom
-    {
-        type            calculated;
-        value           $internalField;
-    }
-
-    hot
-    {
-        type            calculated;
-        value           $internalField;
-    }
-
-    cold
-    {
-        type            calculated;
-        value           $internalField;
     }
 }
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p_rgh

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,30 +20,11 @@ internalField   uniform 1e5;
 
 boundaryField
 {
-    frontAndBack
+    ".*"
     {
-        type            fixedFluxPressure;
+        type            calculated;
-        value           uniform 1e5;
-    }
-
-    topAndBottom
-    {
-        type            fixedFluxPressure;
-        value           uniform 1e5;
-    }
-
-    hot
-    {
-        type            fixedFluxPressure;
-        value           uniform 1e5;
-    }
-
-    cold
-    {
-        type            fixedFluxPressure;
         value           uniform 1e5;
     }
 }
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allclean

@@ -0,0 +1,23 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+rm -rf comms
+rm -rf VTK
+rm -rf constant/cellToRegion constant/polyMesh/sets
+rm -rf 0/bottomWater
+rm -rf 0/topAir
+rm -rf 0/heater
+rm -rf 0/leftSolid
+rm -rf 0/rightSolid
+rm -f 0/cellToRegion
+rm -rf constant/bottomWater/polyMesh
+rm -rf constant/topAir/polyMesh
+rm -rf constant/heater/polyMesh
+rm -rf constant/leftSolid/polyMesh
+rm -rf constant/rightSolid/polyMesh
+
+# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun

@@ -0,0 +1,28 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+./Allrun.pre
+
+#-- Run on single processor
+#runApplication `getApplication` &
+# Simulated external solver
+#runApplication ./externalSolver
+
+# Decompose
+runApplication decomposePar -allRegions
+
+# Run OpenFOAM
+runParallel `getApplication` 4 &
+
+# Simulated external solver
+runApplication ./externalSolver
+
+# Reconstruct
+runApplication reconstructPar -allRegions
+
+
+# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun.pre

@@ -0,0 +1,33 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication topoSet
+runApplication splitMeshRegions -cellZones -overwrite
+
+# remove fluid fields from solid regions (important for post-processing)
+for i in heater leftSolid rightSolid
+do
+   rm -f 0*/$i/{nut,alphat,epsilon,k,U,p_rgh}
+done
+
+
+for i in bottomWater topAir heater leftSolid rightSolid
+do
+   changeDictionary -region $i > log.changeDictionary.$i 2>&1
+done
+
+# Create coupling geometry
+runApplication createExternalCoupledPatchGeometry \
+    -regions '(topAir heater)' coupleGroup
+
+echo
+echo "creating files for paraview post-processing"
+echo
+paraFoam -touchAll
+
+# ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/README.txt

@@ -0,0 +1,3 @@
+Modification of the heatTransfer chtMultiRegionFoam tutorial that demonstrates
+the externalCoupled functionObject in combination with the ./externalSolver
+script to simulate coupling to an external code.

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/g

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       uniformDimensionedVectorField;
-    location    "constant";
     object      g;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -18,5 +17,4 @@ FoamFile
 dimensions      [0 1 -2 0 0 0 0];
 value           (0 -9.81 0);
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/radiationProperties

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/thermophysicalProperties

@@ -9,39 +9,43 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
+    class       dictionary;
-    location    "0";
+    object      thermophysicalProperties;
-    object      k;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -2 0 0 0 0];
+thermoType
-
-internalField   uniform 3.75e-04;
-
-boundaryField
 {
-    frontAndBack
+    type            heRhoThermo;
-    {
+    mixture         pureMixture;
-        type            kqRWallFunction;
+    transport       const;
-        value           uniform 3.75e-04;
+    thermo          hConst;
-    }
+    equationOfState rhoConst;
-    topAndBottom
+    specie          specie;
-    {
+    energy          sensibleEnthalpy;
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
-    hot
-    {
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
-    cold
-    {
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
 }
 
+mixture
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       18;
+    }
+    equationOfState
+    {
+        rho             1000;
+    }
+    thermodynamics
+    {
+        Cp              4181;
+        Hf              0;
+    }
+    transport
+    {
+        mu              959e-6;
+        Pr              6.62;
+    }
+}
 
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/turbulenceProperties

@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/radiationProperties

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation off;
+
+radiationModel  none;
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/thermophysicalProperties

@@ -9,39 +9,45 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
+    class       dictionary;
-    location    "0";
+    object      thermophysicalProperties;
-    object      nut;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -1 0 0 0 0];
+thermoType
-
-internalField   uniform 0;
-
-boundaryField
 {
-    frontAndBack
+    type            heSolidThermo;
-    {
+    mixture         pureMixture;
-        type            nutUWallFunction;
+    transport       constIso;
-        value           uniform 0;
+    thermo          hConst;
-    }
+    equationOfState rhoConst;
-    topAndBottom
+    specie          specie;
-    {
+    energy          sensibleEnthalpy;
-        type            nutUWallFunction;
-        value           uniform 0;
-    }
-    hot
-    {
-        type            nutUWallFunction;
-        value           uniform 0;
-    }
-    cold
-    {
-        type            nutUWallFunction;
-        value           uniform 0;
-    }
 }
 
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   50;
+    }
+
+    transport
+    {
+        kappa   80;
+    }
+
+    thermodynamics
+    {
+        Hf      0;
+        Cp      450;
+    }
+
+    equationOfState
+    {
+        rho     8000;
+    }
+}
 
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/radiationProperties

@@ -0,0 +1 @@
+../heater/radiationProperties

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/thermophysicalProperties

@@ -0,0 +1 @@
+../heater/thermophysicalProperties

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/blockMeshDict

@@ -14,66 +14,77 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-convertToMeters 0.001;
+convertToMeters 1;
 
 vertices
 (
-    ( 0     0  -260)
+    (-0.1 -0.04  -0.05)
-    (76     0  -260)
+    ( 0.1 -0.04  -0.05)
-    (76  2180  -260)
+    ( 0.1  0.04  -0.05)
-    ( 0  2180  -260)
+    (-0.1  0.04  -0.05)
-    ( 0     0   260)
+    (-0.1 -0.04   0.05)
-    (76     0   260)
+    ( 0.1 -0.04   0.05)
-    (76  2180   260)
+    ( 0.1  0.04   0.05)
-    ( 0  2180   260)
+    (-0.1  0.04   0.05)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (30 10 10) simpleGrading (1 1 1)
 );
 
 edges
 (
 );
 
-blocks
-(
-    hex (0 1 2 3 4 5 6 7) (35 150 15) simpleGrading (1 1 1)
-);
-
 boundary
 (
-    frontAndBack
+    maxY
     {
         type wall;
         faces
         (
-            (0 1 5 4)
+            (3 7 6 2)
-            (2 3 7 6)
         );
     }
-
+    minX
-    topAndBottom
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    minY
+    {
+        type wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    minZ
+    {
+        type wall;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    maxZ
     {
         type wall;
         faces
         (
             (4 5 6 7)
-            (3 2 1 0)
-        );
-    }
-
-    hot
-    {
-        type wall;
-        faces
-        (
-            (6 5 1 2)
-        );
-    }
-
-    cold
-    {
-        type wall;
-        faces
-        (
-            (4 7 3 0)
         );
     }
 );
@@ -81,3 +92,5 @@ boundary
 mergePatchPairs
 (
 );
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
+|  \\    /   O peration     | Version:  dev-OpenCFD.feature-externalCoupled   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,35 +15,47 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-4
+6
 (
-    frontAndBack
+    maxY
     {
         type            wall;
         inGroups        1(wall);
-        nFaces          1050;
+        nFaces          300;
-        startFace       228225;
+        startFace       8300;
     }
-    topAndBottom
+    minX
     {
-        type            wall;
+        type            patch;
-        inGroups        1(wall);
+        nFaces          100;
-        nFaces          10500;
+        startFace       8600;
-        startFace       229275;
     }
-    hot
+    maxX
     {
-        type            wall;
+        type            patch;
-        inGroups        1(wall);
+        nFaces          100;
-        nFaces          2250;
+        startFace       8700;
-        startFace       239775;
     }
-    cold
+    minY
     {
         type            wall;
         inGroups        1(wall);
-        nFaces          2250;
+        nFaces          300;
-        startFace       242025;
+        startFace       8800;
+    }
+    minZ
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          300;
+        startFace       9100;
+    }
+    maxZ
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          300;
+        startFace       9400;
     }
 )
 

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/regionProperties

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      regionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+regions
+(
+    fluid       (bottomWater topAir)
+    solid       (heater leftSolid rightSolid)
+);
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/radiationProperties

@@ -0,0 +1 @@
+../heater/radiationProperties

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/thermophysicalProperties

@@ -0,0 +1 @@
+../heater/thermophysicalProperties

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/g

@@ -0,0 +1 @@
+../bottomWater/g

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/radiationProperties

@@ -0,0 +1 @@
+../bottomWater/radiationProperties

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/thermophysicalProperties

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "constant/bottomWater";
     object      thermophysicalProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -31,17 +31,17 @@ mixture
     specie
     {
         nMoles          1;
-        molWeight       28.96;
+        molWeight       28.9;
     }
     thermodynamics
     {
-        Cp              1004.4;
+        Cp              1000;
         Hf              0;
     }
     transport
     {
-        mu              1.831e-05;
+        mu              1.8e-05;
-        Pr              0.705;
+        Pr              0.7;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/turbulenceProperties

@@ -0,0 +1 @@
+../bottomWater/turbulenceProperties

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/externalSolver

@@ -1,20 +1,29 @@
 #!/bin/sh
 #
 # Dummy external solver to communicate with OpenFOAM via externalCoupled
-# boundary conditions
+# functionObject
 #
 # Functionality is hard-coded for this particular test case
 # - patch temperatures increased by 1K on each step
 #
-cd ${0%/*} || exit 1    # run from this directory
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Check for unassigned variables
+set -u
 
 echo "Executing dummy external solver"
 
 commsDir="comms"
+regionGroupName="heater_topAir"
+patchGroupName="coupleGroup"
+fieldName="T"
+
 lockFile="${commsDir}/OpenFOAM.lock"
-dataFile="${commsDir}/data"
+dataFile="${commsDir}/${regionGroupName}/${patchGroupName}/${fieldName}"
 waitSec=1
 timeOut=10
+nSteps=200  # maximum number of time steps. Note: should be more than
+            # number of iterations on the OpenFOAM side
 refGrad=0
 valueFraction=1
 
@@ -27,16 +36,21 @@ init()
 {
     log "initialisation: creating ${dataFile}.in"
 
-    # Hard-coded for 2 patches of size 2250
+    # Hard-coded for patch of size 8 (heater/minY)
-    n=2250
+    n1=8
-    refCold=283
+    refValue1=500
-    refHot=303
     touch "${dataFile}.in"
-    for i in $(seq 1 $n); do
+    log "initialisation: adding $n1 data elements with refValue $refValue1"
-        echo "$refHot $refGrad $valueFraction" >> "${dataFile}.in"
+    for i in $(seq 1 $n1); do
+        echo "$refValue1 $refGrad $valueFraction" >> "${dataFile}.in"
     done
-    for i in $(seq 1 $n); do
+
-        echo "$refCold $refGrad $valueFraction" >> "${dataFile}.in"
+    # Hard-coded for patch of size 40 (topAir/minX)
+    n2=40
+    refValue2=300
+    log "initialisation: adding $n2 data elements with refValue $refValue2"
+    for i in $(seq 1 $n2); do
+        echo "$refValue2 $refGrad $valueFraction" >> "${dataFile}.in"
     done
 
     # create lock file to pass control to OF
@@ -44,6 +58,10 @@ init()
 }
 
 
+# create the comms directory
+mkdir -p ${commsDir}/${regionGroupName}/${patchGroupName}
+
+
 # tutorial case employs the 'initByExternalOption', so we need to provide
 # the initial values
 init
@@ -51,7 +69,7 @@ init
 
 totalWait=0
 step=0
-while [ 1 ]; do
+while [ $step -lt $nSteps ]; do
     if [ -f $lockFile ]; then
         log "found lock file ${lockFile} - waiting"
         totalWait=$(expr $totalWait + $waitSec)
@@ -70,9 +88,10 @@ while [ 1 ]; do
         log "sleeping for $waitSec secs to simulate external process"
         sleep $waitSec
 
-        log "creating ${dataFile}.in"
+        log "updating ${dataFile}.in from ${dataFile}.out"
 
-        awk '{if( $1 != "#" ){print $2+1 " 0 1"}}' ${dataFile}.out > ${dataFile}.in
+        awk '{if( $1 != "#" ){print $1+1 " 0 1"}}' \
+            ${dataFile}.out | tee ${dataFile}.in
 
         log "creating lock file ${lockFile}"
         touch ${lockFile}

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/README

@@ -0,0 +1,3 @@
+fvSolution is used for outer correctors specification.
+fvSchemes is only so that pre-processing activities can proceed
+

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/changeDictionaryDict

@@ -0,0 +1,173 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    U
+    {
+        internalField   uniform (0.001 0 0);
+
+        boundaryField
+        {
+            minX
+            {
+                type            fixedValue;
+                value           uniform (0.001 0 0);
+            }
+
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform (0 0 0);
+            }
+
+            ".*"
+            {
+                type            fixedValue;
+                value           uniform (0 0 0);
+            }
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            minX
+            {
+                type            fixedValue;
+                value           uniform 300;
+            }
+
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 300;
+            }
+
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+
+            "bottomWater_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           fluidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+        }
+    }
+
+    epsilon
+    {
+        internalField   uniform 0.01;
+
+        boundaryField
+        {
+            minX
+            {
+                type            fixedValue;
+                value           uniform 0.01;
+            }
+
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.01;
+            }
+
+            ".*"
+            {
+                type            epsilonWallFunction;
+                value           uniform 0.01;
+            }
+        }
+    }
+
+    k
+    {
+        internalField   uniform 0.1;
+
+        boundaryField
+        {
+            minX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.1;
+            }
+
+            maxX
+            {
+                type            zeroGradient;
+                value           uniform 0.1;
+            }
+
+            ".*"
+            {
+                type            kqRWallFunction;
+                value           uniform 0.1;
+            }
+        }
+    }
+
+    p_rgh
+    {
+        internalField   uniform 0;
+
+        boundaryField
+        {
+            minX
+            {
+                type            zeroGradient;
+                value           uniform 0;
+            }
+
+            maxX
+            {
+                type            fixedValue;
+                value           uniform 0;
+            }
+
+            ".*"
+            {
+                type            fixedFluxPressure;
+                value           uniform 0;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 0;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 0;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/decomposeParDict

@@ -0,0 +1,72 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    location    "system";
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+method          scotch;
+// method          hierarchical;
+// method          simple;
+// method          manual;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+
+//// Is the case distributed
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSchemes

@@ -16,7 +16,7 @@ FoamFile
 
 ddtSchemes
 {
-    default steadyState;
+    default Euler;
 }
 
 gradSchemes
@@ -28,18 +28,19 @@ divSchemes
 {
     default         none;
 
-    div(phi,U)      bounded Gauss limitedLinear 0.2;
+    div(phi,U)      Gauss upwind;
-    div(phi,K)      bounded Gauss limitedLinear 0.2;
+    div(phi,K)      Gauss linear;
-    div(phi,h)      bounded Gauss limitedLinear 0.2;
+    div(phi,h)      Gauss upwind;
-    div(phi,k)      bounded Gauss limitedLinear 0.2;
+    div(phi,k)      Gauss upwind;
-    div(phi,epsilon) bounded Gauss limitedLinear 0.2;
+    div(phi,epsilon) Gauss upwind;
-    div(phi,omega) bounded Gauss limitedLinear 0.2;
+    div(phi,R)      Gauss upwind;
+    div(R)          Gauss linear;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         Gauss linear orthogonal;
+    default        Gauss linear corrected;
 }
 
 interpolationSchemes
@@ -49,13 +50,13 @@ interpolationSchemes
 
 snGradSchemes
 {
-    default         orthogonal;
+    default         corrected;
 }
 
-wallDist
+fluxRequired
 {
-    method meshWave;
+    default         no;
+    p_rgh;
 }
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSolution

@@ -10,20 +10,34 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 solvers
 {
+    rho
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-7;
+        relTol          0;
+    }
+
     p_rgh
     {
         solver           GAMG;
         tolerance        1e-7;
         relTol           0.01;
 
-        smoother         DICGaussSeidel;
+        smoother         GaussSeidel;
 
         cacheAgglomeration  true;
         nCellsInCoarsestLevel 10;
@@ -31,47 +45,43 @@ solvers
         mergeLevels      1;
     }
 
-    "(U|h|k|epsilon|omega)"
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance        1e-7;
+        relTol           0;
+    }
+
+    "(U|h|k|epsilon|R)"
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance       1e-8;
+        tolerance        1e-7;
         relTol           0.1;
     }
+
+    "(U|h|k|epsilon|R)Final"
+    {
+        $U;
+        tolerance        1e-7;
+        relTol           0;
+    }
 }
 
-SIMPLE
+PIMPLE
 {
-    momentumPredictor yes;
+    momentumPredictor   on;
+    nCorrectors         2;
     nNonOrthogonalCorrectors 0;
-    pRefCell        0;
-    pRefValue       0;
-
-    residualControl
-    {
-        p_rgh           1e-2;
-        U               1e-3;
-        h               1e-3;
-
-        // possibly check turbulence fields
-        "(k|epsilon|omega)" 1e-3;
-    }
 }
 
 relaxationFactors
 {
-    fields
-    {
-        rho             1.0;
-        p_rgh           0.7;
-    }
     equations
     {
-        U               0.3;
+        "h.*"           1;
-        h               0.3;
+        "U.*"           1;
-        "(k|epsilon|omega)" 0.7;
     }
 }
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/controlDict

@@ -0,0 +1,96 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Library defines new boundary conditions
+libs            ("libOpenFOAM.so" "libjobControl.so");
+
+application     chtMultiRegionFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+writeInterval   0.1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  8;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+maxCo           0.6;
+
+// Maximum diffusion number
+maxDi           10.0;
+
+adjustTimeStep  yes;
+
+functions
+{
+    externalCoupled
+    {
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libjobControl.so");
+
+        type            externalCoupled;
+
+        // Directory to use for communication
+        commsDir        "${FOAM_CASE}/comms";
+
+        // Does external process start first
+        initByExternal  true;
+
+        // Additional output
+        log             true;
+
+        regions
+        {
+            // Region name (wildcards allowed)
+            "(topAir|heater)"
+            {
+                // In topAir adjust the minX patch (fixedValue)
+
+                // Patch or patchGroup
+                coupleGroup
+                {
+                    // Fields to output in commsDir
+                    writeFields (T);
+                    // Fields to read from commsDir
+                    readFields  (T);
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/decomposeParDict

@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    location    "system";
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+method          scotch;
+// method          hierarchical;
+// method          simple;
+// method          manual;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+
+//// Is the case distributed
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSchemes

@@ -10,19 +10,36 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
+    object      fvSchemes;
-    object      decomposeParDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains 4;
+ddtSchemes
-
+{
-method          simple;
+}
-
+
-simpleCoeffs
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+fluxRequired
 {
-    n               (2 2 1);
-    delta           0.001;
 }
 
 

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSolution

@@ -10,20 +10,13 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      RASProperties;
+    object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType RAS;
+PIMPLE
-
-RAS
 {
-    RASModel            kOmegaSST;
+    nOuterCorrectors 1;
-
-    turbulence          on;
-
-    printCoeffs         on;
 }
 
-
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/changeDictionaryDict

@@ -0,0 +1,77 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    boundary
+    {
+        minY
+        {
+            type            patch;
+            inGroups        (coupleGroup);
+        }
+        minZ
+        {
+            type            patch;
+        }
+        maxZ
+        {
+            type            patch;
+        }
+    }
+
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+            "heater_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+
+            heater_to_leftSolid
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                thicknessLayers (1e-3);
+                kappaLayers     (5e-4);
+                value           uniform 300;
+            }
+
+            minY
+            {
+                type            fixedValue;
+                value           uniform 500;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/decomposeParDict

@@ -10,39 +10,35 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      controlDict;
+    location    "system";
+    object      decomposeParDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     buoyantSimpleFoam;
+numberOfSubdomains  4;
 
-startFrom       startTime;
+method          scotch;
 
-startTime       0;
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
 
-stopAt          endTime;
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
 
-endTime         100;
+scotchCoeffs
-
+{
-deltaT          1;
+}
-
-writeControl    timeStep;
-
-writeInterval   10;
-
-purgeWrite      0;
-
-writeFormat     ascii;
-
-writePrecision  6;
-
-writeCompression off;
-
-timeFormat      general;
-
-timePrecision   6;
-
-runTimeModifiable true;
 
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
 
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSchemes

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default             none;
+    laplacian(alpha,h)  Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSolution

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    h
+    {
+        solver           PCG;
+        preconditioner   DIC;
+        tolerance        1e-06;
+        relTol           0.1;
+    }
+
+    hFinal
+    {
+        $h;
+        tolerance        1e-06;
+        relTol           0;
+    }
+}
+
+PIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+// ************************************************************************* //