Merge commit 'OpenCFD/master' into olesenm

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H

@@ -29,6 +29,36 @@ Description
     Basic sub-grid obstacle flame-wrinking enhancement factor model.
     Details supplied by J Puttock 2/7/06.
 
+    <b> Sub-grid flame area generation <\b>
+
+    \f$ n = N - \hat{\dwea{\vec{U}}}.n_{s}.\hat{\dwea{\vec{U}}} \f$
+    \f$ n_{r} = \sqrt{n} \f$
+
+    where:
+
+        \f$ \hat{\dwea{\vec{U}}} = \dwea{\vec{U}} / \vert \dwea{\vec{U}}
+        \vert \f$
+
+        \f$ b = \hat{\dwea{\vec{U}}}.B.\hat{\dwea{\vec{U}}} / n_{r} \f$
+
+    where:
+
+        \f$ B \f$ is the file "B".
+
+        \f$ N \f$ is the file "N".
+
+        \f$  n_{s} \f$ is the file "ns".
+
+    The flame area enhancement factor \f$ \Xi_{sub} \f$ is expected to
+    approach:
+
+    \f[
+        \Xi_{{sub}_{eq}} =
+            1 + max(2.2 \sqrt{b}, min(0.34 \frac{\vert \dwea{\vec{U}}
+            \vert}{{\vec{U}}^{'}}, 1.6)) \times min(\frac{n}{4}, 1)
+    \f]
+
+
 SourceFiles
     basicSubGrid.C
 

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H

@@ -25,10 +25,28 @@ License
 Class
     basicSubGrid
 
+
 Description
+
     Basic sub-grid obstacle flame-wrinking generation rate coefficient model.
     Details supplied by J Puttock 2/7/06.
 
+    \f$ G_{sub} \f$ denotes the generation coefficient and it is given by
+
+    \f[
+        G_{sub} = k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{obs}}
+                 \frac{/Xi_{{sub}_{eq}}-1}{/Xi_{sub}}
+    \f]
+
+    and the removal:
+
+    \f[ - k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{sub}}
+    \frac{\Xi_{sub}-1}{\Xi_{sub}} \f]
+
+    Finally, \f$ G_{sub} \f$ is added to generation rate \f$ G_{in} \f$
+    due to the turbulence.
+
+
 SourceFiles
     basicSubGrid.C
 

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.H

@@ -29,6 +29,50 @@ Description
     Basic sub-grid obstacle drag model.
     Details supplied by J Puttock 2/7/06.
 
+    <b> Sub-grid drag term <\b>
+
+    The resistance term (force per unit of volume) is given by:
+
+    \f[
+        R = -\frac{1}{2} \rho \vert \dwea{\vec{U}} \vert \dwea{\vec{U}}.D
+    \f]
+
+    where:
+
+        \f$ D \f$ is the tensor field "CR" in \f$ m^{-1} \f$
+
+        This is term is treated implicitly in UEqn.H
+
+    <b> Sub-grid turbulence generation <\b>
+
+    The turbulence source term \f$ G_{R} \f$ occurring in the
+    \f$ \kappa-\epsilon \f$ equations for the generation of turbulence due
+    to interaction with unresolved obstacles :
+
+    \f$ G_{R} = C_{s}\beta_{\nu}
+    \mu_{eff} A_{w}^{2}(\dwea{\vec{U}}-\dwea{\vec{U}_{s}})^2 + \frac{1}{2}
+    \rho \vert \dwea{\vec{U}} \vert \dwea{\vec{U}}.T.\dwea{\vec{U}} \f$
+
+    where:
+
+        \f$ C_{s} \f$ = 1
+
+        \f$ \beta_{\nu} \f$ is the volume porosity (file "betav").
+
+        \f$ \mu_{eff} \f$ is the effective viscosity.
+
+        \f$ A_{w}^{2}\f$ is the obstacle surface area per unit of volume
+        (file "Aw").
+
+        \f$ \dwea{\vec{U}_{s}} \f$ is the slip velocity and is considered
+        \f$ \frac{1}{2}. \dwea{\vec{U}} \f$.
+
+        \f$ T \f$ is a tensor in the file CT.
+
+    The term \f$ G_{R} \f$ is treated explicitly in the \f$ \kappa-\epsilon
+    \f$ Eqs in the PDRkEpsilon.C file.
+
+
 SourceFiles
     basic.C
 
@@ -40,7 +84,6 @@ SourceFiles
 #include "PDRDragModel.H"
 #include "XiEqModel.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H

@@ -26,7 +26,17 @@ Class
     PDRkEpsilon
 
 Description
-    Standard k-epsilon turbulence model.
+    Standard k-epsilon turbulence model with additional source terms
+    corresponding to PDR basic drag model (basic.H)
+
+    The turbulence source term \f$ G_{R} \f$ appears in the
+    \f$ \kappa-\epsilon \f$ equation for the generation of turbulence due to
+    interaction with unresolved obstacles.
+
+    In the \f$ \epsilon  \f$ equation \f$ C_{1} G_{R} \f$ is added as a source
+    term.
+
+    In the \f$ \kappa \f$ equation \f$ G_{R} \f$ is added as a source term.
 
 SourceFiles
     PDRkEpsilon.C

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H

@@ -27,6 +27,57 @@ Class
 
 Description
     Base-class for all Xi models used by the b-Xi combustion model.
+    See Technical Report SH/RE/01R for details on the PDR modelling.
+
+    Xi is given through an algebraic expression (algebraic.H),
+    by solving a transport equation (transport.H) or a fixed value (fixed.H).
+    See report TR/HGW/10 for details on the Weller two equations model.
+
+    In the algebraic and transport methods \f$\Xi_{eq}\f$ is calculated in
+    similar way. In the algebraic approach, \f$\Xi_{eq}\f$ is the value used in
+    the \f$ b \f$ transport equation.
+
+    \f$\Xi_{eq}\f$ is calculated as follows:
+
+    \f$\Xi_{eq} = 1 + (1 + 2\Xi_{coeff}(0.5 - \dwea{b}))(\Xi^* - 1)\f$
+
+    where:
+
+        \f$ \dwea{b} \f$ is the regress variable.
+
+        \f$ \Xi^* \f$ is the total equilibrium wrinkling combining the effects
+        of the flame inestability and turbulence interaction and is given by
+
+        \f[
+            \Xi^* = \frac {R}{R - G_\eta - G_{in}}
+        \f]
+
+    where:
+
+        \f$ G_\eta \f$ is the generation rate of wrinkling due to turbulence
+        interaction.
+
+        \f$ G_{in} = \kappa \rho_{u}/\rho_{b} \f$ is the generation
+         rate due to the flame inestability.
+
+    By adding the removal rates of the two effects:
+
+        \f[
+            R = G_\eta \frac{\Xi_{\eta_{eq}}}{\Xi_{\eta_{eq}} - 1}
+              + G_{in} \frac{\Xi_{{in}_{eq}}}{\Xi_{{in}_{eq}} - 1}
+        \f]
+
+    where:
+
+        \f$ R \f$ is the total removal.
+
+        \f$ G_\eta \f$ is a model constant.
+
+        \f$ \Xi_{\eta_{eq}} \f$ is the flame wrinkling due to turbulence.
+
+        \f$ \Xi_{{in}_{eq}} \f$ is the equilibrium level of the flame wrinkling
+        generated by inestability. It is a constant (default 2.5).
+
 
 SourceFiles
     XiModel.C
@@ -51,6 +102,8 @@ namespace Foam
                           Class XiModel Declaration
 \*---------------------------------------------------------------------------*/
 
+
+
 class XiModel
 {
 

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H

@@ -28,6 +28,33 @@ Class
 Description
     Laminar flame speed obtained from the SCOPE correlation.
 
+    Seven parameters are specified in terms of polynomial functions of
+    stoichiometry. Two polynomials are fitted, covering different parts of the
+    flammable range. If the mixture is outside the fitted range, linear
+    interpolation is used between the extreme of the polynomio and the upper or
+    lower flammable limit with the Markstein number constant.
+
+    Variations of pressure and temperature from the reference values are taken
+    into account through \f$ pexp \f$ and \f$ texp \f$
+
+    The laminar burning velocity fitting polynomial is:
+
+    \f$ Su = a_{0}(1+a_{1}x+K+..a_{i}x^{i}..+a_{6}x^{6}) (p/p_{ref})^{pexp}
+    (T/T_{ref})^{texp} \f$
+
+    where:
+
+        \f$ a_{i} \f$ are the polinomial coefficients.
+
+        \f$ pexp \f$ and \f$ texp \f$ are the pressure and temperature factors
+        respectively.
+
+        \f$ x \f$ is the equivalence ratio.
+
+        \f$ T_{ref} \f$ and \f$ p_{ref} \f$ are the temperature and pressure
+        references for the laminar burning velocity.
+
+
 SourceFiles
     SCOPELaminarFlameSpeed.C
 

applications/test/dictionary/dictionaryTest.C

@@ -45,6 +45,27 @@ int main(int argc, char *argv[])
 
     Info<< testDict << endl;
 
+    {
+        dictionary someDict;
+        someDict.add(keyType("a.*", true), "subdictValue");
+
+        dictionary dict;
+        dict.add("someDict", someDict);
+        dict.add(keyType(".*", true), "parentValue");
+
+        Info<< "dict:" << dict << endl;
+
+        // Wildcard find.
+        Info<< "Wildcard find \"abc\" in top directory : "
+            << dict.lookup("abc") << endl;
+        Info<< "Wildcard find \"abc\" in sub directory : "
+            << dict.subDict("someDict").lookup("abc")
+            << endl;
+        Info<< "Recursive wildcard find \"def\" in sub directory : "
+            << dict.subDict("someDict").lookup("def", true)
+            << endl;
+    }
+
     return 0;
 }
 

applications/utilities/mesh/conversion/polyDualMesh/Make/files

@@ -1,3 +1,4 @@
-polyDualMeshApp.C
+meshDualiser.C
+makePolyDualMesh.C
 
 EXE = $(FOAM_APPBIN)/polyDualMesh

applications/utilities/mesh/conversion/polyDualMesh/Make/options

@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/conversion/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools -lconversion
+    -lfiniteVolume \
+    -ldynamicMesh \
+    -lmeshTools

applications/utilities/mesh/conversion/polyDualMesh/makePolyDualMesh.C

@@ -0,0 +1,515 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Calculate the dual of a polyMesh. Adheres to all the feature&patch edges.
+
+    Feature angle:
+    convex features : point becomes single boundary cell with multiple
+                      boundary faces.
+    concave features: point becomes multiple boundary cells.
+
+    -splitAllFaces:
+    Normally only constructs a single face between two cells. This single face
+    might be too distorted. splitAllFaces will create a single face for every
+    original cell the face passes through. The mesh will thus have
+    multiple faces inbetween two cells! (so is not strictly upper-triangular
+    anymore - checkMesh will complain)
+    -doNotPreserveFaceZones:
+    By default all faceZones are preserved by marking all faces, edges and
+    points on them as features. The -doNotPreserveFaceZones disables this
+    behaviour.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "mathematicalConstants.H"
+#include "polyTopoChange.H"
+#include "mapPolyMesh.H"
+#include "PackedList.H"
+#include "meshTools.H"
+#include "OFstream.H"
+#include "meshDualiser.H"
+#include "ReadFields.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Naive feature detection. All boundary edges with angle > featureAngle become
+// feature edges. All points on feature edges become feature points. All
+// boundary faces become feature faces.
+void simpleMarkFeatures
+(
+    const polyMesh& mesh,
+    const PackedList<1>& isBoundaryEdge,
+    const scalar featureAngle,
+    const bool doNotPreserveFaceZones,
+
+    labelList& featureFaces,
+    labelList& featureEdges,
+    labelList& singleCellFeaturePoints,
+    labelList& multiCellFeaturePoints
+)
+{
+    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Working sets
+    labelHashSet featureEdgeSet;
+    labelHashSet singleCellFeaturePointSet;
+    labelHashSet multiCellFeaturePointSet;
+
+
+    // 1. Mark all edges between patches
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        const labelList& meshEdges = pp.meshEdges();
+
+        // All patch corner edges. These need to be feature points & edges!
+        for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
+        {
+            label meshEdgeI = meshEdges[edgeI];
+            featureEdgeSet.insert(meshEdgeI);
+            singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][0]);
+            singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][1]);
+        }
+    }
+
+
+
+    // 2. Mark all geometric feature edges
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Make distinction between convex features where the boundary point becomes
+    // a single cell and concave features where the boundary point becomes
+    // multiple 'half' cells.
+
+    // Addressing for all outside faces
+    primitivePatch allBoundary
+    (
+        SubList<face>
+        (
+            mesh.faces(),
+            mesh.nFaces()-mesh.nInternalFaces(),
+            mesh.nInternalFaces()
+        ),
+        mesh.points()
+    );
+
+    // Check for non-manifold points (surface pinched at point)
+    allBoundary.checkPointManifold(false, &singleCellFeaturePointSet);
+
+    // Check for non-manifold edges (surface pinched at edge)
+    const labelListList& edgeFaces = allBoundary.edgeFaces();
+    const labelList& meshPoints = allBoundary.meshPoints();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() > 2)
+        {
+            const edge& e = allBoundary.edges()[edgeI];
+
+            //Info<< "Detected non-manifold boundary edge:" << edgeI
+            //    << " coords:"
+            //    << allBoundary.points()[meshPoints[e[0]]]
+            //    << allBoundary.points()[meshPoints[e[1]]] << endl;
+
+            singleCellFeaturePointSet.insert(meshPoints[e[0]]);
+            singleCellFeaturePointSet.insert(meshPoints[e[1]]);
+        }
+    }
+
+    // Check for features.
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() == 2)
+        {
+            label f0 = eFaces[0];
+            label f1 = eFaces[1];
+
+            // check angle
+            const vector& n0 = allBoundary.faceNormals()[f0];
+            const vector& n1 = allBoundary.faceNormals()[f1];
+
+            if ((n0 & n1) < minCos)
+            {
+                const edge& e = allBoundary.edges()[edgeI];
+                label v0 = meshPoints[e[0]];
+                label v1 = meshPoints[e[1]];
+
+                label meshEdgeI = meshTools::findEdge(mesh, v0, v1);
+                featureEdgeSet.insert(meshEdgeI);
+
+                // Check if convex or concave by looking at angle
+                // between face centres and normal
+                vector c1c0
+                (
+                    allBoundary[f1].centre(allBoundary.points())
+                  - allBoundary[f0].centre(allBoundary.points())
+                );
+
+                if ((c1c0 & n0) > SMALL)
+                {
+                    // Found concave edge. Make into multiCell features
+                    Info<< "Detected concave feature edge:" << edgeI
+                        << " cos:" << (c1c0 & n0)
+                        << " coords:"
+                        << allBoundary.points()[v0]
+                        << allBoundary.points()[v1]
+                        << endl;
+
+                    singleCellFeaturePointSet.erase(v0);
+                    multiCellFeaturePointSet.insert(v0);
+                    singleCellFeaturePointSet.erase(v1);
+                    multiCellFeaturePointSet.insert(v1);
+                }
+                else
+                {
+                    // Convex. singleCell feature.
+                    if (!multiCellFeaturePointSet.found(v0))
+                    {
+                        singleCellFeaturePointSet.insert(v0);
+                    }
+                    if (!multiCellFeaturePointSet.found(v1))
+                    {
+                        singleCellFeaturePointSet.insert(v1);
+                    }
+                }
+            }
+        }
+    }
+
+
+    // 3. Mark all feature faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Face centres that need inclusion in the dual mesh
+    labelHashSet featureFaceSet(mesh.nFaces()-mesh.nInternalFaces());
+    // A. boundary faces.
+    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    {
+        featureFaceSet.insert(faceI);
+    }
+
+    // B. face zones.
+    const faceZoneMesh& faceZones = mesh.faceZones();
+
+    if (doNotPreserveFaceZones)
+    {
+        if (faceZones.size() > 0)
+        {
+            WarningIn("simpleMarkFeatures(..)")
+                << "Detected " << faceZones.size()
+                << " faceZones. These will not be preserved."
+                << endl;
+        }
+    }
+    else
+    {
+        if (faceZones.size() > 0)
+        {
+            Info<< "Detected " << faceZones.size()
+                << " faceZones. Preserving these by marking their"
+                << " points, edges and faces as features." << endl;
+        }
+
+        forAll(faceZones, zoneI)
+        {
+            const faceZone& fz = faceZones[zoneI];
+
+            Info<< "Inserting all faces in faceZone " << fz.name()
+                << " as features." << endl;
+
+            forAll(fz, i)
+            {
+                label faceI = fz[i];
+                const face& f = mesh.faces()[faceI];
+                const labelList& fEdges = mesh.faceEdges()[faceI];
+
+                featureFaceSet.insert(faceI);
+                forAll(f, fp)
+                {
+                    // Mark point as multi cell point (since both sides of
+                    // face should have different cells)
+                    singleCellFeaturePointSet.erase(f[fp]);
+                    multiCellFeaturePointSet.insert(f[fp]);
+
+                    // Make sure there are points on the edges.
+                    featureEdgeSet.insert(fEdges[fp]);
+                }
+            }
+        }
+    }
+
+    // Transfer to arguments
+    featureFaces = featureFaceSet.toc();
+    featureEdges = featureEdgeSet.toc();
+    singleCellFeaturePoints = singleCellFeaturePointSet.toc();
+    multiCellFeaturePoints = multiCellFeaturePointSet.toc();
+}
+
+
+// Dump features to .obj files
+void dumpFeatures
+(
+    const polyMesh& mesh,
+    const labelList& featureFaces,
+    const labelList& featureEdges,
+    const labelList& singleCellFeaturePoints,
+    const labelList& multiCellFeaturePoints
+)
+{
+    {
+        OFstream str("featureFaces.obj");
+        Info<< "Dumping centres of featureFaces to obj file " << str.name()
+            << endl;
+        forAll(featureFaces, i)
+        {
+            meshTools::writeOBJ(str, mesh.faceCentres()[featureFaces[i]]);
+        }
+    }
+    {
+        OFstream str("featureEdges.obj");
+        Info<< "Dumping featureEdges to obj file " << str.name() << endl;
+        label vertI = 0;
+
+        forAll(featureEdges, i)
+        {
+            const edge& e = mesh.edges()[featureEdges[i]];
+            meshTools::writeOBJ(str, mesh.points()[e[0]]);
+            vertI++;
+            meshTools::writeOBJ(str, mesh.points()[e[1]]);
+            vertI++;
+            str<< "l " << vertI-1 << ' ' << vertI << nl;
+        }
+    }
+    {
+        OFstream str("singleCellFeaturePoints.obj");
+        Info<< "Dumping featurePoints that become a single cell to obj file "
+            << str.name() << endl;
+        forAll(singleCellFeaturePoints, i)
+        {
+            meshTools::writeOBJ(str, mesh.points()[singleCellFeaturePoints[i]]);
+        }
+    }
+    {
+        OFstream str("multiCellFeaturePoints.obj");
+        Info<< "Dumping featurePoints that become multiple cells to obj file "
+            << str.name() << endl;
+        forAll(multiCellFeaturePoints, i)
+        {
+            meshTools::writeOBJ(str, mesh.points()[multiCellFeaturePoints[i]]);
+        }
+    }
+}
+
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+#   include "addTimeOptions.H"
+    argList::validArgs.append("feature angle[0-180]");
+    argList::validOptions.insert("splitAllFaces", "");
+    argList::validOptions.insert("doNotPreserveFaceZones", "");
+    argList::validOptions.insert("overwrite", "");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+    // Get times list
+    instantList Times = runTime.times();
+#   include "checkTimeOptions.H"
+    runTime.setTime(Times[startTime], startTime);
+
+#   include "createMesh.H"
+
+    // Mark boundary edges and points.
+    // (Note: in 1.4.2 we can use the built-in mesh point ordering
+    //  facility instead)
+    PackedList<1> isBoundaryEdge(mesh.nEdges());
+    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    {
+        const labelList& fEdges = mesh.faceEdges()[faceI];
+
+        forAll(fEdges, i)
+        {
+            isBoundaryEdge.set(fEdges[i], 1);
+        }
+    }
+
+    scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
+
+    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
+
+    Info<< "Feature:" << featureAngle << endl
+        << "minCos :" << minCos << endl
+        << endl;
+
+
+    const bool splitAllFaces = args.options().found("splitAllFaces");
+    const bool overwrite = args.options().found("overwrite");
+    const bool doNotPreserveFaceZones = args.options().found
+    (
+        "doNotPreserveFaceZones"
+    );
+
+    // Face(centre)s that need inclusion in the dual mesh
+    labelList featureFaces;
+    // Edge(centre)s  ,,
+    labelList featureEdges;
+    // Points (that become a single cell) that need inclusion in the dual mesh
+    labelList singleCellFeaturePoints;
+    // Points (that become a mulitple cells)        ,,
+    labelList multiCellFeaturePoints;
+
+    // Sample implementation of feature detection.
+    simpleMarkFeatures
+    (
+        mesh,
+        isBoundaryEdge,
+        featureAngle,
+        doNotPreserveFaceZones,
+
+        featureFaces,
+        featureEdges,
+        singleCellFeaturePoints,
+        multiCellFeaturePoints
+    );
+
+    // If we want to split all polyMesh faces into one dualface per cell
+    // we are passing through we also need a point
+    // at the polyMesh facecentre and edgemid of the faces we want to
+    // split.
+    if (splitAllFaces)
+    {
+        featureEdges = identity(mesh.nEdges());
+        featureFaces = identity(mesh.nFaces());
+    }
+
+    // Write obj files for debugging
+    dumpFeatures
+    (
+        mesh,
+        featureFaces,
+        featureEdges,
+        singleCellFeaturePoints,
+        multiCellFeaturePoints
+    );
+
+
+
+    // Read objects in time directory
+    IOobjectList objects(mesh, runTime.timeName());
+
+    // Read vol fields.
+    PtrList<volScalarField> vsFlds;
+    ReadFields(mesh, objects, vsFlds);
+
+    PtrList<volVectorField> vvFlds;
+    ReadFields(mesh, objects, vvFlds);
+
+    PtrList<volSphericalTensorField> vstFlds;
+    ReadFields(mesh, objects, vstFlds);
+
+    PtrList<volSymmTensorField> vsymtFlds;
+    ReadFields(mesh, objects, vsymtFlds);
+
+    PtrList<volTensorField> vtFlds;
+    ReadFields(mesh, objects, vtFlds);
+
+    // Read surface fields.
+    PtrList<surfaceScalarField> ssFlds;
+    ReadFields(mesh, objects, ssFlds);
+
+    PtrList<surfaceVectorField> svFlds;
+    ReadFields(mesh, objects, svFlds);
+
+    PtrList<surfaceSphericalTensorField> sstFlds;
+    ReadFields(mesh, objects, sstFlds);
+
+    PtrList<surfaceSymmTensorField> ssymtFlds;
+    ReadFields(mesh, objects, ssymtFlds);
+
+    PtrList<surfaceTensorField> stFlds;
+    ReadFields(mesh, objects, stFlds);
+
+
+    // Topo change container
+    polyTopoChange meshMod(mesh.boundaryMesh().size());
+
+    // Mesh dualiser engine
+    meshDualiser dualMaker(mesh);
+
+    // Insert all commands into polyTopoChange to create dual of mesh. This does
+    // all the hard work.
+    dualMaker.setRefinement
+    (
+        splitAllFaces,
+        featureFaces,
+        featureEdges,
+        singleCellFeaturePoints,
+        multiCellFeaturePoints,
+        meshMod
+    );
+
+    // Create mesh, return map from old to new mesh.
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+    // Update fields
+    mesh.updateMesh(map);
+
+    // Optionally inflate mesh
+    if (map().hasMotionPoints())
+    {
+        mesh.movePoints(map().preMotionPoints());
+    }
+
+    if (!overwrite)
+    {
+        runTime++;
+        mesh.setInstance(runTime.timeName());
+    }
+
+    Info<< "Writing dual mesh to " << runTime.timeName() << endl;
+
+    mesh.write();
+ 
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.H

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    meshDualiser
+
+Description
+    Creates dual of polyMesh. Every point becomes a cell (or multiple cells
+    for feature points), a walk around every edge creates faces between them.
+
+    Put all points you want in the final mesh into featurePoints; all edge(mid)s
+    you want in the final mesh into featureEdges; all face(centre)s in
+    faceFaces.
+
+    Usually to preserve boundaries:
+        - all boundary faces are featureFaces
+        - all edges and points inbetween different patches are
+          featureEdges/points.
+
+    In same way you can also preserve internal faces (e.g. faceZones)
+
+SourceFiles
+    meshDualiser.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshDualiser_H
+#define meshDualiser_H
+
+#include "DynamicList.H"
+#include "PackedList.H"
+#include "boolList.H"
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class polyTopoChange;
+
+/*---------------------------------------------------------------------------*\
+                           Class meshDualiser Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshDualiser
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        //- From point on cell to dual cell. Either single entry or
+        //  one entry per pointCells
+        labelListList pointToDualCells_;
+
+        //- From point to dual point (or -1 if not feature point).
+        labelList pointToDualPoint_;
+
+        //- From cell to dual point. All cells become point
+        labelList cellToDualPoint_;
+
+        //- From face to dual point (or -1 if not feature face)
+        labelList faceToDualPoint_;
+
+        //- From edge to dual point (or -1 if not feature edge)
+        labelList edgeToDualPoint_;
+
+
+    // Private Member Functions
+
+        static void checkPolyTopoChange(const polyTopoChange&);
+
+        static void dumpPolyTopoChange(const polyTopoChange&, const fileName&);
+
+        //- Find dual cell given point and cell
+        label findDualCell(const label cellI, const label pointI) const;
+
+        //- Helper function to generate dualpoints on all boundary edges
+        //  emanating from (boundary & feature) point
+        void generateDualBoundaryEdges
+        (
+            const PackedList<1>&,
+            const label pointI,
+            polyTopoChange&
+        );
+
+        //- Check that owner and neighbour of face have same dual cell
+        bool sameDualCell
+        (
+            const label faceI,
+            const label pointI
+        ) const;
+
+        //- Add internal face
+        label addInternalFace
+        (
+            const label masterPointI,
+            const label masterEdgeI,
+            const label masterFaceI,
+
+            const bool edgeOrder,
+            const label dualCell0,
+            const label dualCell1,
+            const DynamicList<label>& verts,
+            polyTopoChange& meshMod
+        ) const;
+
+        //- Add boundary face
+        label addBoundaryFace
+        (
+            const label masterPointI,
+            const label masterEdgeI,
+            const label masterFaceI,
+
+            const label dualCellI,
+            const label patchI,
+            const DynamicList<label>& verts,
+            polyTopoChange& meshMod
+        ) const;
+
+        //- Create internal faces walking around edge
+        void createFacesAroundEdge
+        (
+            const bool splitFace,
+            const PackedList<1>&,
+            const label edgeI,
+            const label startFaceI,
+            polyTopoChange&,
+            boolList& doneEFaces
+        ) const;
+
+        //- Create single internal face from internal face
+        void createFaceFromInternalFace
+        (
+            const label faceI,
+            label& fp,
+            polyTopoChange&
+        ) const;
+
+        //- Creates boundary faces walking around point on patchI.
+        void createFacesAroundBoundaryPoint
+        (
+            const label patchI,
+            const label patchPointI,
+            const label startFaceI,
+            polyTopoChange&,
+            boolList& donePFaces            // pFaces visited
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        meshDualiser(const meshDualiser&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const meshDualiser&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("meshDualiser");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        meshDualiser(const polyMesh&);
+
+
+    // Member Functions
+
+        // Access
+
+            //- From point on cell to dual cell. Either single entry or
+            //  one entry per pointCells.
+            const labelListList& pointToDualCells() const
+            {
+                return pointToDualCells_;
+            }
+
+            //- From point to dual point (or -1 if not feature point).
+            const labelList& pointToDualPoint() const
+            {
+                return pointToDualPoint_;
+            }
+
+            //- From cell to dual point (at cell centre). All cells become
+            //  points.
+            const labelList& cellToDualPoint() const
+            {
+                return cellToDualPoint_;
+            }
+
+            //- From face to dual point (at face centre; or -1 if not
+            //  feature face).
+            const labelList& faceToDualPoint() const
+            {
+                return faceToDualPoint_;
+            }
+
+            //- From edge to dual point (at edge mid; or -1 if not feature
+            //  edge).
+            const labelList& edgeToDualPoint() const
+            {
+                return edgeToDualPoint_;
+            }
+
+
+        // Edit
+
+
+            //- Insert all changes into meshMod to convert the polyMesh into
+            //  its dual.
+            //  featureFaces  : faces where we want a point at the face centre
+            //  featureEdges  : edges           ,,                 edge mid
+            //  featurePoints : points          ,,    point. Two variants:
+            //      singleCellFeaturePoints : point becomes one dualcell.
+            //          Use this for e.g. convex boundary points.
+            //      multiCellFeaturePoints : one dualcell per original cell
+            //          around point. Use this for e.g. concave boundary points
+            //          since it prevents big concave boundary cells.
+            void setRefinement
+            (
+                const bool splitFace,
+                const labelList& featureFaces,
+                const labelList& featureEdges,
+                const labelList& singleCellFeaturePoints,
+                const labelList& multiCellFeaturePoints,
+                polyTopoChange& meshMod
+            );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C

@@ -1,79 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Application
-    polyDualMesh
-
-Description
-    Calculate the dual of a polyMesh. Adheres to all the feature&patch edges
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "Time.H"
-#include "polyDualMesh.H"
-#include "mathematicalConstants.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::noParallel();
-    argList::validArgs.append("feature angle[0-180]");
-    argList::validOptions.insert("overwrite", "");
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-    runTime.functionObjects().off();
-#   include "createPolyMesh.H"
-
-    scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
-    bool overwrite = args.options().found("overwrite");
-
-    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
-
-    Info<< "Feature:" << featureAngle << endl
-        << "minCos :" << minCos << endl
-        << endl;
-
-    polyDualMesh dualMesh(mesh, minCos);
-
-    if (!overwrite)
-    {
-        runTime++;
-    }
-
-    Pout<< "Writing dualMesh to " << runTime.timeName() << endl;
-
-    dualMesh.write();
-
-    Info << "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C

@@ -617,8 +617,10 @@ void Foam::vtkPV3Foam::addPatchNames(vtkRenderer* renderer)
         {
             const labelList& eFaces = edgeFaces[edgeI];
 
-            if (eFaces.size() != 2)
+            if (eFaces.size() == 1)
             {
+                // Note: could also do ones with > 2 faces but this gives
+                // too many zones for baffles
                 featEdge[edgeI] = true;
             }
             else if (mag(n[eFaces[0]] & n[eFaces[1]]) < 0.5)

bin/doxyFilt

@@ -50,9 +50,9 @@ then
    */applications/solvers/*.C | */applications/utilities/*.C )
       awkScript=$WM_PROJECT_DIR/bin/tools/doxyFilt-top.awk
       ;;
-   */applications/solvers/*.H | */applications/utilities/*.H )
+#   */applications/solvers/*.H | */applications/utilities/*.H )
-      awkScript=$WM_PROJECT_DIR/bin/tools/doxyFilt-ignore.awk
+#      awkScript=$WM_PROJECT_DIR/bin/tools/doxyFilt-ignore.awk
-      ;;
+#      ;;
    esac
 
    awk -f $awkScript $1 | \

doc/Doxygen/Doxyfile

@@ -14,6 +14,16 @@
 # Project related configuration options
 #---------------------------------------------------------------------------
 
+
+#--------------------------------
+# PATH FOR OPEN CFD LATEX MACROS
+#-------------------------------
+
+@INLUDE_PATH = $(TEXINPUTS)
+@INLUDE_PATH += $(BIBINPUTS)
+@INLUDE_PATH += $(BSTINPUTS)
+
+
 # This tag specifies the encoding used for all characters in the config file that
 # follow. The default is UTF-8 which is also the encoding used for all text before
 # the first occurrence of this tag. Doxygen uses libiconv (or the iconv built into
@@ -477,9 +487,11 @@ WARN_LOGFILE           =
 # directories like "/usr/src/myproject". Separate the files or directories
 # with spaces.
 
-INPUT                  = $(WM_PROJECT_DIR)/src \
+INPUT              = $(WM_PROJECT_DIR)/src \
-                         $(WM_PROJECT_DIR)/applications/utilities \
+                     $(WM_PROJECT_DIR)/applications/utilities \
-                         $(WM_PROJECT_DIR)/applications/solvers
+                     $(WM_PROJECT_DIR)/applications/solvers
+
+
 
 # This tag can be used to specify the character encoding of the source files that
 # doxygen parses. Internally doxygen uses the UTF-8 encoding, which is also the default
@@ -536,7 +548,9 @@ EXCLUDE_SYMBOLS        =
 # directories that contain example code fragments that are included (see
 # the \include command).
 
-EXAMPLE_PATH           =
+EXAMPLE_PATH           = $(TEXINPUTS) \
+ 						 $(BIBINPUTS) \
+						 $(BSTINPUTS)
 
 # If the value of the EXAMPLE_PATH tag contains directories, you can use the
 # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
@@ -824,7 +838,8 @@ PAPER_TYPE             = a4wide
 # The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
 # packages that should be included in the LaTeX output.
 
-EXTRA_PACKAGES         =
+EXTRA_PACKAGES         = conditionalEqns finiteVolume algorithmic tensorCommon \
+						 tensorOperator tensorEquation 
 
 # The LATEX_HEADER tag can be used to specify a personal LaTeX header for
 # the generated latex document. The header should contain everything until

etc/settings.csh

@@ -187,7 +187,7 @@ case MPICH-GM:
     setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/mpich-gm
     breaksw
 
-case MPICH-GM:
+case HPMPI:
     setenv MPI_HOME /opt/hpmpi
     setenv MPI_ARCH_PATH $MPI_HOME
     setenv MPICH_ROOT=$MPI_ARCH_PATH

src/OSspecific/Unix/regularExpression.H

@@ -37,6 +37,8 @@ SourceFiles
 
 #include <sys/types.h>
 #include <regex.h>
+#include "string.H"
+#include "IOstreams.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +75,7 @@ public:
         {
             preg_ = new regex_t;
 
-            if (regcomp(preg_, s.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
+            if (regcomp(preg_, s.c_str(), REG_EXTENDED) != 0)
             {
                 FatalErrorIn
                 (
@@ -102,10 +104,12 @@ public:
         //- Matches?
         inline bool matches(const string& s) const
         {
-            size_t nmatch = 0;
+            size_t nmatch = 1;
-            regmatch_t *pmatch = NULL;
+            regmatch_t pmatch[1];
 
-            return regexec(preg_, s.c_str(), nmatch, pmatch, 0) == 0;
+            int errVal = regexec(preg_, s.c_str(), nmatch, pmatch, 0);
+
+            return (errVal == 0 && pmatch[0].rm_eo == label(s.size()));
         }
 };
 

src/OpenFOAM/db/dictionary/dictionary.C

@@ -27,6 +27,7 @@ License
 #include "dictionary.H"
 #include "primitiveEntry.H"
 #include "dictionaryEntry.H"
+#include "regularExpression.H"
 
 /* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
 

src/OpenFOAM/db/dictionary/dictionary.H

@@ -60,7 +60,6 @@ SourceFiles
 #include "HashTable.H"
 #include "wordList.H"
 #include "className.H"
-#include "regularExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,7 +67,7 @@ namespace Foam
 {
 
 // Forward declaration of friend functions and operators
-
+class regularExpression;
 class dictionary;
 Istream& operator>>(Istream&, dictionary&);
 Ostream& operator<<(Ostream&, const dictionary&);

src/OpenFOAM/db/dictionary/dictionaryIO.C

@@ -27,6 +27,7 @@ License
 #include "dictionary.H"
 #include "IFstream.H"
 #include "inputModeEntry.H"
+#include "regularExpression.H"
 
 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
 

src/OpenFOAM/fields/Fields/transformField/transformField.C

@@ -26,6 +26,7 @@ License
 
 #include "transformField.H"
 #include "FieldM.H"
+#include "diagTensor.H"
 
 // * * * * * * * * * * * * * * * global functions  * * * * * * * * * * * * * //
 
@@ -75,7 +76,8 @@ void Foam::transform
 {
     vector T = tr.t();
 
-    if (mag(tr.r().w()) > SMALL)
+    // Check if any rotation
+    if (mag(tr.r().R() - I) > SMALL)
     {
         transform(rtf, tr.r(), tf);
 
@@ -90,6 +92,10 @@ void Foam::transform
         {
             TFOR_ALL_F_OP_S_OP_F(vector, rtf, =, vector, T, +, vector, tf);
         }
+        else
+        {
+            rtf = vector::zero;
+        }
     }
 }
 

src/engine/ignition/ignitionSite.C

@@ -37,14 +37,17 @@ namespace Foam
 
 void ignitionSite::findIgnitionCells(const fvMesh& mesh)
 {
+    // Bit tricky: generate C and V before shortcutting if cannot find
+    // cell locally. mesh.C generation uses parallel communication.
+    const volVectorField& centres = mesh.C();
+    const scalarField& vols = mesh.V();
+
     label ignCell = mesh.findCell(location_);
     if (ignCell == -1)
     {
         return;
     }
 
-    const volVectorField& centres = mesh.C();
-    const scalarField& vols = mesh.V();
     scalar radius = diameter_/2.0;
 
     cells_.setSize(1);

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

@@ -49,6 +49,22 @@ void Foam::setRefCell
             if (Pstream::master())
             {
                 refCelli = readLabel(dict.lookup(refCellName));
+
+                if (refCelli < 0 || refCelli >= field.mesh().nCells())
+                {
+                    FatalErrorIn
+                    (
+                        "void Foam::setRefCell"
+                         "("
+                         "    const volScalarField&,"
+                         "    const dictionary&,"
+                         "    label& scalar&,"
+                         "    bool"
+                         ")"
+                    )   << "Illegal master cellID " << refCelli
+                        << ". Should be 0.." << field.mesh().nCells()
+                        << exit(FatalError);
+                }
             }
             else
             {
@@ -75,7 +91,7 @@ void Foam::setRefCell
                 )
                   << "Unable to set reference cell for field " << field.name()
                   << nl << "    Reference point " << refPointName
-                  << " found on multiple domains" << nl << abort(FatalError);
+                  << " found on multiple domains" << nl << exit(FatalError);
             }
         }
         else
@@ -92,7 +108,7 @@ void Foam::setRefCell
             )
               << "Unable to set reference cell for field" << field.name() << nl
               << "    Please supply either " << refCellName
-              << " or " << refPointName << nl << abort(FatalError);
+              << " or " << refPointName << nl << exit(FatalError);
         }
 
         refValue = readScalar(dict.lookup(refValueName));

src/meshTools/searchableSurface/searchableSphere.C

@@ -292,7 +292,8 @@ void Foam::searchableSphere::getNormal
         if (info[i].hit())
         {
             normal[i] = info[i].hitPoint() - centre_;
-            normal[i] /= mag(normal[i]);
+
+            normal[i] /= mag(normal[i])+VSMALL;
         }
         else
         {

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -30,7 +30,8 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurfaceCell.H"
+#include "isoSurface.H"
+//#include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,19 +43,46 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::distanceSurface::createGeometry() const
+void Foam::distanceSurface::createGeometry()
 {
     // Clear any stored topo
     facesPtr_.clear();
 
+
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
     // Distance to cell centres
-    scalarField cellDistance(mesh().nCells());
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    cellDistancePtr_.reset
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "cellDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm,
+            dimensionedScalar("zero", dimless/dimTime, 0)
+        )
+    );
+    volScalarField& cellDistance = cellDistancePtr_();
+
+    // Internal field
     {
+        const pointField& cc = fvm.C();
+        scalarField& fld = cellDistance.internalField();
+
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
         (
-            mesh().cellCentres(),
+            cc,
-            scalarField(mesh().nCells(), GREAT),
+            scalarField(cc.size(), GREAT),
             nearest
         );
 
@@ -63,98 +91,109 @@ void Foam::distanceSurface::createGeometry() const
             vectorField normal;
             surfPtr_().getNormal(nearest, normal);
 
-            forAll(cellDistance, cellI)
+            forAll(nearest, i)
             {
-                vector d(mesh().cellCentres()[cellI]-nearest[cellI].hitPoint());
+                vector d(cc[i]-nearest[i].hitPoint());
 
-                if ((d&normal[cellI]) > 0)
+                if ((d&normal[i]) > 0)
                 {
-                    cellDistance[cellI] = Foam::mag(d);
+                    fld[i] = Foam::mag(d);
                 }
                 else
                 {
-                    cellDistance[cellI] = -Foam::mag(d);
+                    fld[i] = -Foam::mag(d);
                 }
             }
         }
         else
         {
-            forAll(cellDistance, cellI)
+            forAll(nearest, i)
             {
-                cellDistance[cellI] = Foam::mag
+                fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                (
-                    nearest[cellI].hitPoint()
-                  - mesh().cellCentres()[cellI]
-                );
             }
         }
     }
 
+    // Patch fields
+    {
+        forAll(fvm.C().boundaryField(), patchI)
+        {
+            const pointField& cc = fvm.C().boundaryField()[patchI];
+            fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
+
+            List<pointIndexHit> nearest;
+            surfPtr_().findNearest
+            (
+                cc,
+                scalarField(cc.size(), GREAT),
+                nearest
+            );
+
+            if (signed_)
+            {
+                vectorField normal;
+                surfPtr_().getNormal(nearest, normal);
+
+                forAll(nearest, i)
+                {
+                    vector d(cc[i]-nearest[i].hitPoint());
+
+                    if ((d&normal[i]) > 0)
+                    {
+                        fld[i] = Foam::mag(d);
+                    }
+                    else
+                    {
+                        fld[i] = -Foam::mag(d);
+                    }
+                }
+            }
+            else
+            {
+                forAll(nearest, i)
+                {
+                    fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+                }
+            }
+        }
+    }
+
+    // On processor patches the mesh.C() will already be the cell centre
+    // on the opposite side so no need to swap cellDistance.
+
+
     // Distance to points
-    scalarField pointDistance(mesh().nPoints());
+    pointDistance_.setSize(fvm.nPoints());
     {
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
         (
-            mesh().points(),
+            fvm.points(),
-            scalarField(mesh().nPoints(), GREAT),
+            scalarField(fvm.nPoints(), GREAT),
             nearest
         );
-        forAll(pointDistance, pointI)
+        forAll(pointDistance_, pointI)
         {
-            pointDistance[pointI] = Foam::mag
+            pointDistance_[pointI] = Foam::mag
             (
                 nearest[pointI].hitPoint()
-              - mesh().points()[pointI]
+              - fvm.points()[pointI]
             );
         }
     }
 
     //- Direct from cell field and point field.
-    const isoSurfaceCell iso
+    isoSurfPtr_.reset
     (
-        mesh(),
+        new isoSurface
-        cellDistance,
+        (
-        pointDistance,
+            cellDistance,
-        distance_,
+            pointDistance_,
-        regularise_
+            distance_,
+            regularise_
+        )
     );
 
-    ////- From point field and interpolated cell.
-    //scalarField cellAvg(mesh().nCells(), scalar(0.0));
-    //labelField nPointCells(mesh().nCells(), 0);
-    //{
-    //    for (label pointI = 0; pointI < mesh().nPoints(); pointI++)
-    //    {
-    //        const labelList& pCells = mesh().pointCells(pointI);
-    //
-    //        forAll(pCells, i)
-    //        {
-    //            label cellI = pCells[i];
-    //
-    //            cellAvg[cellI] += pointDistance[pointI];
-    //            nPointCells[cellI]++;
-    //        }
-    //    }
-    //}
-    //forAll(cellAvg, cellI)
-    //{
-    //    cellAvg[cellI] /= nPointCells[cellI];
-    //}
-    //
-    //const isoSurface iso
-    //(
-    //    mesh(),
-    //    cellAvg,
-    //    pointDistance,
-    //    distance_,
-    //    regularise_
-    //);
-
-
-    const_cast<distanceSurface&>(*this).triSurface::operator=(iso);
-    meshCells_ = iso.meshCells();
-
     if (debug)
     {
         print(Pout);
@@ -193,9 +232,8 @@ Foam::distanceSurface::distanceSurface
     signed_(readBool(dict.lookup("signed"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
     zoneName_(word::null),
-    facesPtr_(NULL),
+    isoSurfPtr_(NULL),
-    storedTimeIndex_(-1),
+    facesPtr_(NULL)
-    meshCells_(0)
 {
 //    label zoneId = -1;
 //    if (dict.readIfPresent("zone", zoneName_))

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -37,8 +37,9 @@ SourceFiles
 #define distanceSurface_H
 
 #include "sampledSurface.H"
-#include "triSurface.H"
 #include "searchableSurface.H"
+//#include "isoSurfaceCell.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,8 +52,7 @@ namespace Foam
 
 class distanceSurface
 :
-    public sampledSurface,
+    public sampledSurface
-    public triSurface
 {
     // Private data
 
@@ -71,23 +71,24 @@ class distanceSurface
         //- zone name (if restricted to zones)
         word zoneName_;
 
+
+        //- Distance to cell centres
+        autoPtr<volScalarField> cellDistancePtr_;
+
+        //- Distance to points
+        scalarField pointDistance_;
+
+        //- Constructed iso surface
+        autoPtr<isoSurface> isoSurfPtr_;
+
         //- triangles converted to faceList
         mutable autoPtr<faceList> facesPtr_;
 
 
-        // Recreated for every isoSurface
-
-            //- Time at last call
-            mutable label storedTimeIndex_;
-
-            //- For every triangle the original cell in mesh
-            mutable labelList meshCells_;
-
-
     // Private Member Functions
 
-        //- Create iso surface (if time has changed)
+        //- Create iso surface
-        void createGeometry() const;
+        void createGeometry();
 
         //- sample field on faces
         template <class Type>
@@ -129,7 +130,7 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return triSurface::points();
+            return surface().points();
         }
 
         //- Faces of surface
@@ -137,7 +138,7 @@ public:
         {
             if (!facesPtr_.valid())
             {
-                const triSurface& s = *this;
+                const triSurface& s = surface();
 
                 facesPtr_.reset(new faceList(s.size()));
 
@@ -153,6 +154,10 @@ public:
         //- Correct for mesh movement and/or field changes
         virtual void correct(const bool meshChanged);
 
+        const isoSurface& surface() const
+        {
+            return isoSurfPtr_();
+        }
 
         //- sample field on surface
         virtual tmp<scalarField> sample

src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C

@@ -39,7 +39,7 @@ Foam::distanceSurface::sampleField
     const GeometricField<Type, fvPatchField, volMesh>& vField
 ) const
 {
-    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+    return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
 }
 
 
@@ -50,33 +50,25 @@ Foam::distanceSurface::interpolateField
     const interpolation<Type>& interpolator
 ) const
 {
-    // One value per point
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
 
-    boolList pointDone(points().size(), false);
+    // Get fields to sample. Assume volPointInterpolation!
+    const GeometricField<Type, fvPatchField, volMesh>& volFld =
+        interpolator.psi();
 
-    forAll(faces(), cutFaceI)
+    tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld
-    {
+    (
-        const face& f = faces()[cutFaceI];
+        volPointInterpolation::New(fvm).interpolate(volFld)
+    );
 
-        forAll(f, faceVertI)
+    // Sample.
-        {
+    return surface().interpolate
-            label pointI = f[faceVertI];
+    (
-
+        cellDistancePtr_(),
-            if (!pointDone[pointI])
+        pointDistance_,
-            {
+        volFld,
-                values[pointI] = interpolator.interpolate
+        pointFld()
-                (
+    );
-                    points()[pointI],
-                    meshCells_[cutFaceI]
-                );
-                pointDone[pointI] = true;
-            }
-        }
-    }
-
-    return tvalues;
 }
 
 

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -108,17 +108,15 @@ void Foam::isoSurface::calcCutTypes
         const polyPatch& pp = patches[patchI];
 
         label faceI = pp.start();
-        forAll(pp, i)
-        {
-            bool ownLower = (cVals[own[faceI]] < iso_);
-            bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
 
-            if (ownLower != neiLower)
+        if (isA<emptyPolyPatch>(pp))
-            {
+        {
-                faceCutType_[faceI] = CUT;
+            // Assume zero gradient so owner and neighbour/boundary value equal
-            }
+
-            else
+            forAll(pp, i)
             {
+                bool ownLower = (cVals[own[faceI]] < iso_);
+
                 // Mesh edge.
                 const face f = mesh_.faces()[faceI];
 
@@ -129,15 +127,48 @@ void Foam::isoSurface::calcCutTypes
                     (
                         (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
                      || (fpLower != ownLower)
-                     || (fpLower != neiLower)
                     )
                     {
                         faceCutType_[faceI] = CUT;
                         break;
                     }
                 }
+                faceI++;
+            }
+        }
+        else
+        {
+            forAll(pp, i)
+            {
+                bool ownLower = (cVals[own[faceI]] < iso_);
+                bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
+
+                if (ownLower != neiLower)
+                {
+                    faceCutType_[faceI] = CUT;
+                }
+                else
+                {
+                    // Mesh edge.
+                    const face f = mesh_.faces()[faceI];
+
+                    forAll(f, fp)
+                    {
+                        bool fpLower = (pVals[f[fp]] < iso_);
+                        if
+                        (
+                            (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
+                         || (fpLower != ownLower)
+                         || (fpLower != neiLower)
+                        )
+                        {
+                            faceCutType_[faceI] = CUT;
+                            break;
+                        }
+                    }
+                }
+                faceI++;
             }
-            faceI++;
         }
     }
 
@@ -328,8 +359,17 @@ void Foam::isoSurface::getNeighbour
         label patchI = boundaryRegion[bFaceI];
         label patchFaceI = faceI-mesh_.boundaryMesh()[patchI].start();
 
-        nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+        if (isA<emptyPolyPatch>(mesh_.boundaryMesh()[patchI]))
-        nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+        {
+            // Assume zero gradient
+            nbrValue = cVals[own[faceI]];
+            nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+        }
+        else
+        {
+            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+            nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+        }
     }
 }
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

@@ -201,6 +201,34 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
 }
 
 
+void Foam::isoSurfaceCell::calcCutTypes
+(
+    const PackedList<1>& isTet,
+    const scalarField& cVals,
+    const scalarField& pVals
+)
+{
+    cellCutType_.setSize(mesh_.nCells());
+    nCutCells_ = 0;
+    forAll(mesh_.cells(), cellI)
+    {
+        cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI);
+
+        if (cellCutType_[cellI] == CUT)
+        {
+            nCutCells_++;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : detected " << nCutCells_
+            << " candidate cut cells." << endl;
+    }
+}
+
+
+
 // Return the two common points between two triangles
 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
 (
@@ -325,7 +353,6 @@ void Foam::isoSurfaceCell::calcSnappedCc
     const PackedList<1>& isTet,
     const scalarField& cVals,
     const scalarField& pVals,
-    const List<cellCutType>& cellCutType,
 
     DynamicList<point>& snappedPoints,
     labelList& snappedCc
@@ -344,7 +371,7 @@ void Foam::isoSurfaceCell::calcSnappedCc
 
     forAll(mesh_.cells(), cellI)
     {
-        if (cellCutType[cellI] == CUT && isTet.get(cellI) == 0)
+        if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0)
         {
             scalar cVal = cVals[cellI];
 
@@ -604,7 +631,6 @@ void Foam::isoSurfaceCell::calcSnappedPoint
     const PackedList<1>& isTet,
     const scalarField& cVals,
     const scalarField& pVals,
-    const List<cellCutType>& cellCutType,
 
     DynamicList<point>& snappedPoints,
     labelList& snappedPoint
@@ -635,10 +661,10 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
             if
             (
-                cellCutType[mesh_.faceOwner()[faceI]] == CUT
+                cellCutType_[mesh_.faceOwner()[faceI]] == CUT
              || (
                     mesh_.isInternalFace(faceI)
-                 && cellCutType[mesh_.faceNeighbour()[faceI]] == CUT
+                 && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT
                 )
             )
             {
@@ -766,171 +792,6 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 }
 
 
-Foam::point Foam::isoSurfaceCell::generatePoint
-(
-    const DynamicList<point>& snappedPoints,
-
-    const scalar s0,
-    const point& p0,
-    const label p0Index,
-
-    const scalar s1,
-    const point& p1,
-    const label p1Index
-) const
-{
-    scalar d = s1-s0;
-
-    if (mag(d) > VSMALL)
-    {
-        scalar s = (iso_-s0)/d;
-
-        if (s >= 0.5 && s <= 1 && p1Index != -1)
-        {
-            return snappedPoints[p1Index];
-        }
-        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
-        {
-            return snappedPoints[p0Index];
-        }
-        else
-        {
-            return s*p1 + (1.0-s)*p0;
-        }
-    }
-    else
-    {
-        scalar s = 0.4999;
-
-        return s*p1 + (1.0-s)*p0;
-    }
-}
-
-
-void Foam::isoSurfaceCell::generateTriPoints
-(
-    const DynamicList<point>& snapped,
-
-    const scalar s0,
-    const point& p0,
-    const label p0Index,
-
-    const scalar s1,
-    const point& p1,
-    const label p1Index,
-
-    const scalar s2,
-    const point& p2,
-    const label p2Index,
-
-    const scalar s3,
-    const point& p3,
-    const label p3Index,
-
-    DynamicList<point>& points
-) const
-{
-    int triIndex = 0;
-    if (s0 < iso_)
-    {
-        triIndex |= 1;
-    }
-    if (s1 < iso_)
-    {
-        triIndex |= 2;
-    }
-    if (s2 < iso_)
-    {
-        triIndex |= 4;
-    }
-    if (s3 < iso_)
-    {
-        triIndex |= 8;
-    }
-
-    /* Form the vertices of the triangles for each case */
-    switch (triIndex)
-    {
-        case 0x00:
-        case 0x0F:
-        break;
-
-        case 0x0E:
-        case 0x01:
-            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-        break;
-
-        case 0x0D:
-        case 0x02:
-            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-        break;
-
-        case 0x0C:
-        case 0x03:
-        {
-            point tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            point tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
-
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp2);
-            points.append(tp2);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
-        }
-        break;
-
-        case 0x0B:
-        case 0x04:
-        {
-            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
-        }
-        break;
-
-        case 0x0A:
-        case 0x05:
-        {
-            point tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            point tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            points.append(tp0);
-            points.append(tp1);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
-        }
-        break;
-
-        case 0x09:
-        case 0x06:
-        {
-            point tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            point tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp0);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(tp1);
-        }
-        break;
-
-        case 0x07:
-        case 0x08:
-            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
-        break;
-    }
-}
 
 
 Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
