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/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));
     }
 
 );
@@ -249,8 +307,10 @@ surfaces
         isoValue        0.5;
         interpolate     true;
 
-        //zone            ABC;          // Optional: zone only
+        // zone            ABC;         // Optional: zone only
-        //exposedPatchName fixedWalls;  // 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
@@ -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)
     }
@@ -281,8 +344,10 @@ surfaces
         }
         interpolate     true;
 
-        //zone            ABC;          // Optional: zone only
+        // zone             ABC;          // Optional: zone only
-        //exposedPatchName fixedWalls;  // 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
@@ -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;
 

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/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/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&>
                 (
@@ -136,12 +159,13 @@ bool Foam::externalCoupledFunctionObject::readData
                 readColumns
                 (
                     bf[patchI].size(),              // number of lines to read
-                4*pTraits<Type>::nComponents+1, // nColumns: 4*Type + 1*scalar
+                    4*pTraits<Type>::nComponents+1, // nColumns: 4*Type+1*scalar
                     masterFilePtr,
                     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/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/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,52 +122,76 @@ 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
         {
             // Not in registry. Interpolate.
 
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::getIsoFields() : "
-                    << "checking pointField " << pointFldName
-                    << " for same time " << fvm.time().timeName()
-                    << endl;
-            }
-
-            if
-            (
-                storedPointFieldPtr_.empty()
-             || (fvm.time().timeName() != storedPointFieldPtr_().instance())
-            )
-            {
             if (debug)
             {
                 Info<< "sampledIsoSurface::getIsoFields() :"
-                        << " interpolating volField " << volFieldPtr_->name()
+                    << " creating and storing pointField "
-                        << " to get pointField " << pointFldName << endl;
+                    << pointFldName << " for time "
+                    << fvm.time().timeName() << endl;
             }
 
-                storedPointFieldPtr_.reset
+            tmp<pointScalarField> tpfld
             (
-                    volPointInterpolation::New(fvm)
+                volPointInterpolation::New(fvm).interpolate
-                    .interpolate(*volFieldPtr_).ptr()
+                (
+                    *volFieldPtr_,
+                    pointFldName,
+                    false
+                )
             );
-                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);
 }

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;
+    transport       const;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
     {
-        type            kqRWallFunction;
+        nMoles          1;
-        value           uniform 3.75e-04;
+        molWeight       18;
     }
-    topAndBottom
+    equationOfState
     {
-        type            kqRWallFunction;
+        rho             1000;
-        value           uniform 3.75e-04;
     }
-    hot
+    thermodynamics
     {
-        type            kqRWallFunction;
+        Cp              4181;
-        value           uniform 3.75e-04;
+        Hf              0;
     }
-    cold
+    transport
     {
-        type            kqRWallFunction;
+        mu              959e-6;
-        value           uniform 3.75e-04;
+        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;
+    transport       constIso;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
     {
-        type            nutUWallFunction;
+        nMoles      1;
-        value           uniform 0;
+        molWeight   50;
     }
-    topAndBottom
+
+    transport
     {
-        type            nutUWallFunction;
+        kappa   80;
-        value           uniform 0;
     }
-    hot
+
+    thermodynamics
     {
-        type            nutUWallFunction;
+        Hf      0;
-        value           uniform 0;
+        Cp      450;
     }
-    cold
+
+    equationOfState
     {
-        type            nutUWallFunction;
+        rho     8000;
-        value           uniform 0;
     }
 }
 
-
 // ************************************************************************* //

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;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/changeDictionaryDict

@@ -0,0 +1,65 @@
+/*--------------------------------*- 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
+    {
+        minZ
+        {
+            type            patch;
+        }
+        maxZ
+        {
+            type            patch;
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+            "leftSolid_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+
+            leftSolid_to_heater
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                thicknessLayers (1e-3);
+                kappaLayers     (5e-4);
+                value           uniform 300;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/decomposeParDict

@@ -0,0 +1,44 @@
+/*--------------------------------*- 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      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method          scotch;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSchemes

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

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSolution

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

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/changeDictionaryDict

@@ -0,0 +1,54 @@
+/*--------------------------------*- 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
+    {
+        minZ
+        {
+            type            patch;
+        }
+        maxZ
+        {
+            type            patch;
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+            "rightSolid_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/decomposeParDict

@@ -0,0 +1,44 @@
+/*--------------------------------*- 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      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method          scotch;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSchemes

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

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSolution

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

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/changeDictionaryDict

@@ -0,0 +1,179 @@
+/*--------------------------------*- 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
+    {
+        minX
+        {
+            inGroups        (coupleGroup);
+        }
+    }
+
+    U
+    {
+        internalField   uniform (0.1 0 0);
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            fixedValue;
+                value           uniform (0 0 0);
+            }
+            minX
+            {
+                type            fixedValue;
+                value           uniform ( 0.1 0 0 );
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform ( 0 0 0 );
+                value           uniform ( 0.1 0 0 );
+            }
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+            }
+
+            minX
+            {
+                type            fixedValue;
+                value           uniform 300;
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 300;
+                value           uniform 300;
+            }
+
+            "topAir_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           fluidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+        }
+    }
+
+    epsilon
+    {
+        internalField   uniform 0.01;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            epsilonWallFunction;
+                value           uniform 0.01;
+            }
+
+            minX
+            {
+                type            fixedValue;
+                value           uniform 0.01;
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.01;
+                value           uniform 0.01;
+            }
+        }
+    }
+
+    k
+    {
+        internalField   uniform 0.1;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            kqRWallFunction;
+                value           uniform 0.1;
+            }
+
+            minX
+            {
+                type            fixedValue;
+                value           uniform 0.1;
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.1;
+                value           uniform 0.1;
+            }
+        }
+    }
+
+    p_rgh
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            fixedFluxPressure;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            fixedValue;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/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/topAir/fvSchemes

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

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSolution

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

tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topoSetDict

@@ -0,0 +1,178 @@
+/*--------------------------------*- 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      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // Heater
+    {
+        name    heaterCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-0.01001    0 -100 )(0.01001 0.00999 100);
+        }
+    }
+    {
+        name    heaterCellSet;
+        type    cellSet;
+        action  add;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-0.01001 -100 -0.01001)(0.01001 0.00999 0.01001);
+        }
+    }
+    {
+        name    heater;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set heaterCellSet;
+        }
+    }
+
+    // leftSolid
+    {
+        name    leftSolidCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-100 0 -100 )(-0.01001 0.00999 100);
+        }
+    }
+    {
+        name    leftSolid;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set leftSolidCellSet;
+        }
+    }
+
+    // rightSolid
+    {
+        name    rightSolidCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (0.01001 0 -100 )(100 0.00999 100);
+        }
+    }
+    {
+        name    rightSolid;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set rightSolidCellSet;
+        }
+    }
+
+    // topAir
+    {
+        name    topAirCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-100 0.00999 -100 )(100 100 100);
+        }
+    }
+    {
+        name    topAir;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set topAirCellSet;
+        }
+    }
+
+
+    // bottomWater is all the other cells
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        sourceInfo
+        {
+            set heaterCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set leftSolidCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set rightSolidCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set topAirCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  invert;
+    }
+    {
+        name    bottomWater;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set bottomWaterCellSet;
+        }
+    }
+);
+
+
+// ************************************************************************* //