@@ -1589,26 +1450,9 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
 
     // Determine if any cut through cell
+    calcCutTypes(isTet, cVals, pVals);
 
-    List<cellCutType> cellCutType(mesh_.nCells());
+    DynamicList<point> snappedPoints(nCutCells_);
-    label nCutCells = 0;
-    forAll(mesh_.cells(), cellI)
-    {
-        cellCutType[cellI] = calcCutType(isTet, cVals, pVals, cellI);
-
-        if (cellCutType[cellI] == CUT)
-        {
-            nCutCells++;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : detected " << nCutCells
-            << " candidate cut cells." << endl;
-    }
-
-    DynamicList<point> snappedPoints(nCutCells);
 
     // Per cc -1 or a point inside snappedPoints.
     labelList snappedCc;
@@ -1619,7 +1463,6 @@ Foam::isoSurfaceCell::isoSurfaceCell
             isTet,
             cVals,
             pVals,
-            cellCutType,
             snappedPoints,
             snappedCc
         );
@@ -1630,6 +1473,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
         snappedCc = -1;
     }
 
+    snappedPoints.shrink();
+
     if (debug)
     {
         Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
@@ -1648,7 +1493,6 @@ Foam::isoSurfaceCell::isoSurfaceCell
             isTet,
             cVals,
             pVals,
-            cellCutType,
             snappedPoints,
             snappedPoint
         );
@@ -1667,115 +1511,24 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
 
 
-    DynamicList<point> triPoints(nCutCells);
+    DynamicList<point> triPoints(nCutCells_);
-    DynamicList<label> triMeshCells(nCutCells);
+    DynamicList<label> triMeshCells(nCutCells_);
 
-    forAll(mesh_.cells(), cellI)
+    generateTriPoints
-    {
+    (
-        if (cellCutType[cellI] != NOTCUT)
+        cVals,
-        {
+        pVals,
-            label oldNPoints = triPoints.size();
 
-            const cell& cFaces = mesh_.cells()[cellI];
+        mesh_.cellCentres(),
+        mesh_.points(),
 
-            if (isTet.get(cellI) == 1)
+        snappedPoints,
-            {
+        snappedCc,
-                // For tets don't do cell-centre decomposition, just use the
+        snappedPoint,
-                // tet points and values
 
-                const face& f0 = mesh_.faces()[cFaces[0]];
+        triPoints,
-
+        triMeshCells
-                // Get the other point
+    );
-                const face& f1 = mesh_.faces()[cFaces[1]];
-                label oppositeI = -1;
-                forAll(f1, fp)
-                {
-                    oppositeI = f1[fp];
-
-                    if (findIndex(f0, oppositeI) == -1)
-                    {
-                        break;
-                    }
-                }
-
-                generateTriPoints
-                (
-                    snappedPoints,
-                    pVals[f0[0]],
-                    mesh_.points()[f0[0]],
-                    snappedPoint[f0[0]],
-
-                    pVals[f0[1]],
-                    mesh_.points()[f0[1]],
-                    snappedPoint[f0[1]],
-
-                    pVals[f0[2]],
-                    mesh_.points()[f0[2]],
-                    snappedPoint[f0[2]],
-
-                    pVals[oppositeI],
-                    mesh_.points()[oppositeI],
-                    snappedPoint[oppositeI],
-
-                    triPoints
-                );
-            }
-            else
-            {
-                const cell& cFaces = mesh_.cells()[cellI];
-
-                forAll(cFaces, cFaceI)
-                {
-                    label faceI = cFaces[cFaceI];
-                    const face& f = mesh_.faces()[faceI];
-
-                    for (label fp = 1; fp < f.size() - 1; fp++)
-                    {
-                        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-                    //List<triFace> tris(triangulate(f));
-                    //forAll(tris, i)
-                    //{
-                    //    const triFace& tri = tris[i];
-
-
-                        generateTriPoints
-                        (
-                            snappedPoints,
-
-                            pVals[tri[0]],
-                            mesh_.points()[tri[0]],
-                            snappedPoint[tri[0]],
-
-                            pVals[tri[1]],
-                            mesh_.points()[tri[1]],
-                            snappedPoint[tri[1]],
-
-                            pVals[tri[2]],
-                            mesh_.points()[tri[2]],
-                            snappedPoint[tri[2]],
-
-                            cVals[cellI],
-                            mesh_.cellCentres()[cellI],
-                            snappedCc[cellI],
-
-                            triPoints
-                        );
-                    }
-                }
-            }
-
-
-            // Every three triPoints is a cell
-            label nCells = (triPoints.size()-oldNPoints)/3;
-            for (label i = 0; i < nCells; i++)
-            {
-                triMeshCells.append(cellI);
-            }
-        }
-    }
-
-    triPoints.shrink();
-    triMeshCells.shrink();
 
     if (debug)
     {
@@ -1878,123 +1631,123 @@ Foam::isoSurfaceCell::isoSurfaceCell
 }
 
 
-//XXXXXXX
+////XXXXXXX
-// Experimental retriangulation of triangles per cell. Problem is that
+//// Experimental retriangulation of triangles per cell. Problem is that
-// -it is very expensive   -only gets rid of internal points, not of boundary
+//// -it is very expensive   -only gets rid of internal points, not of boundary
-// ones so limited benefit (e.g. 60 v.s. 88 triangles)
+//// ones so limited benefit (e.g. 60 v.s. 88 triangles)
-void Foam::isoSurfaceCell::combineCellTriangles()
+//void Foam::isoSurfaceCell::combineCellTriangles()
-{
+//{
-    if (size() > 0)
+//    if (size() > 0)
-    {
+//    {
-        DynamicList<labelledTri> newTris(size());
+//        DynamicList<labelledTri> newTris(size());
-        DynamicList<label> newTriToCell(size());
+//        DynamicList<label> newTriToCell(size());
-
+//
-        label startTriI = 0;
+//        label startTriI = 0;
-
+//
-        DynamicList<labelledTri> tris;
+//        DynamicList<labelledTri> tris;
-
+//
-        for (label triI = 1; triI <= meshCells_.size(); triI++)
+//        for (label triI = 1; triI <= meshCells_.size(); triI++)
-        {
+//        {
-            if
+//            if
-            (
+//            (
-                triI == meshCells_.size()
+//                triI == meshCells_.size()
-             || meshCells_[triI] != meshCells_[startTriI]
+//             || meshCells_[triI] != meshCells_[startTriI]
-            )
+//            )
-            {
+//            {
-                label nTris = triI-startTriI;
+//                label nTris = triI-startTriI;
-
+//
-                if (nTris == 1)
+//                if (nTris == 1)
-                {
+//                {
-                    newTris.append(operator[](startTriI));
+//                    newTris.append(operator[](startTriI));
-                    newTriToCell.append(meshCells_[startTriI]);
+//                    newTriToCell.append(meshCells_[startTriI]);
-                }
+//                }
-                else
+//                else
-                {
+//                {
-                    // Collect from startTriI to triI in a triSurface
+//                    // Collect from startTriI to triI in a triSurface
-                    tris.clear();
+//                    tris.clear();
-                    for (label i = startTriI; i < triI; i++)
+//                    for (label i = startTriI; i < triI; i++)
-                    {
+//                    {
-                        tris.append(operator[](i));
+//                        tris.append(operator[](i));
-                    }
+//                    }
-                    triSurface cellTris(tris, patches(), points());
+//                    triSurface cellTris(tris, patches(), points());
-                    tris.clear();
+//                    tris.clear();
-
+//
-                    // Get outside
+//                    // Get outside
-                    const labelListList& loops = cellTris.edgeLoops();
+//                    const labelListList& loops = cellTris.edgeLoops();
-
+//
-                    forAll(loops, i)
+//                    forAll(loops, i)
-                    {
+//                    {
-                        // Do proper triangulation of loop
+//                        // Do proper triangulation of loop
-                        face loop(renumber(cellTris.meshPoints(), loops[i]));
+//                        face loop(renumber(cellTris.meshPoints(), loops[i]));
-
+//
-                        faceTriangulation faceTris
+//                        faceTriangulation faceTris
-                        (
+//                        (
-                            points(),
+//                            points(),
-                            loop,
+//                            loop,
-                            true
+//                            true
-                        );
+//                        );
-
+//
-                        // Copy into newTris
+//                        // Copy into newTris
-                        forAll(faceTris, faceTriI)
+//                        forAll(faceTris, faceTriI)
-                        {
+//                        {
-                            const triFace& tri = faceTris[faceTriI];
+//                            const triFace& tri = faceTris[faceTriI];
-
+//
-                            newTris.append
+//                            newTris.append
-                            (
+//                            (
-                                labelledTri
+//                                labelledTri
-                                (
+//                                (
-                                    tri[0],
+//                                    tri[0],
-                                    tri[1],
+//                                    tri[1],
-                                    tri[2],
+//                                    tri[2],
-                                    operator[](startTriI).region()
+//                                    operator[](startTriI).region()
-                                )
+//                                )
-                            );
+//                            );
-                            newTriToCell.append(meshCells_[startTriI]);
+//                            newTriToCell.append(meshCells_[startTriI]);
-                        }
+//                        }
-                    }
+//                    }
-                }
+//                }
-
+//
-                startTriI = triI;
+//                startTriI = triI;
-            }
+//            }
-        }
+//        }
-        newTris.shrink();
+//        newTris.shrink();
-        newTriToCell.shrink();
+//        newTriToCell.shrink();
-
+//
-        // Compact
+//        // Compact
-        pointField newPoints(points().size());
+//        pointField newPoints(points().size());
-        label newPointI = 0;
+//        label newPointI = 0;
-        labelList oldToNewPoint(points().size(), -1);
+//        labelList oldToNewPoint(points().size(), -1);
-
+//
-        forAll(newTris, i)
+//        forAll(newTris, i)
-        {
+//        {
-            labelledTri& tri = newTris[i];
+//            labelledTri& tri = newTris[i];
-            forAll(tri, j)
+//            forAll(tri, j)
-            {
+//            {
-                label pointI = tri[j];
+//                label pointI = tri[j];
-
+//
-                if (oldToNewPoint[pointI] == -1)
+//                if (oldToNewPoint[pointI] == -1)
-                {
+//                {
-                    oldToNewPoint[pointI] = newPointI;
+//                    oldToNewPoint[pointI] = newPointI;
-                    newPoints[newPointI++] = points()[pointI];
+//                    newPoints[newPointI++] = points()[pointI];
-                }
+//                }
-                tri[j] = oldToNewPoint[pointI];
+//                tri[j] = oldToNewPoint[pointI];
-            }
+//            }
-        }
+//        }
-        newPoints.setSize(newPointI);
+//        newPoints.setSize(newPointI);
-
+//
-        triSurface::operator=
+//        triSurface::operator=
-        (
+//        (
-            triSurface
+//            triSurface
-            (
+//            (
-                newTris,
+//                newTris,
-                patches(),
+//                patches(),
-                newPoints,
+//                newPoints,
-                true
+//                true
-            )
+//            )
-        );
+//        );
-        meshCells_.transfer(newTriToCell);
+//        meshCells_.transfer(newTriToCell);
-    }
+//    }
-}
+//}
-//XXXXXXX
+////XXXXXXX
 
 // ************************************************************************* //

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

@@ -91,6 +91,11 @@ class isoSurfaceCell
         //- When to merge points
         const scalar mergeDistance_;
 
+        //- Whether cell might be cut
+        List<cellCutType> cellCutType_;
+
+        //- Estimated number of cut cells
+        label nCutCells_;
 
         //- For every triangle the original cell in mesh
         labelList meshCells_;
@@ -119,6 +124,13 @@ class isoSurfaceCell
             const label
         ) const;
 
+        void calcCutTypes
+        (
+            const PackedList<1>&,
+            const scalarField& cellValues,
+            const scalarField& pointValues
+        );
+
         static labelPair findCommonPoints
         (
             const labelledTri&,
@@ -140,7 +152,6 @@ class isoSurfaceCell
             const PackedList<1>& isTet,
             const scalarField& cVals,
             const scalarField& pVals,
-            const List<cellCutType>& cellCutType,
             DynamicList<point>& triPoints,
             labelList& snappedCc
         ) const;
@@ -173,39 +184,57 @@ class isoSurfaceCell
             const PackedList<1>& isTet,
             const scalarField& cVals,
             const scalarField& pVals,
-            const List<cellCutType>& cellCutType,
             DynamicList<point>& triPoints,
             labelList& snappedPoint
         ) const;
 
         //- Generate single point by interpolation or snapping
-        point generatePoint
+        template<class Type>
+        Type generatePoint
         (
-            const DynamicList<point>& snappedPoints,
+            const DynamicList<Type>& snappedPoints,
             const scalar s0,
-            const point& p0,
+            const Type& p0,
             const label p0Index,
             const scalar s1,
-            const point& p1,
+            const Type& p1,
             const label p1Index
         ) const;
 
+        template<class Type>
         void generateTriPoints
         (
-            const DynamicList<point>& snapped,
+            const DynamicList<Type>& snapped,
             const scalar s0,
-            const point& p0,
+            const Type& p0,
             const label p0Index,
             const scalar s1,
-            const point& p1,
+            const Type& p1,
             const label p1Index,
             const scalar s2,
-            const point& p2,
+            const Type& p2,
             const label p2Index,
             const scalar s3,
-            const point& p3,
+            const Type& p3,
             const label p3Index,
-            DynamicList<point>& points
+            DynamicList<Type>& points
+        ) const;
+
+        template<class Type>
+        void generateTriPoints
+        (
+            const scalarField& cVals,
+            const scalarField& pVals,
+
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords,
+
+            const DynamicList<Type>& snappedPoints,
+            const labelList& snappedCc,
+            const labelList& snappedPoint,
+
+            DynamicList<Type>& triPoints,
+            DynamicList<label>& triMeshCells
         ) const;
 
         triSurface stitchTriPoints
@@ -311,6 +340,18 @@ public:
         {
             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;
 };
 
 
@@ -320,6 +361,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "isoSurfaceCellTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C

@@ -0,0 +1,384 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoSurfaceCell.H"
+#include "polyMesh.H"
+#include "tetMatcher.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::isoSurfaceCell::generatePoint
+(
+    const DynamicList<Type>& snappedPoints,
+
+    const scalar s0,
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index
+) const
+{
+    scalar d = s1-s0;
+
+    if (mag(d) > VSMALL)
+    {
+        scalar s = (iso_-s0)/d;
+
+        if (s >= 0.5 && s <= 1 && p1Index != -1)
+        {
+            return snappedPoints[p1Index];
+        }
+        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
+        {
+            return snappedPoints[p0Index];
+        }
+        else
+        {
+            return s*p1 + (1.0-s)*p0;
+        }
+    }
+    else
+    {
+        scalar s = 0.4999;
+
+        return s*p1 + (1.0-s)*p0;
+    }
+}
+
+
+template<class Type>
+void Foam::isoSurfaceCell::generateTriPoints
+(
+    const DynamicList<Type>& snapped,
+
+    const scalar s0,
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index,
+
+    const scalar s2,
+    const Type& p2,
+    const label p2Index,
+
+    const scalar s3,
+    const Type& p3,
+    const label p3Index,
+
+    DynamicList<Type>& points
+) const
+{
+    int triIndex = 0;
+    if (s0 < iso_)
+    {
+        triIndex |= 1;
+    }
+    if (s1 < iso_)
+    {
+        triIndex |= 2;
+    }
+    if (s2 < iso_)
+    {
+        triIndex |= 4;
+    }
+    if (s3 < iso_)
+    {
+        triIndex |= 8;
+    }
+
+    /* Form the vertices of the triangles for each case */
+    switch (triIndex)
+    {
+        case 0x00:
+        case 0x0F:
+        break;
+
+        case 0x0E:
+        case 0x01:
+            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+        break;
+
+        case 0x0D:
+        case 0x02:
+            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+        break;
+
+        case 0x0C:
+        case 0x03:
+        {
+            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp2);
+            points.append(tp2);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x0B:
+        case 0x04:
+        {
+            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+        }
+        break;
+
+        case 0x0A:
+        case 0x05:
+        {
+            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(tp1);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x09:
+        case 0x06:
+        {
+            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp0);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x07:
+        case 0x08:
+            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+        break;
+    }
+}
+
+
+template<class Type>
+void Foam::isoSurfaceCell::generateTriPoints
+(
+    const scalarField& cVals,
+    const scalarField& pVals,
+
+    const Field<Type>& cCoords,
+    const Field<Type>& pCoords,
+
+    const DynamicList<Type>& snappedPoints,
+    const labelList& snappedCc,
+    const labelList& snappedPoint,
+
+    DynamicList<Type>& triPoints,
+    DynamicList<label>& triMeshCells
+) const
+{
+    tetMatcher tet;
+
+    forAll(mesh_.cells(), cellI)
+    {
+        if (cellCutType_[cellI] != NOTCUT)
+        {
+            label oldNPoints = triPoints.size();
+
+            const cell& cFaces = mesh_.cells()[cellI];
+
+            if (tet.isA(mesh_, cellI))
+            {
+                // For tets don't do cell-centre decomposition, just use the
+                // tet points and values
+
+                const face& f0 = mesh_.faces()[cFaces[0]];
+
+                // Get the other point
+                const face& f1 = mesh_.faces()[cFaces[1]];
+                label oppositeI = -1;
+                forAll(f1, fp)
+                {
+                    oppositeI = f1[fp];
+
+                    if (findIndex(f0, oppositeI) == -1)
+                    {
+                        break;
+                    }
+                }
+
+                generateTriPoints
+                (
+                    snappedPoints,
+                    pVals[f0[0]],
+                    pCoords[f0[0]],
+                    snappedPoint[f0[0]],
+
+                    pVals[f0[1]],
+                    pCoords[f0[1]],
+                    snappedPoint[f0[1]],
+
+                    pVals[f0[2]],
+                    pCoords[f0[2]],
+                    snappedPoint[f0[2]],
+
+                    pVals[oppositeI],
+                    pCoords[oppositeI],
+                    snappedPoint[oppositeI],
+
+                    triPoints
+                );
+            }
+            else
+            {
+                const cell& cFaces = mesh_.cells()[cellI];
+
+                forAll(cFaces, cFaceI)
+                {
+                    label faceI = cFaces[cFaceI];
+                    const face& f = mesh_.faces()[faceI];
+
+                    for (label fp = 1; fp < f.size() - 1; fp++)
+                    {
+                        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+                    //List<triFace> tris(triangulate(f));
+                    //forAll(tris, i)
+                    //{
+                    //    const triFace& tri = tris[i];
+
+
+                        generateTriPoints
+                        (
+                            snappedPoints,
+
+                            pVals[tri[0]],
+                            pCoords[tri[0]],
+                            snappedPoint[tri[0]],
+
+                            pVals[tri[1]],
+                            pCoords[tri[1]],
+                            snappedPoint[tri[1]],
+
+                            pVals[tri[2]],
+                            pCoords[tri[2]],
+                            snappedPoint[tri[2]],
+
+                            cVals[cellI],
+                            cCoords[cellI],
+                            snappedCc[cellI],
+
+                            triPoints
+                        );
+                    }
+                }
+            }
+
+
+            // Every three triPoints is a cell
+            label nCells = (triPoints.size()-oldNPoints)/3;
+            for (label i = 0; i < nCells; i++)
+            {
+                triMeshCells.append(cellI);
+            }
+        }
+    }
+
+    triPoints.shrink();
+    triMeshCells.shrink();
+}
+
+
+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<label> triMeshCells(nCutCells_);
+
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    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;
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -382,6 +382,39 @@ void Foam::isoSurface::generateTriPoints
                 }
             }
         }
+        else if (isA<emptyPolyPatch>(pp))
+        {
+            // Assume zero-gradient.
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                if (faceCutType_[faceI] != NOTCUT)
+                {
+                    generateTriPoints
+                    (
+                        cVals,
+                        pVals,
+
+                        cCoords,
+                        pCoords,
+
+                        snappedPoints,
+                        snappedCc,
+                        snappedPoint,
+                        faceI,
+
+                        cVals[own[faceI]],
+                        cCoords.boundaryField()[patchI][i],
+                        -1, // fc not snapped
+
+                        triPoints,
+                        triMeshCells
+                    );
+                }
+                faceI++;
+            }
+        }
         else
         {
             label faceI = pp.start();
@@ -465,8 +498,12 @@ Foam::isoSurface::interpolate
 
 
     // One value per point
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    tmp<Field<Type> > tvalues
+    (
+        new Field<Type>(points().size(), pTraits<Type>::zero)
+    );
     Field<Type>& values = tvalues();
+    labelList nValues(values.size(), 0);
 
     forAll(triPoints, i)
     {
@@ -474,10 +511,36 @@ Foam::isoSurface::interpolate
 
         if (mergedPointI >= 0)
         {
-            values[mergedPointI] = triPoints[i];
+            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

@@ -65,7 +65,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
         if (debug)
         {
             Info<< "sampledIsoSurface::getIsoField() : reading "
-                << isoField_ << " from time " <<fvm.time().timeName()
+                << isoField_ << " from time " << fvm.time().timeName()
                 << endl;
         }
 
@@ -124,10 +124,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
 
     if (debug)
     {
-        Info<< "sampledIsoSurface::getIsoField() : obtained volField "
+        Info<< "sampledIsoSurface::getIsoField() : volField "
             << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
             << " max:" << max(*volFieldPtr_).value() << endl;
-        Info<< "sampledIsoSurface::getIsoField() : obtained pointField "
+        Info<< "sampledIsoSurface::getIsoField() : pointField "
             << pointFieldPtr_->name()
             << " min:" << gMin(pointFieldPtr_->internalField())
             << " max:" << gMax(pointFieldPtr_->internalField()) << endl;

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

@@ -57,11 +57,16 @@ Foam::sampledIsoSurface::interpolateField
     const GeometricField<Type, fvPatchField, volMesh>& volFld =
         interpolator.psi();
 
-    tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld
+    tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointFld
     (
         volPointInterpolation::New(fvm).interpolate(volFld)
     );
 
+    const GeometricField<Type, pointPatchField, pointMesh>& pointFld =
+        tpointFld();
+
+    // Get pointers to sampling field (both original and interpolated one)
+    getIsoFields();
     // Recreate geometry if time has changed
     createGeometry();
 
@@ -71,7 +76,7 @@ Foam::sampledIsoSurface::interpolateField
         *volFieldPtr_,
         *pointFieldPtr_,
         volFld,
-        pointFld()
+        pointFld
     );
 }
 

src/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.C

@@ -49,7 +49,7 @@ Foam::laplaceFilter::laplaceFilter(const fvMesh& mesh, scalar widthCoeff)
     (
         IOobject
         (
-            "anisotropicFilterCoeff",
+            "laplaceFilterCoeff",
             mesh.time().timeName(),
             mesh
         ),
@@ -70,7 +70,7 @@ Foam::laplaceFilter::laplaceFilter(const fvMesh& mesh, const dictionary& bd)
     (
         IOobject
         (
-            "anisotropicFilterCoeff",
+            "laplaceFilterCoeff",
             mesh.time().timeName(),
             mesh
         ),

src/turbulenceModels/compressible/LES/LESModel/LESModel.C

@@ -38,6 +38,7 @@ namespace compressible
 
 defineTypeNameAndDebug(LESModel, 0);
 defineRunTimeSelectionTable(LESModel, dictionary);
+addToRunTimeSelectionTable(turbulenceModel, LESModel, turbulenceModel);
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 

src/turbulenceModels/compressible/RAS/RASModel/RASModel.C

@@ -39,6 +39,7 @@ namespace compressible
 
 defineTypeNameAndDebug(RASModel, 0);
 defineRunTimeSelectionTable(RASModel, dictionary);
+addToRunTimeSelectionTable(turbulenceModel, RASModel, turbulenceModel);
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //