ENH: tet decomposed particle tracking.

Squashed merge of particleInteractions up to
commit e7cb5bcf0315c359539ef1e715e1d51991343391

applications/solvers/combustion/dieselEngineFoam/Make/options

@@ -2,6 +2,7 @@ EXE_INC = \
     -I../engineFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/dieselSpray/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
@@ -24,6 +25,7 @@ EXE_LIBS = \
     -lreactionThermophysicalModels \
     -lfiniteVolume \
     -llagrangian \
+    -lmeshTools \
     -ldieselSpray \
     -lliquids \
     -lliquidMixture \

applications/solvers/combustion/dieselFoam/Make/options

@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/dieselSpray/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
@@ -21,6 +22,7 @@ EXE_LIBS = \
     -lcompressibleLESModels \
     -lreactionThermophysicalModels \
     -llagrangian \
+    -lmeshTools \
     -ldieselSpray \
     -lliquids \
     -lliquidMixture \

applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options

@@ -9,4 +9,3 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -ldsmc
-

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/Make/files

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

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/Make/options

@@ -0,0 +1,33 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+
+
+EXE_LIBS = \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lthermophysicalFunctions \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lradiation \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsurfaceFilmModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/createFields.H

@@ -0,0 +1,147 @@
+    Info<< "\nReading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar rhoInfValue
+    (
+        transportProperties.lookup("rhoInf")
+    );
+
+    volScalarField rhoInf
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        rhoInfValue
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    const volScalarField nu = laminarTransport.nu();
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+    volScalarField mu
+    (
+        IOobject
+        (
+            "mu",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        nu*rhoInfValue
+    );
+
+    word kinematicCloudName("kinematicCloud");
+    args.optionReadIfPresent("cloudName", kinematicCloudName);
+
+    Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
+    basicKinematicCloud kinematicCloud
+    (
+        kinematicCloudName,
+        rhoInf,
+        U,
+        mu,
+        g
+    );
+
+    IOobject Hheader
+    (
+        "H",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ
+    );
+
+    autoPtr<volVectorField> HPtr_;
+
+    if (Hheader.headerOk())
+    {
+        Info<< "\nReading field H\n" << endl;
+
+        HPtr_.reset
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    "H",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh
+            )
+        );
+    }
+
+    IOobject HdotGradHheader
+    (
+        "HdotGradH",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ
+    );
+
+    autoPtr<volVectorField> HdotGradHPtr_;
+
+    if (HdotGradHheader.headerOk())
+    {
+        Info<< "Reading field HdotGradH" << endl;
+
+        HdotGradHPtr_.reset
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    "HdotGradH",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh
+            )
+        );
+    }

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    uncoupledKinematicParcelFoam
+
+Description
+    Transient solver for the passive transport of a single kinematic
+    particle could.
+
+    Uses a pre-calculated velocity field to evolve the cloud.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+#include "basicKinematicCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "cloudName",
+        "name",
+        "specify alternative cloud name. default is 'kinematicCloud'"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #   include "createDynamicFvMesh.H"
+
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        mesh.update();
+
+        U.correctBoundaryConditions();
+
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+
+        laminarTransport.correct();
+
+        mu = nu*rhoInfValue;
+
+        kinematicCloud.evolve();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/findCell-octree/findCell-octree.C

@@ -29,6 +29,8 @@ License
 #include "IStringStream.H"
 #include "octree.H"
 #include "octreeDataCell.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
 #include "OFstream.H"
 
 using namespace Foam;
@@ -44,11 +46,13 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
+    label nReps = 100000;
+
     const point sample = args.argRead<point>(1);
 
-    treeBoundBox meshBb(mesh.points());
+    treeBoundBox meshBb(mesh.bounds());
 
-    // Calculate typical cell releated size to shift bb by.
+    // Calculate typical cell related size to shift bb by.
     scalar typDim = meshBb.avgDim()/(2.0*Foam::cbrt(scalar(mesh.nCells())));
 
     treeBoundBox shiftedBb
@@ -57,16 +61,13 @@ int main(int argc, char *argv[])
         meshBb.max() + vector(typDim, typDim, typDim)
     );
 
-
     Info<< "Mesh" << endl;
     Info<< "   bounding box           : " << meshBb << endl;
     Info<< "   bounding box (shifted) : " << shiftedBb << endl;
     Info<< "   typical dimension      : " << shiftedBb.typDim() << endl;
 
-
-    /*
-     * Now we have allBb and shiftedBb
-     */
+    Info<< "Initialised mesh in "
+        << runTime.cpuTimeIncrement() << " s" << endl;
 
     // Wrap indices and mesh information into helper object
     octreeDataCell shapes(mesh);
@@ -80,14 +81,52 @@ int main(int argc, char *argv[])
         10.0        // maximum ratio of cubes v.s. cells
     );
 
-    Info<< "Point:" << sample << " is in shape "
-        << oc.find(sample) << nl
-        << "Point:" << sample << " is in cell  "
-        << mesh.findCell(sample) << endl;
+    for (label i = 0; i < nReps - 1 ; i++)
+    {
+        oc.find(sample);
+    }
 
+    Info<< "Point:" << sample << " is in shape "
+        << oc.find(sample) << endl;
 
     oc.printStats(Info);
 
+    Info<< "Found in octree " << nReps << " times in "
+        << runTime.cpuTimeIncrement() << " s" << endl;
+
+    indexedOctree<treeDataCell> ioc
+    (
+        treeDataCell(true, mesh),
+        shiftedBb,
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+    for (label i = 0; i < nReps - 1 ; i++)
+    {
+        ioc.findInside(sample);
+    }
+
+    Info<< "Point:" << sample << " is in shape "
+        << ioc.findInside(sample)
+        << ", where the possible cells were:" << nl
+        << ioc.findIndices(sample)
+        << endl;
+
+    Info<< "Found in indexedOctree " << nReps << " times in "
+        << runTime.cpuTimeIncrement() << " s" << endl;
+
+    for (label i = 0; i < nReps - 1 ; i++)
+    {
+        mesh.findCell(sample);
+    }
+
+    Info<< "Point:" << sample << " is in cell  "
+        << mesh.findCell(sample) << endl;
+
+    Info<< "Found in mesh.findCell " << nReps << " times in "
+        << runTime.cpuTimeIncrement() << " s" << endl;
 
     Info<< "End\n" << endl;
 

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -341,11 +341,14 @@ meshQualityControls
     //  Set to very negative number (e.g. -1E30) to disable.
     minVol 1e-13;
 
-    //- Minimum tet volume. Is absolute volume of the tet formed by the
-    //  face-centre decomposition triangle and the cell centre.
-    //  Set to a sensible fraction of the smallest cell volume expected.
-    //  Set to very negative number (e.g. -1E30) to disable.
-    minTetVol 1e-20;
+    //- Minimum quality of the tet formed by the face-centre
+    //  and variable base point minimum decomposition triangles and
+    //  the cell centre.  Set to very negative number (e.g. -1E30) to
+    //  disable.
+    //     <0 = inside out tet,
+    //      0 = flat tet
+    //      1 = regular tet
+    minTetQuality 1e-9;
 
     //- Minimum face area. Set to <0 to disable.
     minArea -1;

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

@@ -6,6 +6,7 @@
 #include "EdgeMap.H"
 #include "wedgePolyPatch.H"
 #include "unitConversion.H"
+#include "polyMeshTetDecomposition.H"
 
 
 // Find wedge with opposite orientation. Note: does not actually check that
@@ -413,7 +414,6 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
-
     {
         faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
         if (mesh.checkFacePyramids(true, -SMALL, &faces))
@@ -433,25 +433,6 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
-    {
-        faceSet faces(mesh, "wrongOrientedTriangleFaces", mesh.nFaces()/100+1);
-        if (mesh.checkFaceTets(true, 0, &faces))
-        {
-            noFailedChecks++;
-
-            label nFaces = returnReduce(faces.size(), sumOp<label>());
-
-            if (nFaces > 0)
-            {
-                Info<< "  <<Writing " << nFaces
-                    << " faces with incorrectly orientated triangles to set "
-                    << faces.name() << endl;
-                faces.instance() = mesh.pointsInstance();
-                faces.write();
-            }
-        }
-    }
-
     {
         faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
         if (mesh.checkFaceSkewness(true, &faces))
@@ -470,6 +451,35 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
+    if (allGeometry)
+    {
+        faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
+        if
+        (
+            polyMeshTetDecomposition::checkFaceTets
+            (
+                mesh,
+                polyMeshTetDecomposition::minTetQuality,
+                true,
+                &faces
+            )
+        )
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
+            {
+                Info<< "  <<Writing " << nFaces
+                    << " faces with low quality or negative volume "
+                    << "decomposition tets to set " << faces.name() << endl;
+                faces.instance() = mesh.pointsInstance();
+                faces.write();
+            }
+        }
+    }
+
     if (allGeometry)
     {
         // Note use of nPoints since don't want edge construction.

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -558,7 +558,6 @@ int main(int argc, char *argv[])
                 lagrangianScalarFieldFields
             );
 
-
             lagrangianFieldDecomposer::readFields
             (
                 cloudI,
@@ -687,6 +686,19 @@ int main(int argc, char *argv[])
             )
         );
 
+        labelIOList faceProcAddressing
+        (
+            IOobject
+            (
+                "faceProcAddressing",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
         labelIOList cellProcAddressing
         (
             IOobject
@@ -728,19 +740,6 @@ int main(int argc, char *argv[])
          || surfaceTensorFields.size()
         )
         {
-            labelIOList faceProcAddressing
-            (
-                IOobject
-                (
-                    "faceProcAddressing",
-                    procMesh.facesInstance(),
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                )
-            );
-
             fvFieldDecomposer fieldDecomposer
             (
                 mesh,
@@ -814,6 +813,7 @@ int main(int argc, char *argv[])
                 (
                     mesh,
                     procMesh,
+                    faceProcAddressing,
                     cellProcAddressing,
                     cloudDirs[cloudI],
                     lagrangianPositions[cloudI],

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C

@@ -35,6 +35,7 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
 (
     const polyMesh& mesh,
     const polyMesh& procMesh,
+    const labelList& faceProcAddressing,
     const labelList& cellProcAddressing,
     const word& cloudName,
     const Cloud<indexedParticle>& lagrangianPositions,
@@ -47,6 +48,14 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
 {
     label pi = 0;
 
+    // faceProcAddressing not required currently
+    // labelList decodedProcFaceAddressing(faceProcAddressing.size());
+
+    // forAll(faceProcAddressing, i)
+    // {
+    //     decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
+    // }
+
     forAll(cellProcAddressing, procCelli)
     {
         label celli = cellProcAddressing[procCelli];
@@ -60,6 +69,33 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
                 const indexedParticle& ppi = *iter();
                 particleIndices_[pi++] = ppi.index();
 
+                // label mappedTetFace = findIndex
+                // (
+                //     decodedProcFaceAddressing,
+                //     ppi.tetFace()
+                // );
+
+                // if (mappedTetFace == -1)
+                // {
+                //     FatalErrorIn
+                //     (
+                //         "Foam::lagrangianFieldDecomposer"
+                //         "::lagrangianFieldDecomposer"
+                //         "("
+                //             "const polyMesh& mesh, "
+                //             "const polyMesh& procMesh, "
+                //             "const labelList& faceProcAddressing, "
+                //             "const labelList& cellProcAddressing, "
+                //             "const word& cloudName, "
+                //             "const Cloud<indexedParticle>& "
+                //             "lagrangianPositions, "
+                //             "const List<SLList<indexedParticle*>*>& "
+                //             "cellParticles"
+                //         ")"
+                //     )   << "Face lookup failure." << nl
+                //         << abort(FatalError);
+                // }
+
                 positions_.append
                 (
                     new passiveParticle

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H

@@ -84,6 +84,7 @@ public:
         (
             const polyMesh& mesh,
             const polyMesh& procMesh,
+            const labelList& faceProcAddressing,
             const labelList& cellProcAddressing,
             const word& cloudName,
             const Cloud<indexedParticle>& lagrangianPositions,

applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
     /* -DFULLDEBUG -g -O0 */ \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
     -lgenericPatchFields \
     -llagrangian
 

applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/conversion/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
+    -lmeshTools \
     -lgenericPatchFields \
     -lconversion

applications/utilities/postProcessing/dataConversion/foamToFieldview9/Make/options

@@ -1,8 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
     -lgenericPatchFields \
     -llagrangian

applications/utilities/postProcessing/dataConversion/foamToGMV/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/lagrangian/lnInclude \
@@ -14,5 +15,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
     -lgenericPatchFields \
     -llagrangian

applications/utilities/postProcessing/graphics/ensightFoamReader/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/browser/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
@@ -7,6 +8,7 @@ EXE_INC = \
 LIB_LIBS = \
     -lOpenFOAM \
     -lfiniteVolume \
+    -lmeshTools \
     -lgenericPatchFields \
     -llagrangian \
     $(PROJECT_LIBS)

applications/utilities/postProcessing/lagrangian/particleTracks/Make/options

@@ -1,8 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
     -lgenericPatchFields \
     -llagrangian

applications/utilities/preProcessing/dsmcInitialise/Make/options

@@ -9,4 +9,3 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -ldsmc
-

applications/utilities/preProcessing/mapFields/MapLagrangianFields.H

@@ -37,6 +37,7 @@ Description
 #include "GeometricField.H"
 #include "meshToMesh.H"
 #include "IOobjectList.H"
+#include "IOFieldField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,11 +57,13 @@ void MapLagrangianFields
 {
     const fvMesh& meshTarget = meshToMeshInterp.toMesh();
 
+    {
         IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
 
         forAllIter(IOobjectList, fields, fieldIter)
         {
-        Info<< "    mapping lagrangian field " << fieldIter()->name() << endl;
+            Info<< "    mapping lagrangian field "
+                << fieldIter()->name() << endl;
 
             // Read field (does not need mesh)
             IOField<Type> fieldSource(*fieldIter());
@@ -80,6 +83,7 @@ void MapLagrangianFields
                 ),
                 addParticles.size()
             );
+
             forAll(addParticles, i)
             {
                 fieldTarget[i] = fieldSource[addParticles[i]];
@@ -90,6 +94,84 @@ void MapLagrangianFields
         }
     }
 
+    {
+        IOobjectList fieldFields =
+            objects.lookupClass(IOField<Field<Type> >::typeName);
+
+        forAllIter(IOobjectList, fieldFields, fieldIter)
+        {
+            Info<< "    mapping lagrangian fieldField "
+                << fieldIter()->name() << endl;
+
+            // Read field (does not need mesh)
+            IOField<Field<Type> > fieldSource(*fieldIter());
+
+            // Map - use IOFieldField to automatically write in
+            // compact form for binary format.
+            IOFieldField<Field<Type>, Type> fieldTarget
+            (
+                IOobject
+                (
+                    fieldIter()->name(),
+                    meshTarget.time().timeName(),
+                    cloud::prefix/cloudName,
+                    meshTarget,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                addParticles.size()
+            );
+
+            forAll(addParticles, i)
+            {
+                fieldTarget[i] = fieldSource[addParticles[i]];
+            }
+
+            // Write field
+            fieldTarget.write();
+        }
+    }
+
+    {
+        IOobjectList fieldFields =
+            objects.lookupClass(IOFieldField<Field<Type>, Type>::typeName);
+
+        forAllIter(IOobjectList, fieldFields, fieldIter)
+        {
+            Info<< "    mapping lagrangian fieldField "
+                << fieldIter()->name() << endl;
+
+            // Read field (does not need mesh)
+            IOFieldField<Field<Type>, Type> fieldSource(*fieldIter());
+
+            // Map
+            IOFieldField<Field<Type>, Type> fieldTarget
+            (
+                IOobject
+                (
+                    fieldIter()->name(),
+                    meshTarget.time().timeName(),
+                    cloud::prefix/cloudName,
+                    meshTarget,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                addParticles.size()
+            );
+
+            forAll(addParticles, i)
+            {
+                fieldTarget[i] = fieldSource[addParticles[i]];
+            }
+
+            // Write field
+            fieldTarget.write();
+        }
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/utilities/preProcessing/mapFields/mapLagrangian.C

@@ -38,12 +38,13 @@ static const scalar perturbFactor = 1E-6;
 
 // Special version of findCell that generates a cell guaranteed to be
 // compatible with tracking.
-static label findCell(const meshSearch& meshSearcher, const point& pt)
+static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
 {
-    const polyMesh& mesh = meshSearcher.mesh();
+    label cellI = -1;
+    label tetFaceI = -1;
+    label tetPtI = -1;
 
-    // Use tracking to find cell containing pt
-    label cellI = meshSearcher.findCell(pt);
+    cloud.findCellFacePt(pt, cellI, tetFaceI, tetPtI);
 
     if (cellI >= 0)
     {
@@ -54,16 +55,24 @@ static label findCell(const meshSearch& meshSearcher, const point& pt)
         // See if particle on face by finding nearest face and shifting
         // particle.
 
+        const polyMesh& mesh = cloud.pMesh();
+
+        meshSearch meshSearcher(mesh, false);
+
         label faceI = meshSearcher.findNearestBoundaryFace(pt);
 
         if (faceI >= 0)
         {
             const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
+
             const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
 
-            return meshSearcher.findCell(perturbPt);
+            cloud.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI);
+
+            return cellI;
         }
     }
+
     return -1;
 }
 
@@ -207,8 +216,6 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp)
 
             if (unmappedSource.size())
             {
-                meshSearch targetSearcher(meshTarget, false);
-
                 sourceParticleI = 0;
 
                 forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
@@ -216,7 +223,7 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp)
                     if (unmappedSource.found(sourceParticleI))
                     {
                         label targetCell =
-                            findCell(targetSearcher, iter().position());
+                            findCell(targetParcels, iter().position());
 
                         if (targetCell >= 0)
                         {

applications/utilities/surface/surfaceSubset/surfaceSubset.C

@@ -35,7 +35,6 @@ Description
 #include "IOdictionary.H"
 #include "boundBox.H"
 #include "indexedOctree.H"
-#include "octree.H"
 #include "treeDataTriSurface.H"
 #include "Random.H"
 

etc/controlDict

@@ -780,6 +780,7 @@ DebugSwitches
     syringePressure     0;
     tensorAverageField  0;
     tensorField         0;
+    tetDecomposedPolyMesh 0;
     thermoCloud         0;
     thermophysicalFunction 0;
     time                0;

src/Allwmake

@@ -20,8 +20,6 @@ Pstream/Allwmake
 OSspecific/$WM_OSTYPE/Allwmake
 wmake libso OpenFOAM
 
-wmake libso lagrangian/basic
-
 wmake libso fileFormats
 wmake libso edgeMesh
 wmake libso surfMesh
@@ -33,6 +31,7 @@ parallel/decompose/AllwmakeLnInclude
 dummyThirdParty/Allwmake
 
 wmake libso meshTools
+wmake libso lagrangian/basic
 wmake libso finiteVolume
 wmake libso genericPatchFields
 

src/OpenFOAM/Make/files

@@ -365,6 +365,8 @@ $(globalMeshData)/globalPoints.C
 $(globalMeshData)/globalIndex.C
 
 $(polyMesh)/syncTools/syncTools.C
+$(polyMesh)/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+$(polyMesh)/polyMeshTetDecomposition/tetIndices.C
 
 zone = $(polyMesh)/zones/zone
 $(zone)/zone.C

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -33,6 +33,7 @@ License
 #include "processorPolyPatch.H"
 #include "OSspecific.H"
 #include "demandDrivenData.H"
+#include "polyMeshTetDecomposition.H"
 
 #include "pointMesh.H"
 
@@ -204,6 +205,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
     bounds_(points_),
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
+    tetBasePtIsPtr_(NULL),
     pointZones_
     (
         IOobject
@@ -392,6 +394,7 @@ Foam::polyMesh::polyMesh
     bounds_(points_, syncPar),
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
+    tetBasePtIsPtr_(NULL),
     pointZones_
     (
         IOobject
@@ -548,6 +551,7 @@ Foam::polyMesh::polyMesh
     bounds_(points_, syncPar),
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
+    tetBasePtIsPtr_(NULL),
     pointZones_
     (
         IOobject
@@ -847,6 +851,28 @@ Foam::label Foam::polyMesh::nSolutionD() const
 }
 
 
+const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
+{
+    if (!tetBasePtIsPtr_)
+    {
+        if (debug)
+        {
+            WarningIn("const labelList& polyMesh::tetBasePtIs() const")
+                << "Tet base point indices not available.  "
+                << "Forcing storage of base points."
+                << endl;
+        }
+
+        tetBasePtIsPtr_ = new labelList
+        (
+            polyMeshTetDecomposition::findFaceBasePts(*this)
+        );
+    }
+
+    return *tetBasePtIsPtr_;
+}
+
+
 // Add boundary patches. Constructor helper
 void Foam::polyMesh::addPatches
 (

src/OpenFOAM/meshes/polyMesh/polyMesh.H

@@ -62,6 +62,7 @@ namespace Foam
 
 class globalMeshData;
 class mapPolyMesh;
+class polyMeshTetDecomposition;
 
 /*---------------------------------------------------------------------------*\
                           Class polyMesh Declaration
@@ -127,6 +128,9 @@ private:
             //  defined according to the presence of empty patches
             mutable Vector<label> solutionD_;
 
+            //- Base point for face decomposition into tets
+            mutable labelList* tetBasePtIsPtr_;
+
 
         // Zoning information
 
@@ -356,6 +360,9 @@ public:
             //- Return the number of valid solved-for dimensions in the mesh
             label nSolutionD() const;
 
+            //- Return the tetBasePtIs
+            const labelList& tetBasePtIs() const;
+
             //- Return point zone mesh
             const pointZoneMesh& pointZones() const
             {

src/OpenFOAM/meshes/polyMesh/polyMeshClear.C

@@ -71,6 +71,9 @@ void Foam::polyMesh::clearGeom()
     geometricD_ = Vector<label>::zero;
     solutionD_ = Vector<label>::zero;
 
+    // Remove the stored tet base points
+    deleteDemandDrivenData(tetBasePtIsPtr_);
+
     pointMesh::Delete(*this);
 }
 

src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C

@@ -503,6 +503,7 @@ Foam::polyMesh::polyMesh
     bounds_(points_, syncPar),
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
+    tetBasePtIsPtr_(NULL),
     pointZones_
     (
         IOobject
@@ -775,6 +776,7 @@ Foam::polyMesh::polyMesh
     bounds_(points_, syncPar),
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
+    tetBasePtIsPtr_(NULL),
     pointZones_
     (
         IOobject

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -0,0 +1,627 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMeshTetDecomposition.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = 1e-9;
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
+(
+    const polyMesh& mesh,
+    label fI,
+    const point& nCc,
+    scalar tol,
+    bool report
+)
+{
+    const faceList& pFaces = mesh.faces();
+    const pointField& pPts = mesh.points();
+    const vectorField& pC = mesh.cellCentres();
+    const labelList& pOwner = mesh.faceOwner();
+
+    const face& f = pFaces[fI];
+
+    label oCI = pOwner[fI];
+
+    const point& oCc = pC[oCI];
+
+    forAll(f, faceBasePtI)
+    {
+        scalar thisBaseMinTetQuality = VGREAT;
+
+        const point& tetBasePt = pPts[f[faceBasePtI]];
+
+        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+        {
+            label facePtI = (tetPtI + faceBasePtI) % f.size();
+            label otherFacePtI = f.fcIndex(facePtI);
+
+            List<scalar> tetQualities(2, 0.0);
+
+            {
+                // owner cell tet
+                label ptAI = f[facePtI];
+                label ptBI = f[otherFacePtI];
+
+                const point& pA = pPts[ptAI];
+                const point& pB = pPts[ptBI];
+
+                tetPointRef tet(oCc, tetBasePt, pA, pB);
+
+                tetQualities[0] = tet.quality();
+            }
+
+            {
+                // neighbour cell tet
+                label ptAI = f[otherFacePtI];
+                label ptBI = f[facePtI];
+
+                const point& pA = pPts[ptAI];
+                const point& pB = pPts[ptBI];
+
+                tetPointRef tet(nCc, tetBasePt, pA, pB);
+
+                tetQualities[1] = tet.quality();
+            }
+
+            if (min(tetQualities) < thisBaseMinTetQuality)
+            {
+                thisBaseMinTetQuality = min(tetQualities);
+            }
+        }
+
+        if (thisBaseMinTetQuality > tol)
+        {
+            return faceBasePtI;
+        }
+
+    }
+
+    // If a base point hasn't triggered a return by now, then there is
+    // non that can produce a good decomposition
+    return -1;
+}
+
+
+Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
+(
+    const polyMesh& mesh,
+    label fI,
+    scalar tol,
+    bool report
+)
+{
+    return findSharedBasePoint
+    (
+        mesh,
+        fI,
+        mesh.cellCentres()[mesh.faceNeighbour()[fI]],
+        tol,
+        report
+    );
+}
+
+
+Foam::label Foam::polyMeshTetDecomposition::findBasePoint
+(
+    const polyMesh& mesh,
+    label fI,
+    scalar tol,
+    bool report
+)
+{
+    const faceList& pFaces = mesh.faces();
+    const pointField& pPts = mesh.points();
+    const vectorField& pC = mesh.cellCentres();
+    const labelList& pOwner = mesh.faceOwner();
+
+    const face& f = pFaces[fI];
+
+    label cI = pOwner[fI];
+
+    bool own = (pOwner[fI] == cI);
+
+    const point& cC = pC[cI];
+
+    forAll(f, faceBasePtI)
+    {
+        scalar thisBaseMinTetQuality = VGREAT;
+
+        const point& tetBasePt = pPts[f[faceBasePtI]];
+
+        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+        {
+            label facePtI = (tetPtI + faceBasePtI) % f.size();
+            label otherFacePtI = f.fcIndex(facePtI);
+
+            label ptAI = -1;
+            label ptBI = -1;
+
+            if (own)
+            {
+                ptAI = f[facePtI];
+                ptBI = f[otherFacePtI];
+            }
+            else
+            {
+                ptAI = f[otherFacePtI];
+                ptBI = f[facePtI];
+            }
+
+            const point& pA = pPts[ptAI];
+            const point& pB = pPts[ptBI];
+
+            tetPointRef tet(cC, tetBasePt, pA, pB);
+
+            scalar tetQuality = tet.quality();
+
+            if (tetQuality < thisBaseMinTetQuality)
+            {
+                thisBaseMinTetQuality = tetQuality;
+            }
+        }
+
+        if (thisBaseMinTetQuality > tol)
+        {
+            return faceBasePtI;
+        }
+    }
+
+    // If a base point hasn't triggered a return by now, then there is
+    // non that can produce a good decomposition
+    return -1;
+}
+
+
+Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
+(
+    const polyMesh& mesh,
+    scalar tol,
+    bool report
+)
+{
+    const labelList& pOwner = mesh.faceOwner();
+
+    // Find a suitable base point for each face, considering both
+    // cells for interface faces or those on coupled patches
+
+    labelList tetBasePtIs(mesh.nFaces(), -1);
+
+    label nInternalFaces = mesh.nInternalFaces();
+
+    for (label fI = 0; fI < nInternalFaces; fI++)
+    {
+        tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
+    }
+
+    pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
+
+    for(label faceI = nInternalFaces; faceI < mesh.nFaces(); faceI++)
+    {
+        neighbourCellCentres[faceI - nInternalFaces] =
+            mesh.cellCentres()[pOwner[faceI]];
+    }
+
+    syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    SubList<label> boundaryFaceTetBasePtIs
+    (
+        tetBasePtIs,
+        mesh.nFaces() - nInternalFaces,
+        nInternalFaces
+    );
+
+    for
+    (
+        label fI = nInternalFaces, bFI = 0;
+        fI < mesh.nFaces();
+        fI++, bFI++
+    )
+    {
+        label patchI =
+            mesh.boundaryMesh().patchID()[bFI];
+
+        if (patches[patchI].coupled())
+        {
+            const coupledPolyPatch& pp =
+                refCast<const coupledPolyPatch>(patches[patchI]);
+
+            if (pp.owner())
+            {
+                boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
+                (
+                    mesh,
+                    fI,
+                    neighbourCellCentres[bFI],
+                    tol,
+                    report
+                );
+            }
+            else
+            {
+                // Assign -2, to distinguish from a failed base point
+                // find, which returns -1.
+                boundaryFaceTetBasePtIs[bFI] = -2;
+            }
+        }
+        else
+        {
+            boundaryFaceTetBasePtIs[bFI] = findBasePoint
+            (
+                mesh,
+                fI,
+                tol,
+                report
+            );
+        }
+    }
+
+    // maxEqOp will replace the -2 values on the neighbour patches
+    // with the result from the owner base point find.
+
+    syncTools::syncBoundaryFaceList
+    (
+        mesh,
+        boundaryFaceTetBasePtIs,
+        maxEqOp<label>()
+    );
+
+    for
+    (
+        label fI = nInternalFaces, bFI = 0;
+        fI < mesh.nFaces();
+        fI++, bFI++
+    )
+    {
+        label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
+
+        if (bFTetBasePtI == -2)
+        {
+            FatalErrorIn
+            (
+                "Foam::labelList"
+                "Foam::polyMeshTetDecomposition::findFaceBasePts"
+                "("
+                    "const polyMesh& mesh, "
+                    "scalar tol, "
+                    "bool report"
+                ")"
+            )
+                << "Coupled face base point exchange failure for face "
+                << fI
+                << abort(FatalError);
+        }
+
+        if (bFTetBasePtI < 1)
+        {
+            // If the base point is -1, it should be left as such to
+            // indicate a problem, if it is 0, then no action is required.
+
+            continue;
+        }
+
+        label patchI = mesh.boundaryMesh().patchID()[bFI];
+
+        if (patches[patchI].coupled())
+        {
+            const coupledPolyPatch& pp =
+                refCast<const coupledPolyPatch>(patches[patchI]);
+
+            // Calculated base points on coupled faces are those of
+            // the owner patch face. They need to be reindexed to for
+            // the non-owner face, which has the opposite order.
+
+            // So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
+
+            // i.e.:
+
+            // owner coupledPolyPatch face
+            // face    (a b c d e f)
+            // fPtI     0 1 2 3 4 5
+            //            +
+            //              #
+
+            // neighbour coupledPolyPatch face
+            // face    (a f e d c b)
+            // fPtI     0 1 2 3 4 5
+            //                    +
+            //                  #
+            // +: 6 - 1 = 5
+            // #: 6 - 2 = 4
+
+            if (!pp.owner())
+            {
+                bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
+            }
+        }
+    }
+
+    return tetBasePtIs;
+}
+
+
+bool Foam::polyMeshTetDecomposition::checkFaceTets
+(
+    const polyMesh& mesh,
+    scalar tol,
+    const bool report,
+    labelHashSet* setPtr
+)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    const vectorField& cc = mesh.cellCentres();
+    const vectorField& fc = mesh.faceCentres();
+
+    // Calculate coupled cell centre
+    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+
+    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    {
+        neiCc[faceI - mesh.nInternalFaces()] = cc[own[faceI]];
+    }
+
+    syncTools::swapBoundaryFacePositions(mesh, neiCc);
+
+    const faceList& fcs = mesh.faces();
+
+    const pointField& p = mesh.points();
+
+    label nErrorTets = 0;
+
+    forAll(fcs, faceI)
+    {
+        const face& f = fcs[faceI];
+
+        forAll(f, fPtI)
+        {
+            scalar tetQual = tetPointRef
+            (
+                p[f[fPtI]],
+                p[f.nextLabel(fPtI)],
+                fc[faceI],
+                cc[own[faceI]]
+            ).quality();
+
+            if (tetQual > -tol)
+            {
+                if (setPtr)
+                {
+                    setPtr->insert(faceI);
+                }
+
+                nErrorTets++;
+                break;              // no need to check other tets
+            }
+        }
+
+        if (mesh.isInternalFace(faceI))
+        {
+            // Create the neighbour tet - it will have positive volume
+            const face& f = fcs[faceI];
+
+            forAll(f, fPtI)
+            {
+                scalar tetQual = tetPointRef
+                (
+                    p[f[fPtI]],
+                    p[f.nextLabel(fPtI)],
+                    fc[faceI],
+                    cc[nei[faceI]]
+                ).quality();
+
+                if (tetQual < tol)
+                {
+                    if (setPtr)
+                    {
+                        setPtr->insert(faceI);
+                    }
+
+                    nErrorTets++;
+                    break;
+                }
+            }
+
+            if (findSharedBasePoint(mesh, faceI, tol, report) == -1)
+            {
+                if (setPtr)
+                {
+                    setPtr->insert(faceI);
+                }
+
+                nErrorTets++;
+            }
+        }
+        else
+        {
+            label patchI = patches.patchID()[faceI - mesh.nInternalFaces()];
+
+            if (patches[patchI].coupled())
+            {
+                if
+                (
+                    findSharedBasePoint
+                    (
+                        mesh,
+                        faceI,
+                        neiCc[faceI - mesh.nInternalFaces()],
+                        tol,
+                        report
+                    ) == -1
+                )
+                {
+                    if (setPtr)
+                    {
+                        setPtr->insert(faceI);
+                    }
+
+                    nErrorTets++;
+                }
+            }
+            else
+            {
+                if (findBasePoint(mesh, faceI, tol, report) == -1)
+                {
+                    if (setPtr)
+                    {
+                        setPtr->insert(faceI);
+                    }
+
+                    nErrorTets++;
+                }
+            }
+        }
+    }
+
+    reduce(nErrorTets, sumOp<label>());
+
+    if (nErrorTets > 0)
+    {
+        if (report)
+        {
+            Info<< " ***Error in face tets: "
+                << nErrorTets << " faces with low quality or negative volume "
+                << "decomposition tets." << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (report)
+        {
+            Info<< "    Face tets OK." << endl;
+        }
+
+        return false;
+    }
+}
+
+
+Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
+(
+    const polyMesh& mesh,
+    label fI,
+    label cI
+)
+{
+    const faceList& pFaces = mesh.faces();
+    const labelList& pOwner = mesh.faceOwner();
+
+    const face& f = pFaces[fI];
+
+    label nTets = f.size() - 2;
+
+    List<tetIndices> faceTets(nTets);
+
+    bool own = (pOwner[fI] == cI);
+
+    label tetBasePtI = mesh.tetBasePtIs()[fI];
+
+    if (tetBasePtI == -1)
+    {
+        FatalErrorIn
+        (
+            "Foam::List<Foam::FixedList<Foam::label, 4> >"
+            "Foam::Cloud<ParticleType>::"
+            "faceTetIndices(label fI, label cI) const"
+        )
+        << "No base point for face " << fI << ", " << f
+            << ", produces a valid tet decomposition."
+            << abort(FatalError);
+    }
+
+    for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+    {
+        tetIndices& faceTetIs = faceTets[tetPtI - 1];
+
+        label facePtI = (tetPtI + tetBasePtI) % f.size();
+        label otherFacePtI = f.fcIndex(facePtI);
+
+        faceTetIs.cell() = cI;
+
+        faceTetIs.face() = fI;
+
+        faceTetIs.faceBasePt() = tetBasePtI;
+
+        if (own)
+        {
+            faceTetIs.facePtA() = facePtI;
+            faceTetIs.facePtB() = otherFacePtI;
+        }
+        else
+        {
+            faceTetIs.facePtA() = otherFacePtI;
+            faceTetIs.facePtB() = facePtI;
+        }
+
+        faceTetIs.tetPt() = tetPtI;
+    }
+
+    return faceTets;
+}
+
+
+Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
+(
+    const polyMesh& mesh,
+    label cI
+)
+{
+    const faceList& pFaces = mesh.faces();
+    const cellList& pCells = mesh.cells();
+
+    const cell& thisCell = pCells[cI];
+
+    label nTets = 0;
+
+    forAll(thisCell, cFI)
+    {
+        nTets += pFaces[thisCell[cFI]].size() - 2;
+    }
+
+    DynamicList<tetIndices> cellTets(nTets);
+
+    forAll(thisCell, cFI)
+    {
+        label fI = thisCell[cFI];
+
+        cellTets.append(faceTetIndices(mesh, fI, cI));
+    }
+
+    return cellTets;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyMeshTetDecomposition
+
+Description
+    Tools for performing the minimum decomposition of faces of the
+    mesh into triangles so that the cells may be tet decomposed.
+    Includes functions for finding variable face starting (base)
+    points on each face to avoid the decomposition of cells into tets
+    that have negative or zero volume.
+
+SourceFiles
+    polyMeshTetDecomposition.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyMeshTetDecomposition_H
+#define polyMeshTetDecomposition_H
+
+#include "polyMesh.H"
+#include "coupledPolyPatch.H"
+#include "syncTools.H"
+#include "tetPointRef.H"
+#include "tetIndices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class polyMeshTetDecomposition Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyMeshTetDecomposition
+{
+
+public:
+
+    // Static data members
+
+        //- Minimum tetrahedron quality
+        static const scalar minTetQuality;
+
+
+    // Member Functions
+
+        //- Find the first suitable base point to use for a minimum
+        //  triangle decomposition of the face, suiting owner and
+        //  neighbour cells.  Finds the first base point on the face
+        //  whose worst quality tet from either cell is better than
+        //  tolerance.  Neighbour cell centre supplied.  For coupled
+        //  patches.
+        static label findSharedBasePoint
+        (
+            const polyMesh& mesh,
+            label fI,
+            const point& nCc,
+            scalar tol,
+            bool report = false
+        );
+
+        //- As for findSharedBasePoint, but using neighbour cell
+        //  centre from the mesh.  For internal faces.
+        static label findSharedBasePoint
+        (
+            const polyMesh& mesh,
+            label fI,
+            scalar tol,
+            bool report = false
+        );
+
+        //- Find the base point to use for a minimum triangle
+        //  decomposition of the face, using only the owner
+        //  information.  For non-coupled boundary faces.
+        static label findBasePoint
+        (
+            const polyMesh& mesh,
+            label fI,
+            scalar tol,
+            bool report = false
+        );
+
+        //- Find a suitable base point for each face for decomposition
+        //  into tets
+        static labelList findFaceBasePts
+        (
+            const polyMesh& mesh,
+            scalar tol = minTetQuality,
+            bool report = false
+        );
+
+        //- Check face-decomposition tet volume
+        static bool checkFaceTets
+        (
+            const polyMesh& mesh,
+            scalar tol = minTetQuality,
+            const bool report = false,
+            labelHashSet* setPtr = NULL
+        );
+
+        //- Return the tet decomposition of the given face, with
+        //  respect to the given cell
+        static List<tetIndices> faceTetIndices
+        (
+            const polyMesh& mesh,
+            label fI,
+            label cI
+        );
+
+        //- Return the tet decomposition of the given cell, see
+        //  findFacePt for the meaning of the indices
+        static List<tetIndices> cellTetIndices
+        (
+            const polyMesh& mesh,
+            label cI
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.C

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetIndices.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::tetIndices::tetIndices()
+:
+    cellI_(-1),
+    faceI_(-1),
+    faceBasePtI_(-1),
+    facePtAI_(-1),
+    facePtBI_(-1),
+    tetPtI_(-1)
+{}
+
+
+Foam::tetIndices::tetIndices
+(
+    label cellI,
+    label faceI,
+    label faceBasePtI,
+    label facePtAI,
+    label facePtBI,
+    label tetPtI
+)
+:
+    cellI_(cellI),
+    faceI_(faceI),
+    faceBasePtI_(faceBasePtI),
+    facePtAI_(facePtAI),
+    facePtBI_(facePtBI),
+    tetPtI_(tetPtI)
+{}
+
+
+Foam::tetIndices::tetIndices
+(
+    label cellI,
+    label faceI,
+    label tetPtI,
+    const polyMesh& mesh
+)
+:
+    cellI_(cellI),
+    faceI_(faceI),
+    faceBasePtI_(-1),
+    facePtAI_(-1),
+    facePtBI_(-1),
+    tetPtI_(tetPtI)
+{
+    const faceList& pFaces = mesh.faces();
+    const labelList& pOwner = mesh.faceOwner();
+
+    const Foam::face& f = pFaces[faceI_];
+
+    bool own = (pOwner[faceI_] == cellI_);
+
+    faceBasePtI_ = mesh.tetBasePtIs()[faceI_];
+
+    label facePtI = (tetPtI_ + faceBasePtI_) % f.size();
+    label otherFacePtI = f.fcIndex(facePtI);
+
+    if (own)
+    {
+        facePtAI_ = facePtI;
+        facePtBI_ = otherFacePtI;
+    }
+    else
+    {
+        facePtAI_ = otherFacePtI;
+        facePtBI_ = facePtI;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::tetIndices::~tetIndices()
+{}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, tetIndices& tI)
+{
+    is  >> tI.cell()
+        >> tI.face()
+        >> tI.faceBasePt()
+        >> tI.facePtA()
+        >> tI.facePtB()
+        >> tI.tetPt();
+
+    // Check state of Istream
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::tetIndices&)"
+    );
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const tetIndices& tI)
+{
+    os  << tI.cell() << token::SPACE
+        << tI.face() << token::SPACE
+        << tI.faceBasePt() << token::SPACE
+        << tI.facePtA() << token::SPACE
+        << tI.facePtB() << token::SPACE
+        << tI.tetPt() << token::SPACE
+        << endl;
+
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
+        "const Foam::tetIndices&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H

@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::tetIndices
+
+Description
+    Storage and named access for the indices of a tet which is part of
+    the decomposition of a cell.
+
+SourceFiles
+    tetIndicesI.H
+    tetIndices.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetIndices_H
+#define tetIndices_H
+
+#include "label.H"
+#include "tetPointRef.H"
+#include "triPointRef.H"
+#include "polyMesh.H"
+#include "face.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class tetIndices Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetIndices
+{
+    // Private data
+
+    //- cell that this is a decomposed tet of
+    label cellI_;
+
+    //- face that holds this decomposed tet
+    label faceI_;
+
+    //- base point on the face
+    label faceBasePtI_;
+
+    //- point on the face such that the right-hand circulation
+    //      {faceBasePtI_, facePtIA_, facePtBI_}
+    //  forms a triangle that points out of the tet
+    label facePtAI_;
+
+    //- point on the face such that the right-hand circulation
+    //      {faceBasePtI_, facePtIA_, facePtBI_}
+    //  forms a triangle that points out of the tet
+    label facePtBI_;
+
+    //- point on the face, *relative to the base point*, which
+    //  characterises this tet on the face.
+    label tetPtI_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        tetIndices();
+
+        //- Construct from components
+        tetIndices
+        (
+            label cellI,
+            label faceI,
+            label faceBasePtI,
+            label facePtAI,
+            label facePtBI,
+            label tetPtI
+        );
+
+        //- Construct from cell, face, pt description of tet
+        tetIndices
+        (
+            label cellI,
+            label faceI,
+            label tetPtI,
+            const polyMesh& mesh
+        );
+
+
+    //- Destructor
+    ~tetIndices();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the cell
+            inline label cell() const;
+
+            //- Return the face
+            inline label face() const;
+
+            //- Return the face base point
+            inline label faceBasePt() const;
+
+            //- Return face point A
+            inline label facePtA() const;
+
+            //- Return face point B
+            inline label facePtB() const;
+
+            //- Return the characterising tetPtI
+            inline label tetPt() const;
+
+            //- Return the geometry corresponding to this tet from the
+            //  supplied mesh
+            inline tetPointRef tet(const polyMesh& mesh) const;
+
+            //- Return the geometry corresponding to this tet from the
+            //  supplied mesh using the old positions
+            inline tetPointRef oldTet(const polyMesh& mesh) const;
+
+            //- Return the geometry corresponding to the tri on the
+            //  mesh face for this tet from the supplied mesh
+            inline triPointRef faceTri(const polyMesh& mesh) const;
+
+            //- Return the geometry corresponding to the tri on the
+            //  mesh face for this tet from the supplied mesh using
+            //  the old position
+            inline triPointRef oldFaceTri(const polyMesh& mesh) const;
+
+
+        // Edit
+
+            //- Return non-const access to the cell
+            inline label& cell();
+
+            //- Return non-const access to the face
+            inline label& face();
+
+            //- Return non-const access to the face base point
+            inline label& faceBasePt();
+
+            //- Return non-const access to face point A
+            inline label& facePtA();
+
+            //- Return non-const access to face point B
+            inline label& facePtB();
+
+            //- Return non-const access to the characterising tetPtI
+            inline label& tetPt();
+
+
+    // Member Operators
+
+        inline bool operator==(const tetIndices&) const;
+        inline bool operator!=(const tetIndices&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, tetIndices&);
+        friend Ostream& operator<<(Ostream&, const tetIndices&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "tetIndicesI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H

@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::tetIndices::cell() const
+{
+    return cellI_;
+}
+
+
+Foam::label Foam::tetIndices::face() const
+{
+    return faceI_;
+}
+
+
+Foam::label Foam::tetIndices::faceBasePt() const
+{
+    return faceBasePtI_;
+}
+
+
+Foam::label Foam::tetIndices::facePtA() const
+{
+    return facePtAI_;
+}
+
+
+Foam::label Foam::tetIndices::facePtB() const
+{
+    return facePtBI_;
+}
+
+
+Foam::label Foam::tetIndices::tetPt() const
+{
+    return tetPtI_;
+}
+
+
+Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const
+{
+    const pointField& pPts = mesh.points();
+    const faceList& pFaces = mesh.faces();
+    const vectorField& pC = mesh.cellCentres();
+
+    const Foam::face& f = pFaces[faceI_];
+
+    return tetPointRef
+    (
+        pC[cellI_],
+        pPts[f[faceBasePtI_]],
+        pPts[f[facePtAI_]],
+        pPts[f[facePtBI_]]
+    );
+}
+
+
+Foam::tetPointRef Foam::tetIndices::oldTet(const polyMesh& mesh) const
+{
+    const pointField& oldPPts = mesh.oldPoints();
+    const faceList& pFaces = mesh.faces();
+
+    // We need to reconstruct the old Cc from oldPoints (it isn't
+    // stored)
+    point oldC = mesh.cells()[cellI_].centre
+    (
+        oldPPts,
+        pFaces
+    );
+
+    const Foam::face& f = pFaces[faceI_];
+
+    return tetPointRef
+    (
+        oldC,
+        oldPPts[f[faceBasePtI_]],
+        oldPPts[f[facePtAI_]],
+        oldPPts[f[facePtBI_]]
+    );
+}
+
+
+Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const
+{
+    const pointField& pPts = mesh.points();
+    const faceList& pFaces = mesh.faces();
+
+    const Foam::face& f = pFaces[faceI_];
+
+    return triPointRef
+    (
+        pPts[f[faceBasePtI_]],
+        pPts[f[facePtAI_]],
+        pPts[f[facePtBI_]]
+    );
+}
+
+
+Foam::triPointRef Foam::tetIndices::oldFaceTri(const polyMesh& mesh) const
+{
+    const pointField& oldPPts = mesh.oldPoints();
+    const faceList& pFaces = mesh.faces();
+
+    const Foam::face& f = pFaces[faceI_];
+
+    return triPointRef
+    (
+        oldPPts[f[faceBasePtI_]],
+        oldPPts[f[facePtAI_]],
+        oldPPts[f[facePtBI_]]
+    );
+}
+
+
+Foam::label& Foam::tetIndices::cell()
+{
+    return cellI_;
+}
+
+
+Foam::label& Foam::tetIndices::face()
+{
+    return faceI_;
+}
+
+
+Foam::label& Foam::tetIndices::faceBasePt()
+{
+    return faceBasePtI_;
+}
+
+
+Foam::label& Foam::tetIndices::facePtA()
+{
+    return facePtAI_;
+}
+
+
+Foam::label& Foam::tetIndices::facePtB()
+{
+    return facePtBI_;
+}
+
+
+Foam::label& Foam::tetIndices::tetPt()
+{
+    return tetPtI_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::tetIndices::operator==(const Foam::tetIndices& rhs) const
+{
+    return
+    (
+        cell() == rhs.cell()
+     && face() == rhs.face()
+     && faceBasePt() == rhs.faceBasePt()
+     && facePtA() == rhs.facePtA()
+     && facePtB() == rhs.facePtB()
+     && tetPt() == rhs.tetPt()
+    );
+}
+
+
+inline bool Foam::tetIndices::operator!=(const Foam::tetIndices& rhs) const
+{
+    return !(*this == rhs);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H

@@ -229,6 +229,19 @@ public:
 
     // Member functions
 
+        //- Return true only if this is a parallel run
+        virtual bool coupled() const
+        {
+            if (Pstream::parRun())
+            {
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+
         //- Return processor number
         int myProcNo() const
         {

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -66,6 +66,7 @@ SourceFiles
 #include "HashSet.H"
 #include "Map.H"
 #include "EdgeMap.H"
+#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -577,14 +578,6 @@ public:
                     labelHashSet* setPtr = NULL
                 ) const;
 
-                //- Check face-decomposition tet volume
-                bool checkFaceTets
-                (
-                    const bool report = false,
-                    const scalar minTetVol = 0,
-                    labelHashSet* setPtr = NULL
-                ) const;
-
                 //- Check face skewness
                 bool checkFaceSkewness
                 (

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -618,114 +618,6 @@ bool Foam::primitiveMesh::checkFacePyramids
 }
 
 
-bool Foam::primitiveMesh::checkFaceTets
-(
-    const bool report,
-    const scalar minTetVol,
-    labelHashSet* setPtr
-) const
-{
-    if (debug)
-    {
-        Info<< "bool primitiveMesh::checkFaceTets("
-            << "const bool, const scalar, labelHashSet*) const: "
-            << "checking face orientation" << endl;
-    }
-
-    // check whether face area vector points to the cell with higher label
-    const vectorField& cc = cellCentres();
-    const vectorField& fc = faceCentres();
-
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
-
-    const faceList& fcs = faces();
-
-    const pointField& p = points();
-
-    label nErrorPyrs = 0;
-
-    forAll(fcs, faceI)
-    {
-        // Create the owner tets - they will have negative volume
-        const face& f = fcs[faceI];
-
-        forAll(f, fp)
-        {
-            scalar tetVol = tetPointRef
-            (
-                p[f[fp]],
-                p[f.nextLabel(fp)],
-                fc[faceI],
-                cc[own[faceI]]
-            ).mag();
-
-            if (tetVol > -minTetVol)
-            {
-                if (setPtr)
-                {
-                    setPtr->insert(faceI);
-                }
-
-                nErrorPyrs++;
-                break;              // no need to check other tets
-            }
-        }
-
-        if (isInternalFace(faceI))
-        {
-            // Create the neighbour tet - it will have positive volume
-            const face& f = fcs[faceI];
-
-            forAll(f, fp)
-            {
-                scalar tetVol = tetPointRef
-                (
-                    p[f[fp]],
-                    p[f.nextLabel(fp)],
-                    fc[faceI],
-                    cc[nei[faceI]]
-                ).mag();
-
-                if (tetVol < minTetVol)
-                {
-                    if (setPtr)
-                    {
-                        setPtr->insert(faceI);
-                    }
-
-                    nErrorPyrs++;
-                    break;
-                }
-            }
-        }
-    }
-
-    reduce(nErrorPyrs, sumOp<label>());
-
-    if (nErrorPyrs > 0)
-    {
-        if (debug || report)
-        {
-            Info<< " ***Error in face tets: "
-                << nErrorPyrs << " faces have incorrectly oriented face"
-                << " decomposition triangles." << endl;
-        }
-
-        return true;
-    }
-    else
-    {
-        if (debug || report)
-        {
-            Info<< "    Face tets OK." << endl;
-        }
-
-        return false;
-    }
-}
-
-
 bool Foam::primitiveMesh::checkFaceSkewness
 (
     const bool report,

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C

@@ -32,32 +32,14 @@ License
 // Is the point in the cell bounding box
 bool Foam::primitiveMesh::pointInCellBB(const point& p, label celli) const
 {
-    const pointField& points = this->points();
-    const faceList& f = faces();
-    const vectorField& centres = cellCentres();
-    const cellList& cf = cells();
-
-    labelList cellVertices = cf[celli].labels(f);
-
-    vector bbmax = -GREAT*vector::one;
-    vector bbmin = GREAT*vector::one;
-
-    forAll(cellVertices, vertexI)
-    {
-        bbmax = max(bbmax, points[cellVertices[vertexI]]);
-        bbmin = min(bbmin, points[cellVertices[vertexI]]);
-    }
-
-    scalar distance = mag(centres[celli] - p);
-
-    if ((distance - mag(bbmax - bbmin)) < SMALL)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
+    return boundBox
+    (
+        cells()[celli].points
+        (
+            faces(),
+            points()
+        )
+    ).contains(p);
 }
 
 

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C

@@ -30,7 +30,6 @@ Description
 #include "triPointRef.H"
 #include "scalarField.H"
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // (Probably very inefficient) minimum containment sphere calculation.

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H

@@ -42,6 +42,7 @@ SourceFiles
 #include "point.H"
 #include "primitiveFieldsFwd.H"
 #include "pointHit.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -133,16 +134,42 @@ public:
 
             inline vector Sd() const;
 
+            //- Return centre (centroid)
+            inline Point centre() const;
 
             //- Return volume
             inline scalar mag() const;
 
             //- Return circum-centre
-            inline vector circumCentre() const;
+            inline Point circumCentre() const;
 
             //- Return circum-radius
             inline scalar circumRadius() const;
 
+            //- Return quality: Ratio of tetrahedron and circum-sphere
+            //  volume, scaled so that a regular tetrahedron has a
+            //  quality of 1
+            inline scalar quality() const;
+
+            //- Return a random point in the tetrahedron from a
+            //  uniform distribution
+            inline Point randomPoint(Random& rndGen) const;
+
+            //- Calculate the barycentric coordinates of the given
+            //  point, in the same order as a, b, c, d.  Returns the
+            //  determinant of the solution.
+            inline scalar barycentric
+            (
+                const point& pt,
+                List<scalar>& bary
+            ) const;
+
+            //- Return nearest point to p on tetrahedron
+            inline pointHit nearestPoint
+            (
+                const point& p
+            ) const;
+
             //- Return (min)containment sphere, i.e. the smallest sphere with
             //  all points inside. Returns pointHit with:
             //  - hit         : if sphere is equal to circumsphere

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -125,6 +125,13 @@ inline vector tetrahedron<Point, PointRef>::Sd() const
 }
 
 
+template<class Point, class PointRef>
+inline Point tetrahedron<Point, PointRef>::centre() const
+{
+    return 0.25*(a_ + b_ + c_ + d_);
+}
+
+
 template<class Point, class PointRef>
 inline scalar tetrahedron<Point, PointRef>::mag() const
 {
@@ -133,19 +140,19 @@ inline scalar tetrahedron<Point, PointRef>::mag() const
 
 
 template<class Point, class PointRef>
-inline vector tetrahedron<Point, PointRef>::circumCentre() const
+inline Point tetrahedron<Point, PointRef>::circumCentre() const
 {
     vector a = b_ - a_;
     vector b = c_ - a_;
     vector c = d_ - a_;
 
-    scalar lamda = magSqr(c) - (a&c);
+    scalar lambda = magSqr(c) - (a & c);
     scalar mu = magSqr(b) - (a & b);
 
     vector ba = b ^ a;
     vector ca = c ^ a;
 
-    return a_ + 0.5*(a + (lamda*ba - mu*ca)/(c&ba));
+    return a_ + 0.5*(a + (lambda*ba - mu*ca)/((c & ba) + ROOTVSMALL));
 }
 
 
@@ -156,6 +163,196 @@ inline scalar tetrahedron<Point, PointRef>::circumRadius() const
 }
 
 
+template<class Point, class PointRef>
+inline scalar tetrahedron<Point, PointRef>::quality() const
+{
+    // Note: 8/(9*sqrt(3)) = 0.5132002393
+
+    return mag()/(0.5132002393*pow3(min(circumRadius(), GREAT)) + ROOTVSMALL);
+}
+
+
+template<class Point, class PointRef>
+inline Point tetrahedron<Point, PointRef>::randomPoint(Random& rndGen) const
+{
+    // Adapted from
+    // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
+
+    scalar s = rndGen.scalar01();
+    scalar t = rndGen.scalar01();
+    scalar u = rndGen.scalar01();
+
+    if (s + t > 1.0)
+    {
+        s = 1.0 - s;
+        t = 1.0 - t;
+    }
+
+    if (t + u > 1.0)
+    {
+        scalar tmp = u;
+        u = 1.0 - s - t;
+        t = 1.0 - tmp;
+    }
+    else if (s + t + u > 1.0)
+    {
+        scalar tmp = u;
+        u = s + t + u - 1.0;
+        s = 1.0 - t - tmp;
+    }
+
+    return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
+}
+
+
+template<class Point, class PointRef>
+scalar tetrahedron<Point, PointRef>::barycentric
+(
+    const point& pt,
+    List<scalar>& bary
+) const
+{
+    // From:
+    // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics)
+
+    vector e0(a_ - d_);
+    vector e1(b_ - d_);
+    vector e2(c_ - d_);
+
+    tensor t
+    (
+        e0.x(), e1.x(), e2.x(),
+        e0.y(), e1.y(), e2.y(),
+        e0.z(), e1.z(), e2.z()
+    );
+
+    scalar detT = det(t);
+
+    if (Foam::mag(detT) < SMALL)
+    {
+        WarningIn
+        (
+            "List<scalar> tetrahedron<Point, PointRef>::barycentric"
+            "("
+                "const point& pt"
+            ") const"
+        )
+            << "Degenerate triangle - returning 1/4 barycentric coordinates."
+            << endl;
+
+        bary = List<scalar>(4, 0.25);
+
+        return detT;
+    }
+
+    vector res = inv(t, detT) & (pt - d_);
+
+    bary.setSize(4);
+
+    bary[0] = res.x();
+    bary[1] = res.y();
+    bary[2] = res.z();
+    bary[3] = (1.0 - res.x() - res.y() - res.z());
+
+    return detT;
+}
+
+
+template<class Point, class PointRef>
+inline pointHit tetrahedron<Point, PointRef>::nearestPoint
+(
+    const point& p
+) const
+{
+    // Adapted from:
+    // Real-time collision detection, Christer Ericson, 2005, p142-144
+
+    // Assuming initially that the point is inside all of the
+    // halfspaces, so closest to itself.
+
+    point closestPt = p;
+
+    scalar minOutsideDistance = VGREAT;
+
+    bool inside = true;
+
+    if (((p - b_) & Sa()) >= 0)
+    {
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
+
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
+        {
+            closestPt = info.rawPoint();
+
+            minOutsideDistance = info.distance();
+        }
+    }
+
+    if (((p - a_) & Sb()) >= 0)
+    {
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
+
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
+        {
+            closestPt = info.rawPoint();
+
+            minOutsideDistance = info.distance();
+        }
+    }
+
+    if (((p - a_) & Sc()) >= 0)
+    {
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
+
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
+        {
+            closestPt = info.rawPoint();
+
+            minOutsideDistance = info.distance();
+        }
+    }
+
+    if (((p - a_) & Sd()) >= 0)
+    {
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
+
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
+        {
+            closestPt = info.rawPoint();
+
+            minOutsideDistance = info.distance();
+        }
+    }
+
+    // If the point is inside, then the distance to the closest point
+    // is zero
+    if (inside)
+    {
+        minOutsideDistance = 0;
+    }
+
+    return pointHit
+    (
+        inside,
+        closestPt,
+        minOutsideDistance,
+        !inside
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class point, class pointRef>

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

@@ -39,7 +39,7 @@ SourceFiles
 #include "vector.H"
 #include "tensor.H"
 #include "pointHit.H"
-
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -141,12 +141,14 @@ public:
             inline vector normal() const;
 
             //- Return circum-centre
-            inline vector circumCentre() const;
+            inline Point circumCentre() const;
 
             //- Return circum-radius
             inline scalar circumRadius() const;
 
-            //- Return quality: Ratio triangle and circum-circle area
+            //- Return quality: Ratio of triangle and circum-circle
+            //  area, scaled so that an equilateral triangle has a
+            //  quality of 1
             inline scalar quality() const;
 
             //- Return swept-volume
@@ -160,6 +162,19 @@ public:
                 scalar density = 1.0
             ) const;
 
+            //- Return a random point on the triangle from a uniform
+            //  distribution
+            inline Point randomPoint(Random& rndGen) const;
+
+            //- Calculate the barycentric coordinates of the given
+            //  point, in the same order as a, b, c.  Returns the
+            //  determinant of the solution.
+            inline scalar barycentric
+            (
+                const point& pt,
+                List<scalar>& bary
+            ) const;
+
             //- Return point intersection with a ray.
             //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.

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

@@ -259,7 +259,7 @@ inline vector triangle<Point, PointRef>::normal() const
 
 
 template<class Point, class PointRef>
-inline vector triangle<Point, PointRef>::circumCentre() const
+inline Point triangle<Point, PointRef>::circumCentre() const
 {
     scalar d1 = (c_ - a_)&(b_ - a_);
     scalar d2 = -(c_ - b_)&(b_ - a_);
@@ -303,12 +303,14 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
 template<class Point, class PointRef>
 inline scalar triangle<Point, PointRef>::quality() const
 {
+    // Note: 3*sqr(3)/(4*pi) = 0.4134966716
+
     return
         mag()
       / (
             constant::mathematical::pi
            *Foam::sqr(circumRadius())
-           *0.413497
+           *0.4134966716
           + VSMALL
         );
 }
@@ -366,6 +368,72 @@ inline tensor triangle<Point, PointRef>::inertia
    *density;
 }
 
+
+template<class Point, class PointRef>
+inline Point triangle<Point, PointRef>::randomPoint(Random& rndGen) const
+{
+    // Generating Random Points in Triangles
+    // by Greg Turk
+    // from "Graphics Gems", Academic Press, 1990
+    // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
+
+    scalar s = rndGen.scalar01();
+
+    scalar t = sqrt(rndGen.scalar01());
+
+    return (1 - t)*a_ + (1 - s)*t*b_ + s*t*c_;
+}
+
+
+template<class Point, class PointRef>
+scalar triangle<Point, PointRef>::barycentric
+(
+    const point& pt,
+    List<scalar>& bary
+) const
+{
+    // From:
+    // Real-time collision detection, Christer Ericson, 2005, p47-48
+
+    vector v0 = b_ - a_;
+    vector v1 = c_ - a_;
+    vector v2 = pt - a_;
+
+    scalar d00 = v0 & v0;
+    scalar d01 = v0 & v1;
+    scalar d11 = v1 & v1;
+    scalar d20 = v2 & v0;
+    scalar d21 = v2 & v1;
+
+    scalar denom = d00*d11 - d01*d01;
+
+    if (Foam::mag(denom) < SMALL)
+    {
+        WarningIn
+        (
+            "List<scalar> triangle<Point, PointRef>::barycentric"
+            "("
+                "const point& pt"
+            ") const"
+        )
+            << "Degenerate triangle - returning 1/3 barycentric coordinates."
+            << endl;
+
+        bary = List<scalar>(3, 1.0/3.0);
+
+        return denom;
+    }
+
+    bary.setSize(3);
+
+    bary[0] = (d11*d20 - d01*d21)/denom;
+    bary[1] = (d00*d21 - d01*d20)/denom;
+    bary[2] = 1.0 - bary[0] - bary[1];
+
+    return denom;
+}
+
+
 template<class Point, class PointRef>
 inline pointHit triangle<Point, PointRef>::ray
 (
@@ -775,4 +843,3 @@ inline Ostream& operator<<(Ostream& os, const triangle<Point, PointRef>& t)
 } // End namespace Foam
 
 // ************************************************************************* //
-

src/OpenFOAM/primitives/Tensor/TensorI.H

@@ -481,7 +481,7 @@ inline Tensor<Cmpt> cof(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the inverse of a tensor give the determinant
+//- Return the inverse of a tensor given the determinant
 template <class Cmpt>
 inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt dett)
 {

src/dynamicMesh/motionSmoother/motionSmootherCheck.C

@@ -68,9 +68,9 @@ bool Foam::motionSmoother::checkMesh
     (
         readScalar(dict.lookup("minVol", true))
     );
-    const scalar minTetVol
+    const scalar minTetQuality
     (
-        readScalar(dict.lookup("minTetVol", true))
+        readScalar(dict.lookup("minTetQuality", true))
     );
     const scalar maxConcave
     (
@@ -160,12 +160,12 @@ bool Foam::motionSmoother::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
-    if (minTetVol > -GREAT)
+    if (minTetQuality > -GREAT)
     {
         polyMeshGeometry::checkFaceTets
         (
             report,
-            minTetVol,
+            minTetQuality,
             mesh,
             mesh.cellCentres(),
             mesh.faceCentres(),
@@ -177,8 +177,8 @@ bool Foam::motionSmoother::checkMesh
 
         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 
-        Info<< "    faces with face-decomposition tet volume < "
-            << setw(5) << minTetVol << "        : "
+        Info<< "    faces with face-decomposition tet quality < "
+            << setw(5) << minTetQuality << "      : "
             << nNewWrongFaces-nWrongFaces << endl;
 
         nWrongFaces = nNewWrongFaces;
@@ -443,9 +443,9 @@ bool Foam::motionSmoother::checkMesh
     (
         readScalar(dict.lookup("minVol", true))
     );
-    const scalar minTetVol
+    const scalar minTetQuality
     (
-        readScalar(dict.lookup("minTetVol", true))
+        readScalar(dict.lookup("minTetQuality", true))
     );
     const scalar maxConcave
     (
@@ -531,12 +531,12 @@ bool Foam::motionSmoother::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
-    if (minTetVol > -GREAT)
+    if (minTetQuality > -GREAT)
     {
         meshGeom.checkFaceTets
         (
             report,
-            minTetVol,
+            minTetQuality,
             meshGeom.mesh().points(),
             checkFaces,
             baffles,
@@ -545,8 +545,8 @@ bool Foam::motionSmoother::checkMesh
 
         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 
-        Info<< "    faces with face-decomposition tet volume < "
-            << setw(5) << minTetVol << "                 : "
+        Info<< "    faces with face-decomposition tet quality < "
+            << setw(5) << minTetQuality << "                : "
             << nNewWrongFaces-nWrongFaces << endl;
 
         nWrongFaces = nNewWrongFaces;

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "polyMeshGeometry.H"
+#include "polyMeshTetDecomposition.H"
 #include "pyramidPointFaceRef.H"
 #include "tetPointRef.H"
 #include "syncTools.H"
@@ -313,7 +314,7 @@ bool Foam::polyMeshGeometry::checkFaceTet
 (
     const polyMesh& mesh,
     const bool report,
-    const scalar minTetVol,
+    const scalar minTetQuality,
     const pointField& p,
     const label faceI,
     const point& fc,    // face centre
@@ -326,15 +327,15 @@ bool Foam::polyMeshGeometry::checkFaceTet
 
     forAll(f, fp)
     {
-        scalar tetVol = tetPointRef
+        scalar tetQual = tetPointRef
         (
             p[f[fp]],
             p[f.nextLabel(fp)],
             fc,
             cc
-        ).mag();
+        ).quality();
 
-        if (tetVol < minTetVol)
+        if (tetQual < minTetQuality)
         {
             if (report)
             {
@@ -345,7 +346,7 @@ bool Foam::polyMeshGeometry::checkFaceTet
                     << "face " << faceI
                     << " has a triangle that points the wrong way."
                      << endl
-                    << "Tet volume: " << tetVol
+                    << "Tet quality: " << tetQual
                     << " Face " << faceI
                     << endl;
             }
@@ -423,7 +424,6 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
     // Severe nonorthogonality threshold
     const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
 
-
     // Calculate coupled cell centre
     pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
 
@@ -431,8 +431,8 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
     {
         neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
     }
-    syncTools::swapBoundaryFacePositions(mesh, neiCc);
 
+    syncTools::swapBoundaryFacePositions(mesh, neiCc);
 
     scalar minDDotS = GREAT;
 
@@ -754,7 +754,7 @@ bool Foam::polyMeshGeometry::checkFacePyramids
                 "polyMeshGeometry::checkFacePyramids("
                 "const bool, const scalar, const pointField&"
                 ", const labelList&, labelHashSet*)"
-            )   << "Error in face pyramids: faces pointing the wrong way!"
+            )   << "Error in face pyramids: faces pointing the wrong way."
                 << endl;
         }
 
@@ -775,7 +775,7 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 bool Foam::polyMeshGeometry::checkFaceTets
 (
     const bool report,
-    const scalar minTetVol,
+    const scalar minTetQuality,
     const polyMesh& mesh,
     const vectorField& cellCentres,
     const vectorField& faceCentres,
@@ -785,11 +785,23 @@ bool Foam::polyMeshGeometry::checkFaceTets
     labelHashSet* setPtr
 )
 {
-    // check whether face area vector points to the cell with higher label
+    // check whether decomposing each cell into tets results in
+    // positive volume, non-flat tets
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    label nErrorPyrs = 0;
+    // Calculate coupled cell centre
+    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+
+    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    {
+        neiCc[faceI - mesh.nInternalFaces()] = cellCentres[own[faceI]];
+    }
+
+    syncTools::swapBoundaryFacePositions(mesh, neiCc);
+
+    label nErrorTets = 0;
 
     forAll(checkFaces, i)
     {
@@ -801,7 +813,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
         (
             mesh,
             report,
-            minTetVol,
+            minTetQuality,
             p,
             faceI,
             cellCentres[own[faceI]],    // face centre
@@ -811,26 +823,94 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
         if (tetError)
         {
-            nErrorPyrs++;
+            nErrorTets++;
         }
 
         if (mesh.isInternalFace(faceI))
         {
-            // Create the neighbour pyramid - it will have positive volume
+            // Create the neighbour tets - they will have positive volume
             bool tetError = checkFaceTet
             (
                 mesh,
                 report,
-                minTetVol,
+                minTetQuality,
                 p,
                 faceI,
                 faceCentres[faceI],         // face centre
                 cellCentres[nei[faceI]],    // cell centre
                 setPtr
             );
+
             if (tetError)
             {
-                nErrorPyrs++;
+                nErrorTets++;
+            }
+
+            if
+            (
+                polyMeshTetDecomposition::findSharedBasePoint
+                (
+                    mesh,
+                    faceI,
+                    minTetQuality,
+                    report
+                ) == -1
+            )
+            {
+                if (setPtr)
+                {
+                    setPtr->insert(faceI);
+                }
+
+                nErrorTets++;
+            }
+        }
+        else
+        {
+            label patchI = patches.whichPatch(faceI);
+
+            if (patches[patchI].coupled())
+            {
+                if
+                (
+                    polyMeshTetDecomposition::findSharedBasePoint
+                    (
+                        mesh,
+                        faceI,
+                        neiCc[faceI - mesh.nInternalFaces()],
+                        minTetQuality,
+                        report
+                    ) == -1
+                )
+                {
+                    if (setPtr)
+                    {
+                        setPtr->insert(faceI);
+                    }
+
+                    nErrorTets++;
+                }
+            }
+            else
+            {
+                if
+                (
+                    polyMeshTetDecomposition::findBasePoint
+                    (
+                        mesh,
+                        faceI,
+                        minTetQuality,
+                        report
+                    ) == -1
+                )
+                {
+                    if (setPtr)
+                    {
+                        setPtr->insert(faceI);
+                    }
+
+                    nErrorTets++;
+                }
             }
         }
     }
@@ -844,7 +924,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
         (
             mesh,
             report,
-            minTetVol,
+            minTetQuality,
             p,
             face0,
             cellCentres[own[face0]],    // face centre
@@ -854,15 +934,15 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
         if (tetError)
         {
-            nErrorPyrs++;
+            nErrorTets++;
         }
 
-        // Create the neighbour pyramid - it will have positive volume
+        // Create the neighbour tets - they will have positive volume
         tetError = checkFaceTet
         (
             mesh,
             report,
-            minTetVol,
+            minTetQuality,
             p,
             face0,
             faceCentres[face0],         // face centre
@@ -872,13 +952,33 @@ bool Foam::polyMeshGeometry::checkFaceTets
 
         if (tetError)
         {
-            nErrorPyrs++;
+            nErrorTets++;
+        }
+
+        if
+        (
+            polyMeshTetDecomposition::findSharedBasePoint
+            (
+                mesh,
+                face0,
+                cellCentres[own[face1]],
+                minTetQuality,
+                report
+            ) == -1
+        )
+        {
+            if (setPtr)
+            {
+                setPtr->insert(face0);
+            }
+
+            nErrorTets++;
         }
     }
 
-    reduce(nErrorPyrs, sumOp<label>());
+    reduce(nErrorTets, sumOp<label>());
 
-    if (nErrorPyrs > 0)
+    if (nErrorTets > 0)
     {
         if (report)
         {
@@ -887,7 +987,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
                 "polyMeshGeometry::checkFaceTets("
                 "const bool, const scalar, const pointField&, const pointField&"
                 ", const labelList&, labelHashSet*)"
-            )   << "Error in face pyramids: faces pointing the wrong way!"
+            )   << "Error in face decomposition: negative tets."
                 << endl;
         }
 
@@ -2135,7 +2235,7 @@ bool Foam::polyMeshGeometry::checkFacePyramids
 bool Foam::polyMeshGeometry::checkFaceTets
 (
     const bool report,
-    const scalar minTetVol,
+    const scalar minTetQuality,
     const pointField& p,
     const labelList& checkFaces,
     const List<labelPair>& baffles,
@@ -2145,7 +2245,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
     return checkFaceTets
     (
         report,
-        minTetVol,
+        minTetQuality,
         mesh_,
         cellCentres_,
         faceCentres_,

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H

@@ -113,7 +113,7 @@ class polyMeshGeometry
         (
             const polyMesh&,
             const bool report,
-            const scalar minTetVol,
+            const scalar minTetQuality,
             const pointField& p,
             const label faceI,
             const point& fc,    // face centre
@@ -121,6 +121,7 @@ class polyMeshGeometry
             labelHashSet* setPtr
         );
 
+
 public:
 
     ClassName("polyMeshGeometry");
@@ -350,7 +351,7 @@ public:
             bool checkFaceTets
             (
                 const bool report,
-                const scalar minTetVol,
+                const scalar minTetQuality,
                 const pointField& p,
                 const labelList& checkFaces,
                 const List<labelPair>& baffles,

src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H

@@ -38,6 +38,7 @@ Description
 #include "typeInfo.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
+#include "tetIndices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -130,9 +131,23 @@ public:
         virtual Type interpolate
         (
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         ) const = 0;
+
+        //- Interpolate field to the given point in the tetrahedron
+        //  defined by the given indices.  Calls interpolate function
+        //  above here execpt where overridden by derived
+        //  interpolation types.
+        virtual Type interpolate
+        (
+            const vector& position,
+            const tetIndices& tetIs,
+            const label faceI = -1
+        ) const
+        {
+            return interpolate(position, tetIs.cell(), faceI);
+        }
 };
 
 

src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.C

@@ -49,11 +49,11 @@ template<class Type>
 Type interpolationCell<Type>::interpolate
 (
     const vector&,
-    const label celli,
+    const label cellI,
     const label
 ) const
 {
-    return this->psi_[celli];
+    return this->psi_[cellI];
 }
 
 

src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H

@@ -72,8 +72,8 @@ public:
         Type interpolate
         (
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         ) const;
 };
 

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C

@@ -25,10 +25,13 @@ License
 
 #include "cellPointWeight.H"
 #include "polyMesh.H"
+#include "tetPointRef.H"
+#include "polyMeshTetDecomposition.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 int Foam::cellPointWeight::debug(debug::debugSwitch("cellPointWeight", 0));
+
 Foam::scalar Foam::cellPointWeight::tol(SMALL);
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
@@ -37,148 +40,101 @@ void Foam::cellPointWeight::findTetrahedron
 (
     const polyMesh& mesh,
     const vector& position,
-    const label cellIndex
+    const label cellI
 )
 {
     if (debug)
     {
-        Pout<< "\nFoam::cellPointWeight::findTetrahedron" << nl
+        Pout<< nl << "Foam::cellPointWeight::findTetrahedron" << nl
             << "position = " << position << nl
-            << "cellIndex = " << cellIndex << endl;
+            << "cellI = " << cellI << endl;
     }
 
-    // Initialise closest triangle variables
-    scalar minUVWClose = VGREAT;
-    label pointIClose = 0;
-    label faceClose = 0;
-
-    const vector& P0 = mesh.cellCentres()[cellIndex];
-    const labelList& cellFaces = mesh.cells()[cellIndex];
-    const scalar cellVolume = mesh.cellVolumes()[cellIndex];
-
-    // Find the tet that the point occupies
-    forAll(cellFaces, faceI)
-    {
-        // Decompose each face into triangles, making a tet when
-        // augmented by the cell centre
-        const labelList& facePoints = mesh.faces()[cellFaces[faceI]];
-
-        label pointI = 1;
-        while ((pointI + 1) < facePoints.size())
-        {
-            // Cartesian co-ordinates of the triangle vertices
-            const vector& P1 = mesh.points()[facePoints[0]];
-            const vector& P2 = mesh.points()[facePoints[pointI]];
-            const vector& P3 = mesh.points()[facePoints[pointI + 1]];
-
-            // Edge vectors
-            const vector e1 = P1 - P0;
-            const vector e2 = P2 - P0;
-            const vector e3 = P3 - P0;
-
-            // Solve for interpolation weighting factors
-
-            // Source term
-            const vector rhs = position - P0;
-
-            // Determinant of coefficients matrix
-            // Note: if det(A) = 0 the tet is degenerate
-            const scalar detA =
-                e1.x()*e2.y()*e3.z() + e2.x()*e3.y()*e1.z()
-              + e3.x()*e1.y()*e2.z() - e1.x()*e3.y()*e2.z()
-              - e2.x()*e1.y()*e3.z() - e3.x()*e2.y()*e1.z();
-
-            if (mag(detA/cellVolume) > tol)
-            {
-                // Solve using Cramers' rule
-                const scalar u =
+    List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
     (
-                   rhs.x()*e2.y()*e3.z() + e2.x()*e3.y()*rhs.z()
-                  + e3.x()*rhs.y()*e2.z() - rhs.x()*e3.y()*e2.z()
-                  - e2.x()*rhs.y()*e3.z() - e3.x()*e2.y()*rhs.z()
-                )/detA;
+        mesh,
+        cellI
+    );
 
-                const scalar v =
-                (
-                    e1.x()*rhs.y()*e3.z() + rhs.x()*e3.y()*e1.z()
-                  + e3.x()*e1.y()*rhs.z() - e1.x()*e3.y()*rhs.z()
-                  - rhs.x()*e1.y()*e3.z() - e3.x()*rhs.y()*e1.z()
-                )/detA;
+    const faceList& pFaces = mesh.faces();
+    const scalar cellVolume = mesh.cellVolumes()[cellI];
 
-                const scalar w =
-                (
-                    e1.x()*e2.y()*rhs.z() + e2.x()*rhs.y()*e1.z()
-                  + rhs.x()*e1.y()*e2.z() - e1.x()*rhs.y()*e2.z()
-                  - e2.x()*e1.y()*rhs.z() - rhs.x()*e2.y()*e1.z()
-                )/detA;
+    forAll(cellTets, tetI)
+    {
+        const tetIndices& tetIs = cellTets[tetI];
+
+        const face& f = pFaces[tetIs.face()];
+
+        // Barycentric coordinates of the position
+        scalar det = tetIs.tet(mesh).barycentric(position, weights_);
+
+        if (mag(det/cellVolume) > tol)
+        {
+            const scalar& u = weights_[0];
+            const scalar& v = weights_[1];
+            const scalar& w = weights_[2];
 
-                // Check if point is in tet
-                // value = 0 indicates position lies on a tet face
             if
             (
-                   (u + tol > 0) && (v + tol > 0) && (w + tol > 0)
+                (u + tol > 0)
+             && (v + tol > 0)
+             && (w + tol > 0)
              && (u + v + w < 1 + tol)
             )
             {
-                    faceVertices_[0] = facePoints[0];
-                    faceVertices_[1] = facePoints[pointI];
-                    faceVertices_[2] = facePoints[pointI + 1];
-
-                    weights_[0] = u;
-                    weights_[1] = v;
-                    weights_[2] = w;
-                    weights_[3] = 1.0 - (u + v + w);
+                faceVertices_[0] = f[tetIs.faceBasePt()];
+                faceVertices_[1] = f[tetIs.facePtA()];;
+                faceVertices_[2] = f[tetIs.facePtB()];;
 
                 return;
             }
-                else
-                {
-                    scalar minU = mag(u);
-                    scalar minV = mag(v);
-                    scalar minW = mag(w);
-                    if (minU > 1.0)
-                    {
-                        minU -= 1.0;
-                    }
-                    if (minV > 1.0)
-                    {
-                        minV -= 1.0;
-                    }
-                    if (minW > 1.0)
-                    {
-                        minW -= 1.0;
-                    }
-                    const scalar minUVW = mag(minU + minV + minW);
-
-                    if (minUVW < minUVWClose)
-                    {
-                        minUVWClose = minUVW;
-                        pointIClose = pointI;
-                        faceClose = faceI;
-                    }
         }
     }
 
-            pointI++;
+    // A suitable point in a tetrahedron was not found, find the
+    // nearest.
+
+    scalar minNearDist = VGREAT;
+
+    label nearestTetI = -1;
+
+    forAll(cellTets, tetI)
+    {
+        const tetIndices& tetIs = cellTets[tetI];
+
+        scalar nearDist = tetIs.tet(mesh).nearestPoint(position).distance();
+
+        if (nearDist < minNearDist)
+        {
+            minNearDist = nearDist;
+
+            nearestTetI = tetI;
         }
     }
 
     if (debug)
     {
         Pout<< "cellPointWeight::findTetrahedron" << nl
-            << "    Tetrahedron search failed; using closest tet values to "
-            << "point " << nl << "    cell: " << cellIndex << nl << endl;
+            << "    Tetrahedron search failed; using closest tet to point "
+            << position << nl
+            << "    cell: "
+            << cellI << nl
+            << endl;
     }
 
-    const labelList& facePointsClose = mesh.faces()[cellFaces[faceClose]];
-    faceVertices_[0] = facePointsClose[0];
-    faceVertices_[1] = facePointsClose[pointIClose];
-    faceVertices_[2] = facePointsClose[pointIClose + 1];
 
-    weights_[0] = 0.25;
-    weights_[1] = 0.25;
-    weights_[2] = 0.25;
-    weights_[3] = 0.25;
+    const tetIndices& tetIs = cellTets[nearestTetI];
+
+    const face& f = pFaces[tetIs.face()];
+
+    // Barycentric coordinates of the position, ignoring if the
+    // determinant is suitable.  If not, the return from barycentric
+    // to weights_ is safe.
+    tetIs.tet(mesh).barycentric(position, weights_);
+
+    faceVertices_[0] = f[tetIs.faceBasePt()];
+    faceVertices_[1] = f[tetIs.facePtA()];
+    faceVertices_[2] = f[tetIs.facePtB()];
 }
 
 
@@ -186,119 +142,112 @@ void Foam::cellPointWeight::findTriangle
 (
     const polyMesh& mesh,
     const vector& position,
-    const label faceIndex
+    const label faceI
 )
 {
     if (debug)
     {
         Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
             << "position = " << position << nl
-            << "faceIndex = " << faceIndex << endl;
+            << "faceI = " << faceI << endl;
     }
 
-    // Initialise closest triangle variables
-    scalar minUVClose = VGREAT;
-    label pointIClose = 0;
+    List<tetIndices> faceTets = polyMeshTetDecomposition::faceTetIndices
+    (
+        mesh,
+        mesh.faceOwner()[faceI],
+        faceI
+    );
 
-    // Decompose each face into triangles, making a tet when
-    // augmented by the cell centre
-    const labelList& facePoints = mesh.faces()[faceIndex];
+    const scalar faceAreaSqr = magSqr(mesh.faceAreas()[faceI]);
 
-    const scalar faceArea2 = magSqr(mesh.faceAreas()[faceIndex]);
+    const face& f =  mesh.faces()[faceI];
 
-    label pointI = 1;
-    while ((pointI + 1) < facePoints.size())
+    forAll(faceTets, tetI)
     {
-        // Cartesian co-ordinates of the triangle vertices
-        const vector& P1 = mesh.points()[facePoints[0]];
-        const vector& P2 = mesh.points()[facePoints[pointI]];
-        const vector& P3 = mesh.points()[facePoints[pointI + 1]];
+        const tetIndices& tetIs = faceTets[tetI];
 
-        // Direction vectors
-        vector v1 = position - P1;
-        const vector v2 = P2 - P1;
-        const vector v3 = P3 - P1;
+        List<scalar> triWeights(3);
 
-        // Plane normal
-        vector n = v2 ^ v3;
-        n /= mag(n);
+        // Barycentric coordinates of the position
+        scalar det = tetIs.faceTri(mesh).barycentric(position, triWeights);
 
-        // Remove any offset to plane
-        v1 -= (n & v1)*v1;
-
-        // Helper variables
-        const scalar d12 = v1 & v2;
-        const scalar d13 = v1 & v3;
-        const scalar d22 = v2 & v2;
-        const scalar d23 = v2 & v3;
-        const scalar d33 = v3 & v3;
-
-        // Determinant of coefficients matrix
-        // Note: if det(A) = 0 the triangle is degenerate
-        const scalar detA = d22*d33 - d23*d23;
-
-        if (0.25*detA/faceArea2 > tol)
+        if (0.25*mag(det)/faceAreaSqr > tol)
         {
-            // Solve using Cramers' rule
-            const scalar u = (d12*d33 - d23*d13)/detA;
-            const scalar v = (d22*d13 - d12*d23)/detA;
+            const scalar& u = triWeights[0];
+            const scalar& v = triWeights[1];
 
-            // Check if point is in triangle
-            if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
+            if
+            (
+                (u + tol > 0)
+             && (v + tol > 0)
+             && (u + v < 1 + tol)
+            )
             {
-                // Indices of the cell vertices making up the triangle
-                faceVertices_[0] = facePoints[0];
-                faceVertices_[1] = facePoints[pointI];
-                faceVertices_[2] = facePoints[pointI + 1];
+                // Weight[0] is for the cell centre.
+                weights_[0] = 0;
+                weights_[1] = triWeights[0];
+                weights_[2] = triWeights[1];
+                weights_[3] = triWeights[2];
 
-                weights_[0] = u;
-                weights_[1] = v;
-                weights_[2] = 1.0 - (u + v);
-                weights_[3] = 0.0;
+                faceVertices_[0] = f[tetIs.faceBasePt()];
+                faceVertices_[1] = f[tetIs.facePtA()];;
+                faceVertices_[2] = f[tetIs.facePtB()];;
 
                 return;
             }
-            else
-            {
-                scalar minU = mag(u);
-                scalar minV = mag(v);
-                if (minU > 1.0)
-                {
-                    minU -= 1.0;
-                }
-                if (minV > 1.0)
-                {
-                    minV -= 1.0;
-                }
-                const scalar minUV = mag(minU + minV);
-
-                if (minUV < minUVClose)
-                {
-                    minUVClose = minUV;
-                    pointIClose = pointI;
-                }
         }
     }
 
-        pointI++;
+    // A suitable point in a triangle was not found, find the nearest.
+
+    scalar minNearDist = VGREAT;
+
+    label nearestTetI = -1;
+
+    forAll(faceTets, tetI)
+    {
+        const tetIndices& tetIs = faceTets[tetI];
+
+        scalar nearDist = tetIs.faceTri(mesh).nearestPoint(position).distance();
+
+        if (nearDist < minNearDist)
+        {
+            minNearDist = nearDist;
+
+            nearestTetI = tetI;
+        }
     }
 
     if (debug)
     {
-        Pout<< "Foam::cellPointWeight::findTriangle"
-            << "Triangle search failed; using closest triangle to point" << nl
-            << "    cell face: " << faceIndex << nl << endl;
+        Pout<< "cellPointWeight::findTriangle" << nl
+            << "    Triangle search failed; using closest tri to point "
+            << position << nl
+            << "    face: "
+            << faceI << nl
+            << endl;
     }
 
-    // Indices of the cell vertices making up the triangle
-    faceVertices_[0] = facePoints[0];
-    faceVertices_[1] = facePoints[pointIClose];
-    faceVertices_[2] = facePoints[pointIClose + 1];
+    const tetIndices& tetIs = faceTets[nearestTetI];
 
-    weights_[0] = 1.0/3.0;
-    weights_[1] = 1.0/3.0;
-    weights_[2] = 1.0/3.0;
-    weights_[3] = 0.0;
+    // Barycentric coordinates of the position, ignoring if the
+    // determinant is suitable.  If not, the return from barycentric
+    // to triWeights is safe.
+
+    List<scalar> triWeights(3);
+
+    tetIs.faceTri(mesh).barycentric(position, triWeights);
+
+    // Weight[0] is for the cell centre.
+    weights_[0] = 0;
+    weights_[1] = triWeights[0];
+    weights_[2] = triWeights[1];
+    weights_[3] = triWeights[2];
+
+    faceVertices_[0] = f[tetIs.faceBasePt()];
+    faceVertices_[1] = f[tetIs.facePtA()];
+    faceVertices_[2] = f[tetIs.facePtB()];
 }
 
 
@@ -308,21 +257,23 @@ Foam::cellPointWeight::cellPointWeight
 (
     const polyMesh& mesh,
     const vector& position,
-    const label cellIndex,
-    const label faceIndex
+    const label cellI,
+    const label faceI
 )
 :
-    cellIndex_(cellIndex)
+    cellI_(cellI),
+    weights_(4),
+    faceVertices_(3)
 {
-    if (faceIndex < 0)
+    if (faceI < 0)
     {
         // Face data not supplied
-        findTetrahedron(mesh, position, cellIndex);
+        findTetrahedron(mesh, position, cellI);
     }
     else
     {
         // Face data supplied
-        findTriangle(mesh, position, faceIndex);
+        findTriangle(mesh, position, faceI);
     }
 }
 

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.H

@@ -55,13 +55,13 @@ protected:
     // Protected data
 
        //- Cell index
-       const label cellIndex_;
+       const label cellI_;
 
        //- Weights applied to tet vertices
-       FixedList<scalar, 4> weights_;
+       List<scalar> weights_;
 
        //- Face vertex indices
-       FixedList<label, 3> faceVertices_;
+       List<label> faceVertices_;
 
 
     // Protected Member Functions
@@ -70,14 +70,14 @@ protected:
         (
             const polyMesh& mesh,
             const vector& position,
-            const label cellIndex
+            const label cellI
         );
 
         void findTriangle
         (
             const polyMesh& mesh,
             const vector& position,
-            const label faceIndex
+            const label faceI
         );
 
 
@@ -98,8 +98,8 @@ public:
         (
             const polyMesh& mesh,
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         );
 
 
@@ -108,17 +108,17 @@ public:
         //- Cell index
         inline label cell() const
         {
-            return cellIndex_;
+            return cellI_;
         }
 
         //- interpolation weights
-        inline const FixedList<scalar, 4>& weights() const
+        inline const List<scalar>& weights() const
         {
             return weights_;
         }
 
         //- interpolation addressing for points on face
-        inline const FixedList<label, 3>& faceVertices() const
+        inline const List<label>& faceVertices() const
         {
             return faceVertices_;
         }

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H

@@ -82,8 +82,17 @@ public:
         inline Type interpolate
         (
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
+        ) const;
+
+        //- Interpolate field to the given point in the tetrahedron
+        //  defined by the given indices.
+        inline Type interpolate
+        (
+            const vector& position,
+            const tetIndices& tetIs,
+            const label faceI = -1
         ) const;
 };
 

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H

@@ -31,13 +31,13 @@ inline Type Foam::interpolationCellPoint<Type>::interpolate
     const cellPointWeight& cpw
 ) const
 {
-    const FixedList<scalar, 4>& weights = cpw.weights();
-    const FixedList<label, 3>& faceVertices = cpw.faceVertices();
+    const List<scalar>& weights = cpw.weights();
+    const List<label>& faceVertices = cpw.faceVertices();
 
-    Type t = psip_[faceVertices[0]]*weights[0];
-    t += psip_[faceVertices[1]]*weights[1];
-    t += psip_[faceVertices[2]]*weights[2];
-    t += this->psi_[cpw.cell()]*weights[3];
+    Type t = this->psi_[cpw.cell()]*weights[0];
+    t += psip_[faceVertices[0]]*weights[1];
+    t += psip_[faceVertices[1]]*weights[2];
+    t += psip_[faceVertices[2]]*weights[3];
 
     return t;
 }
@@ -47,11 +47,66 @@ template<class Type>
 inline Type Foam::interpolationCellPoint<Type>::interpolate
 (
     const vector& position,
-    const label celli,
-    const label facei
+    const label cellI,
+    const label faceI
 ) const
 {
-    return interpolate(cellPointWeight(this->pMesh_, position, celli, facei));
+    return interpolate(cellPointWeight(this->pMesh_, position, cellI, faceI));
+}
+
+
+template<class Type>
+inline Type Foam::interpolationCellPoint<Type>::interpolate
+(
+    const vector& position,
+    const tetIndices& tetIs,
+    const label faceI
+) const
+{
+    // Assumes that the position is consistent with the supplied
+    // tetIndices.  Does not pay attention to whether or not faceI is
+    // supplied or not - the result will be essentially the same.
+    // Performs a consistency check, however.
+
+    if (faceI >= 0)
+    {
+        if (faceI != tetIs.face())
+        {
+            FatalErrorIn
+            (
+                "inline Type Foam::interpolationCellPoint<Type>::interpolate"
+                "("
+                    "const vector& position, "
+                    "const tetIndices& tetIs, "
+                    "const label faceI"
+                ") const"
+            )
+                << "specified face " << faceI << " inconsistent with the face "
+                << "stored by tetIndices: " << tetIs.face()
+                << exit(FatalError);
+        }
+    }
+
+    List<scalar> weights;
+
+    tetIs.tet(this->pMesh_).barycentric(position, weights);
+
+    const faceList& pFaces = this->pMesh_.faces();
+
+    const face& f = pFaces[tetIs.face()];
+
+    // Order of weights is the same as that of the vertices of the tet, i.e.
+    // cellCentre, faceBasePt, facePtA, facePtB.
+
+    Type t = this->psi_[tetIs.cell()]*weights[0];
+
+    t += psip_[f[tetIs.faceBasePt()]]*weights[1];
+
+    t += psip_[f[tetIs.facePtA()]]*weights[2];
+
+    t += psip_[f[tetIs.facePtB()]]*weights[3];
+
+    return t;
 }
 
 

src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C

@@ -56,8 +56,8 @@ template<class Type>
 Type interpolationCellPointFace<Type>::interpolate
 (
     const vector& position,
-    const label nCell,
-    const label facei
+    const label cellI,
+    const label faceI
 ) const
 {
     Type ts[4];
@@ -68,10 +68,10 @@ Type interpolationCellPointFace<Type>::interpolate
     Type t = pTraits<Type>::zero;
 
     // only use face information when the position is on a face
-    if (facei < 0)
+    if (faceI < 0)
     {
-        const vector& cellCentre = this->pMesh_.cellCentres()[nCell];
-        const labelList& cellFaces = this->pMesh_.cells()[nCell];
+        const vector& cellCentre = this->pMesh_.cellCentres()[cellI];
+        const labelList& cellFaces = this->pMesh_.cells()[cellI];
 
         vector projection = position - cellCentre;
         tetPoints[3] = cellCentre;
@@ -85,9 +85,9 @@ Type interpolationCellPointFace<Type>::interpolate
         label closestFace = -1;
         scalar minDistance = GREAT;
 
-        forAll(cellFaces, facei)
+        forAll(cellFaces, faceI)
         {
-            label nFace = cellFaces[facei];
+            label nFace = cellFaces[faceI];
 
             vector normal = this->pMeshFaceAreas_[nFace];
             normal /= mag(normal);
@@ -160,10 +160,10 @@ Type interpolationCellPointFace<Type>::interpolate
         {
             minDistance = GREAT;
 
-            label facei = 0;
-            while (facei < cellFaces.size() && !foundTet)
+            label faceI = 0;
+            while (faceI < cellFaces.size() && !foundTet)
             {
-                label nFace = cellFaces[facei];
+                label nFace = cellFaces[faceI];
                 if (nFace < this->pMeshFaceAreas_.size())
                 {
                     foundTet = findTet
@@ -179,7 +179,7 @@ Type interpolationCellPointFace<Type>::interpolate
                         minDistance
                     );
                 }
-                facei++;
+                faceI++;
             }
         }
 
@@ -217,16 +217,16 @@ Type interpolationCellPointFace<Type>::interpolate
             }
             else
             {
-                label patchi =
+                label patchI =
                     this->pMesh_.boundaryMesh().whichPatch(closestFace);
 
                 // If the boundary patch is not empty use the face value
                 // else use the cell value
-                if (this->psi_.boundaryField()[patchi].size())
+                if (this->psi_.boundaryField()[patchI].size())
                 {
-                    ts[2] = this->psi_.boundaryField()[patchi]
+                    ts[2] = this->psi_.boundaryField()[patchI]
                     [
-                        this->pMesh_.boundaryMesh()[patchi].whichFace
+                        this->pMesh_.boundaryMesh()[patchI].whichFace
                         (
                             closestFace
                         )
@@ -234,11 +234,11 @@ Type interpolationCellPointFace<Type>::interpolate
                 }
                 else
                 {
-                    ts[2] = this->psi_[nCell];
+                    ts[2] = this->psi_[cellI];
                 }
             }
 
-            ts[3] = this->psi_[nCell];
+            ts[3] = this->psi_[cellI];
 
             for (label n=0; n<4; n++)
             {
@@ -251,9 +251,9 @@ Type interpolationCellPointFace<Type>::interpolate
         else
         {
             Info<< "interpolationCellPointFace<Type>::interpolate"
-                << "(const vector&, const label nCell) const : "
+                << "(const vector&, const label cellI) const : "
                 << "search failed; using closest cellFace value" << endl
-                << "cell number " << nCell << tab
+                << "cell number " << cellI << tab
                 << "position " << position << endl;
 
             if (closestFace < psis_.size())
@@ -262,16 +262,16 @@ Type interpolationCellPointFace<Type>::interpolate
             }
             else
             {
-                label patchi =
+                label patchI =
                     this->pMesh_.boundaryMesh().whichPatch(closestFace);
 
                 // If the boundary patch is not empty use the face value
                 // else use the cell value
-                if (this->psi_.boundaryField()[patchi].size())
+                if (this->psi_.boundaryField()[patchI].size())
                 {
-                    t = this->psi_.boundaryField()[patchi]
+                    t = this->psi_.boundaryField()[patchI]
                     [
-                        this->pMesh_.boundaryMesh()[patchi].whichFace
+                        this->pMesh_.boundaryMesh()[patchI].whichFace
                         (
                             closestFace
                         )
@@ -279,7 +279,7 @@ Type interpolationCellPointFace<Type>::interpolate
                 }
                 else
                 {
-                    t = this->psi_[nCell];
+                    t = this->psi_[cellI];
                 }
             }
         }
@@ -289,7 +289,7 @@ Type interpolationCellPointFace<Type>::interpolate
         bool foundTriangle = findTriangle
         (
             position,
-            facei,
+            faceI,
             tetPointLabels,
             phi
         );
@@ -304,48 +304,48 @@ Type interpolationCellPointFace<Type>::interpolate
             }
 
             // ... and the face value
-            if (facei < psis_.size())
+            if (faceI < psis_.size())
             {
-                t += phi[2]*psis_[facei];
+                t += phi[2]*psis_[faceI];
             }
             else
             {
-                label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
+                label patchI = this->pMesh_.boundaryMesh().whichPatch(faceI);
 
                 // If the boundary patch is not empty use the face value
                 // else use the cell value
-                if (this->psi_.boundaryField()[patchi].size())
+                if (this->psi_.boundaryField()[patchI].size())
                 {
-                    t += phi[2]*this->psi_.boundaryField()[patchi]
-                        [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
+                    t += phi[2]*this->psi_.boundaryField()[patchI]
+                        [this->pMesh_.boundaryMesh()[patchI].whichFace(faceI)];
                 }
                 else
                 {
-                    t += phi[2]*this->psi_[nCell];
+                    t += phi[2]*this->psi_[cellI];
                 }
             }
         }
         else
         {
             // use face value only
-            if (facei < psis_.size())
+            if (faceI < psis_.size())
             {
-                t = psis_[facei];
+                t = psis_[faceI];
             }
             else
             {
-                label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
+                label patchI = this->pMesh_.boundaryMesh().whichPatch(faceI);
 
                 // If the boundary patch is not empty use the face value
                 // else use the cell value
-                if (this->psi_.boundaryField()[patchi].size())
+                if (this->psi_.boundaryField()[patchI].size())
                 {
-                    t = this->psi_.boundaryField()[patchi]
-                        [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
+                    t = this->psi_.boundaryField()[patchI]
+                        [this->pMesh_.boundaryMesh()[patchI].whichFace(faceI)];
                 }
                 else
                 {
-                    t = this->psi_[nCell];
+                    t = this->psi_[cellI];
                 }
             }
         }

src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.H

@@ -99,8 +99,8 @@ public:
         Type interpolate
         (
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         ) const;
 };
 

src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.C

@@ -24,9 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cellPointWeightWallModified.H"
-#include "wallPolyPatch.H"
-#include "polyMesh.H"
-#include "polyBoundaryMesh.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -34,36 +31,30 @@ Foam::cellPointWeightWallModified::cellPointWeightWallModified
 (
     const polyMesh& mesh,
     const vector& position,
-    const label cellIndex,
-    const label faceIndex
+    const label cellI,
+    const label faceI
 )
 :
-    cellPointWeight(mesh, position, cellIndex, faceIndex)
+    cellPointWeight(mesh, position, cellI, faceI)
 {
-    if (faceIndex < 0)
-    {
-        findTetrahedron(mesh, position, cellIndex);
-    }
-    else
+    // findTetrahedron or findTriangle will already have been called
+    // by the cellPointWeight constructor
+
+    if (faceI >= 0)
     {
         const polyBoundaryMesh& bm = mesh.boundaryMesh();
-        label patchI = bm.whichPatch(faceIndex);
+        label patchI = bm.whichPatch(faceI);
         if (patchI != -1)
         {
             if (isA<wallPolyPatch>(bm[patchI]))
             {
                 // Apply cell centre value wall faces
-                weights_[0] = 0.0;
+                weights_[0] = 1.0;
                 weights_[1] = 0.0;
                 weights_[2] = 0.0;
-                weights_[3] = 1.0;
+                weights_[3] = 0.0;
             }
         }
-        else
-        {
-            // Interpolate
-            findTriangle(mesh, position, faceIndex);
-        }
     }
 }
 

src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.H

@@ -36,6 +36,9 @@ SourceFiles
 #define cellPointWeightWallModified_H
 
 #include "cellPointWeight.H"
+#include "wallPolyPatch.H"
+#include "polyMesh.H"
+#include "polyBoundaryMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,8 +64,8 @@ public:
         (
             const polyMesh& mesh,
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         );
 };
 

src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H

@@ -74,8 +74,17 @@ public:
         inline Type interpolate
         (
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
+        ) const;
+
+        //- Interpolate field to the given point in the tetrahedron
+        //  defined by the given indices.
+        inline Type interpolate
+        (
+            const vector& position,
+            const tetIndices& tetIs,
+            const label faceI = -1
         ) const;
 };
 

src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModifiedI.H

@@ -31,13 +31,13 @@ inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
     const cellPointWeightWallModified& cpw
 ) const
 {
-    const FixedList<scalar, 4>& weights = cpw.weights();
-    const FixedList<label, 3>& faceVertices = cpw.faceVertices();
+    const List<scalar>& weights = cpw.weights();
+    const List<label>& faceVertices = cpw.faceVertices();
 
-    Type t = this->psip_[faceVertices[0]]*weights[0];
-    t += this->psip_[faceVertices[1]]*weights[1];
-    t += this->psip_[faceVertices[2]]*weights[2];
-    t += this->psi_[cpw.cell()]*weights[3];
+    Type t = this->psi_[cpw.cell()]*weights[0];
+    t += this->psip_[faceVertices[0]]*weights[1];
+    t += this->psip_[faceVertices[1]]*weights[2];
+    t += this->psip_[faceVertices[2]]*weights[3];
 
     return t;
 }
@@ -47,22 +47,74 @@ template<class Type>
 inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
 (
     const vector& position,
-    const label celli,
-    const label facei
+    const label cellI,
+    const label faceI
 ) const
 {
-    return
-        interpolate
+    return interpolate
     (
         cellPointWeightWallModified
         (
             this->pMesh_,
             position,
-                celli,
-                facei
+            cellI,
+            faceI
         )
     );
 }
 
 
+template<class Type>
+inline Type Foam::interpolationCellPointWallModified<Type>::interpolate
+(
+    const vector& position,
+    const tetIndices& tetIs,
+    const label faceI
+) const
+{
+    if (faceI >= 0)
+    {
+        if (faceI != tetIs.face())
+        {
+            FatalErrorIn
+            (
+                "inline Type "
+                "Foam::interpolationCellPointWallModifie<Type>::interpolate"
+                "("
+                    "const vector& position, "
+                    "const tetIndices& tetIs, "
+                    "const label faceI"
+                ") const"
+            )
+                << "specified face " << faceI << " inconsistent with the face "
+                << "stored by tetIndices: " << tetIs.face()
+                << exit(FatalError);
+        }
+
+        const polyBoundaryMesh& bm = this->pMesh_.boundaryMesh();
+        label patchI = bm.whichPatch(faceI);
+
+        if (patchI != -1)
+        {
+            if (isA<wallPolyPatch>(bm[patchI]))
+            {
+                Type t = this->psi_[tetIs.cell()];
+
+                return t;
+            }
+        }
+    }
+
+    // If the wall face selection did not return, then use the normal
+    // interpolate method
+
+    return interpolationCellPoint<Type>::interpolate
+    (
+        position,
+        tetIs,
+        faceI
+    );
+}
+
+
 // ************************************************************************* //

src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.H

@@ -81,8 +81,8 @@ public:
         inline Type interpolate
         (
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         ) const;
 };
 

src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPointI.H

@@ -39,13 +39,13 @@ template<class Type>
 inline Type Foam::interpolationPoint<Type>::interpolate
 (
     const vector& position,
-    const label celli,
-    const label facei
+    const label cellI,
+    const label faceI
 ) const
 {
     return interpolate
     (
-        pointMVCWeight(this->pMesh_, position, celli, facei)
+        pointMVCWeight(this->pMesh_, position, cellI, faceI)
     );
 }
 

src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.H

@@ -119,8 +119,8 @@ public:
         (
             const polyMesh& mesh,
             const vector& position,
-            const label nCell,
-            const label facei = -1
+            const label cellI,
+            const label faceI = -1
         );
 
 

src/lagrangian/basic/Cloud/Cloud.C

@@ -30,6 +30,41 @@ License
 #include "mapPolyMesh.H"
 #include "Time.H"
 #include "OFstream.H"
+#include "wallPolyPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParticleType>
+const Foam::scalar Foam::Cloud<ParticleType>::trackingCorrectionTol = 1e-5;
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::calcCellWallFaces() const
+{
+    cellWallFacesPtr_.reset(new PackedBoolList(pMesh().nCells(), false));
+
+    PackedBoolList& cellWallFaces = cellWallFacesPtr_();
+
+    const polyBoundaryMesh& patches = pMesh().boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            const polyPatch& patch = patches[patchI];
+
+            const labelList& pFaceCells = patch.faceCells();
+
+            forAll(pFaceCells, pFCI)
+            {
+                cellWallFaces[pFaceCells[pFCI]] = true;
+            }
+        }
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -43,7 +78,11 @@ Foam::Cloud<ParticleType>::Cloud
     cloud(pMesh),
     IDLList<ParticleType>(),
     polyMesh_(pMesh),
-    particleCount_(0)
+    particleCount_(0),
+    labels_(),
+    cellTree_(),
+    nTrackingRescues_(),
+    cellWallFacesPtr_()
 {
     IDLList<ParticleType>::operator=(particles);
 }
@@ -60,7 +99,11 @@ Foam::Cloud<ParticleType>::Cloud
     cloud(pMesh, cloudName),
     IDLList<ParticleType>(),
     polyMesh_(pMesh),
-    particleCount_(0)
+    particleCount_(0),
+    labels_(),
+    cellTree_(),
+    nTrackingRescues_(),
+    cellWallFacesPtr_()
 {
     IDLList<ParticleType>::operator=(particles);
 }
@@ -68,6 +111,250 @@ Foam::Cloud<ParticleType>::Cloud
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::findCellFacePt
+(
+    const point& pt,
+    label& cellI,
+    label& tetFaceI,
+    label& tetPtI
+) const
+{
+    cellI = -1;
+    tetFaceI = -1;
+    tetPtI = -1;
+
+    const indexedOctree<treeDataCell>& tree = cellTree();
+
+    // Find nearest cell to the point
+
+    pointIndexHit info = tree.findNearest(pt, sqr(GREAT));
+
+    if (info.hit())
+    {
+        label nearestCellI = tree.shapes().cellLabels()[info.index()];
+
+        // Check the nearest cell to see if the point is inside.
+        findFacePt(nearestCellI, pt, tetFaceI, tetPtI);
+
+        if (tetFaceI != -1)
+        {
+            // Point was in the nearest cell
+
+            cellI = nearestCellI;
+
+            return;
+        }
+        else
+        {
+            // Check the other possible cells that the point may be in
+
+            labelList testCells = tree.findIndices(pt);
+
+            forAll(testCells, pCI)
+            {
+                label testCellI = tree.shapes().cellLabels()[testCells[pCI]];
+
+                if (testCellI == nearestCellI)
+                {
+                    // Don't retest the nearest cell
+
+                    continue;
+                }
+
+                // Check the test cell to see if the point is inside.
+                findFacePt(testCellI, pt, tetFaceI, tetPtI);
+
+                if (tetFaceI != -1)
+                {
+                    // Point was in the test cell
+
+                    cellI = testCellI;
+
+                    return;
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void Foam::Cloud<ParticleType>::findCellFacePt"
+            "("
+                "const point& pt, "
+                "label& cellI, "
+                "label& tetFaceI, "
+                "label& tetPtI"
+            ") const"
+        )   << "Did not find nearest cell in search tree."
+            << abort(FatalError);
+    }
+}
+
+
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::findFacePt
+(
+    label cellI,
+    const point& pt,
+    label& tetFaceI,
+    label& tetPtI
+) const
+{
+    tetFaceI = -1;
+    tetPtI = -1;
+
+    List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
+    (
+        polyMesh_,
+        cellI
+    );
+
+    forAll(cellTets, tetI)
+    {
+        const tetIndices& cellTetIs = cellTets[tetI];
+
+        if (inTet(pt, cellTetIs.tet(polyMesh_)))
+        {
+            tetFaceI = cellTetIs.face();
+            tetPtI = cellTetIs.tetPt();
+
+            return;
+        }
+    }
+}
+
+
+template<class ParticleType>
+bool Foam::Cloud<ParticleType>::inTet
+(
+    const point& pt,
+    const tetPointRef& tet
+) const
+{
+    // For robustness, assuming that the point is in the tet unless
+    // "definitively" shown otherwise by obtaining a positive dot
+    // product greater than a tolerance of SMALL.
+
+    // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal
+    // vectors and base points for the half-space planes are:
+    // area[0] = tet.Sa();
+    // area[1] = tet.Sb();
+    // area[2] = tet.Sc();
+    // area[3] = tet.Sd();
+    // planeBase[0] = tetBasePt = tet.b()
+    // planeBase[1] = ptA       = tet.c()
+    // planeBase[2] = tetBasePt = tet.b()
+    // planeBase[3] = tetBasePt = tet.b()
+
+    vector n = vector::zero;
+
+    {
+        // 0, a
+        const point& basePt = tet.b();
+
+        n = tet.Sa();
+        n /= (mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
+        {
+            return false;
+        }
+    }
+
+    {
+        // 1, b
+        const point& basePt = tet.c();
+
+        n = tet.Sb();
+        n /= (mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
+        {
+            return false;
+        }
+    }
+
+    {
+        // 2, c
+        const point& basePt = tet.b();
+
+        n = tet.Sc();
+        n /= (mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
+        {
+            return false;
+        }
+    }
+
+    {
+        // 3, d
+        const point& basePt = tet.b();
+
+        n = tet.Sd();
+        n /= (mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+template<class ParticleType>
+const Foam::indexedOctree<Foam::treeDataCell>&
+Foam::Cloud<ParticleType>::cellTree() const
+{
+    if (cellTree_.empty())
+    {
+        treeBoundBox overallBb(polyMesh_.points());
+
+        Random rndGen(261782);
+
+        overallBb = overallBb.extend(rndGen, 1E-4);
+        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        cellTree_.reset
+        (
+            new indexedOctree<treeDataCell>
+            (
+                treeDataCell
+                (
+                    false,      // not cache bb
+                    polyMesh_
+                ),
+                overallBb,
+                8,              // maxLevel
+                10,             // leafsize
+                3.0             // duplicity
+            )
+        );
+    }
+
+    return cellTree_();
+}
+
+
+template<class ParticleType>
+const Foam::PackedBoolList& Foam::Cloud<ParticleType>::cellHasWallFaces()
+const
+{
+    if (!cellWallFacesPtr_.valid())
+    {
+        calcCellWallFaces();
+    }
+
+    return cellWallFacesPtr_();
+}
+
+
+
 template<class ParticleType>
 Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
 {
@@ -131,6 +418,9 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
         pIter().stepFraction() = 0;
     }
 
+    // Reset nTrackingRescues
+    nTrackingRescues_ = 0;
+
     // While there are particles to transfer
     while (true)
     {
@@ -162,29 +452,29 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
             {
                 // If we are running in parallel and the particle is on a
                 // boundary face
-                if (Pstream::parRun() && p.facei_ >= pMesh().nInternalFaces())
+                if (Pstream::parRun() && p.faceI_ >= pMesh().nInternalFaces())
                 {
-                    label patchi = pbm.whichPatch(p.facei_);
+                    label patchI = pbm.whichPatch(p.faceI_);
 
                     // ... and the face is on a processor patch
                     // prepare it for transfer
-                    if (procPatchIndices[patchi] != -1)
+                    if (procPatchIndices[patchI] != -1)
                     {
                         label n = neighbourProcIndices
                         [
                             refCast<const processorPolyPatch>
                             (
-                                pbm[patchi]
+                                pbm[patchI]
                             ).neighbProcNo()
                         ];
 
-                        p.prepareForParallelTransfer(patchi, td);
+                        p.prepareForParallelTransfer(patchI, td);
 
                         particleTransferLists[n].append(this->remove(&p));
 
                         patchIndexTransferLists[n].append
                         (
-                            procPatchNeighbours[patchi]
+                            procPatchNeighbours[patchI]
                         );
                     }
                 }
@@ -270,15 +560,22 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
                 {
                     ParticleType& newp = newpIter();
 
-                    label patchi = procPatches[receivePatchIndex[pI++]];
+                    label patchI = procPatches[receivePatchIndex[pI++]];
 
-                    newp.correctAfterParallelTransfer(patchi, td);
+                    newp.correctAfterParallelTransfer(patchI, td);
 
                     addParticle(newParticles.remove(&newp));
                 }
             }
         }
     }
+
+    reduce(nTrackingRescues_, sumOp<label>());
+
+    if (nTrackingRescues_ > 0)
+    {
+        Info<< nTrackingRescues_ << " tracking rescue corrections" << endl;
+    }
 }
 
 
@@ -294,24 +591,30 @@ void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
     const labelList& reverseCellMap = mapper.reverseCellMap();
     const labelList& reverseFaceMap = mapper.reverseFaceMap();
 
+    // Reset stored data that relies on the mesh
+    cellTree_.clear();
+    cellWallFacesPtr_.clear();
+
     forAllIter(typename Cloud<ParticleType>, *this, pIter)
     {
-        if (reverseCellMap[pIter().celli_] >= 0)
+        if (reverseCellMap[pIter().cellI_] >= 0)
         {
-            pIter().celli_ = reverseCellMap[pIter().celli_];
+            pIter().cellI_ = reverseCellMap[pIter().cellI_];
 
-            if (pIter().facei_ >= 0 && reverseFaceMap[pIter().facei_] >= 0)
+            if (pIter().faceI_ >= 0 && reverseFaceMap[pIter().faceI_] >= 0)
             {
-                pIter().facei_ = reverseFaceMap[pIter().facei_];
+                pIter().faceI_ = reverseFaceMap[pIter().faceI_];
             }
             else
             {
-                pIter().facei_ = -1;
+                pIter().faceI_ = -1;
             }
+
+            pIter().initCellFacePt();
         }
         else
         {
-            label trackStartCell = mapper.mergedCell(pIter().celli_);
+            label trackStartCell = mapper.mergedCell(pIter().cellI_);
 
             if (trackStartCell < 0)
             {
@@ -319,9 +622,14 @@ void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
             }
 
             vector p = pIter().position();
+
             const_cast<vector&>(pIter().position()) =
                 polyMesh_.cellCentres()[trackStartCell];
+
             pIter().stepFraction() = 0;
+
+            pIter().initCellFacePt();
+
             pIter().track(p);
         }
     }

src/lagrangian/basic/Cloud/Cloud.H

@@ -40,6 +40,11 @@ SourceFiles
 #include "IOField.H"
 #include "IOFieldField.H"
 #include "polyMesh.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "tetPointRef.H"
+#include "polyMeshTetDecomposition.H"
+#include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,15 +83,28 @@ class Cloud
         //- Overall count of particles ever created. Never decreases.
         mutable label particleCount_;
 
-        //- Temporary storage for addressing. Used in findFaces.
+        //- Temporary storage for addressing. Used in findTris.
         mutable DynamicList<label> labels_;
 
+        //- Search tree to allow spatial tet searching
+        mutable autoPtr<indexedOctree<treeDataCell> > cellTree_;
+
+        //- Count of how many tracking rescue corrections have been
+        //  applied
+        mutable label nTrackingRescues_;
+
+        //- Does the cell have wall faces
+        mutable autoPtr<PackedBoolList> cellWallFacesPtr_;
+
 
     // Private Member Functions
 
         //- Initialise cloud on IO constructor
         void initCloud(const bool checkClass);
 
+        //- Find all cells which have wall faces
+        void calcCellWallFaces() const;
+
         //- Read cloud properties dictionary
         void readCloudUniformProperties();
 
@@ -115,6 +133,10 @@ public:
         //- Name of cloud properties dictionary
         static word cloudPropertiesName;
 
+        //- Fraction of distance to tet centre to move a particle to
+        // 'rescue' it from a tracking problem
+        static const scalar trackingCorrectionTol;
+
 
     // Constructors
 
@@ -163,27 +185,27 @@ public:
             }
 
             //- Is this global face an internal face?
-            bool internalFace(const label facei) const
+            bool internalFace(const label faceI) const
             {
-                return polyMesh_.isInternalFace(facei);
+                return polyMesh_.isInternalFace(faceI);
             }
 
             //- Is this global face a boundary face?
-            bool boundaryFace(const label facei) const
+            bool boundaryFace(const label faceI) const
             {
-                return !internalFace(facei);
+                return !internalFace(faceI);
             }
 
             //- Which patch is this global face on
-            label facePatch(const label facei) const
+            label facePatch(const label faceI) const
             {
-                return polyMesh_.boundaryMesh().whichPatch(facei);
+                return polyMesh_.boundaryMesh().whichPatch(faceI);
             }
 
             //- Which face of this patch is this global face
-            label patchFace(const label patchi, const label facei) const
+            label patchFace(const label patchI, const label faceI) const
             {
-                return polyMesh_.boundaryMesh()[patchi].whichFace(facei);
+                return polyMesh_.boundaryMesh()[patchI].whichFace(faceI);
             }
 
             label size() const
@@ -191,6 +213,61 @@ public:
                 return IDLList<ParticleType>::size();
             };
 
+            //- Find the cell, tetFaceI and tetPtI for the given
+            //  position
+            void findCellFacePt
+            (
+                const point& pt,
+                label& cellI,
+                label& tetFaceI,
+                label& tetPtI
+            ) const;
+
+            //- Find the tetFaceI and tetPtI for the given position in
+            //  the supplied cell, tetFaceI and tetPtI = -1 if not
+            //  found
+            void findFacePt
+            (
+                label cellI,
+                const point& pt,
+                label& tetFaceI,
+                label& tetPtI
+            ) const;
+
+            //- Test if the given position is inside the give tet
+            bool inTet
+            (
+                const point& pt,
+                const tetPointRef& tet
+            ) const;
+
+            //- Build (if necessary) and return the cell search tree
+            const indexedOctree<treeDataCell>& cellTree() const;
+
+            //- Return nTrackingRescues
+            label nTrackingRescues() const
+            {
+                return nTrackingRescues_;
+            }
+
+            //- Increment the nTrackingRescues counter
+            void trackingRescue() const
+            {
+                nTrackingRescues_++;
+            }
+
+            //- Whether each cell has any wall faces (demand driven data)
+            const PackedBoolList& cellHasWallFaces() const;
+
+            //- Switch to specify if particles of the cloud can return
+            //  non-zero wall distance values.  By default, assume
+            //  that they can't (default for wallImpactDistance in
+            //  Particle is 0.0).
+            virtual bool hasWallImpactDistance() const
+            {
+                return false;
+            }
+
 
             // Iterators
 

src/lagrangian/basic/Cloud/CloudIO.C

@@ -126,6 +126,13 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
             << "    " << ioP.path() << nl
             << "    assuming the initial cloud contains 0 particles." << endl;
     }
+
+    forAllIter(typename Cloud<ParticleType>, *this, pIter)
+    {
+        ParticleType& p = pIter();
+
+        p.initCellFacePt();
+    }
 }
 
 
@@ -140,7 +147,11 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh),
     polyMesh_(pMesh),
-    particleCount_(0)
+    particleCount_(0),
+    labels_(),
+    cellTree_(),
+    nTrackingRescues_(),
+    cellWallFacesPtr_()
 {
     initCloud(checkClass);
 }
@@ -156,7 +167,11 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh, cloudName),
     polyMesh_(pMesh),
-    particleCount_(0)
+    particleCount_(0),
+    labels_(),
+    cellTree_(),
+    nTrackingRescues_(),
+    cellWallFacesPtr_()
 {
     initCloud(checkClass);
 }

src/lagrangian/basic/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lmeshTools

src/lagrangian/basic/Particle/Particle.C

@@ -29,73 +29,20 @@ License
 #include "symmetryPolyPatch.H"
 #include "cyclicPolyPatch.H"
 #include "processorPolyPatch.H"
-#include "wallPolyPatch.H"
 #include "transform.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class ParticleType>
-void Foam::Particle<ParticleType>::findFaces
-(
-    const vector& position,
-    DynamicList<label>& faceList
-) const
-{
-    const polyMesh& mesh = cloud_.polyMesh_;
-    const labelList& faces = mesh.cells()[celli_];
-    const vector& C = mesh.cellCentres()[celli_];
-
-    faceList.clear();
-    forAll(faces, i)
-    {
-        label facei = faces[i];
-        scalar lam = lambda(C, position, facei);
-
-        if ((lam > 0) && (lam < 1.0))
-        {
-            faceList.append(facei);
-        }
-    }
-}
-
-
-template<class ParticleType>
-void Foam::Particle<ParticleType>::findFaces
-(
-    const vector& position,
-    const label celli,
-    const scalar stepFraction,
-    DynamicList<label>& faceList
-) const
-{
-    const polyMesh& mesh = cloud_.pMesh();
-    const labelList& faces = mesh.cells()[celli];
-    const vector& C = mesh.cellCentres()[celli];
-
-    faceList.clear();
-    forAll(faces, i)
-    {
-        label facei = faces[i];
-        scalar lam = lambda(C, position, facei, stepFraction);
-
-        if ((lam > 0) && (lam < 1.0))
-        {
-            faceList.append(facei);
-        }
-    }
-}
-
-
 template<class ParticleType>
 template<class TrackData>
 void Foam::Particle<ParticleType>::prepareForParallelTransfer
 (
-    const label patchi,
+    const label patchI,
     TrackData& td
 )
 {
     // Convert the face index to be local to the processor patch
-    facei_ = patchFace(patchi, facei_);
+    faceI_ = patchFace(patchI, faceI_);
 }
 
 
@@ -103,15 +50,15 @@ template<class ParticleType>
 template<class TrackData>
 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
 (
-    const label patchi,
+    const label patchI,
     TrackData& td
 )
 {
     const processorPolyPatch& ppp =
         refCast<const processorPolyPatch>
-        (cloud_.pMesh().boundaryMesh()[patchi]);
+        (cloud_.pMesh().boundaryMesh()[patchI]);
 
-    celli_ = ppp.faceCells()[facei_];
+    cellI_ = ppp.faceCells()[faceI_];
 
     if (!ppp.parallel())
     {
@@ -123,7 +70,7 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
         }
         else
         {
-            const tensor& T = ppp.forwardT()[facei_];
+            const tensor& T = ppp.forwardT()[faceI_];
             transformPosition(T);
             static_cast<ParticleType&>(*this).transformProperties(T);
         }
@@ -140,23 +87,48 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
         }
         else
         {
-            position_ -= ppp.separation()[facei_];
+            position_ -= ppp.separation()[faceI_];
             static_cast<ParticleType&>(*this).transformProperties
             (
-                -ppp.separation()[facei_]
+                -ppp.separation()[faceI_]
             );
         }
     }
 
+    tetFaceI_ = faceI_ + ppp.start();
+
+    // Faces either side of a coupled patch have matched base indices,
+    // tetPtI is specified relative to the base point, already and
+    // opposite circulation directions by design, so if the vertices
+    // are:
+    // source:
+    // face    (a b c d e f)
+    // fPtI     0 1 2 3 4 5
+    //            +
+    // destination:
+    // face    (a f e d c b)
+    // fPtI     0 1 2 3 4 5
+    //                  +
+    // where a is the base point of the face are matching , and we
+    // have fPtI = 1 on the source processor face, i.e. vertex b, then
+    // this because of the face circulation direction change, vertex c
+    // is the characterising point on the destination processor face,
+    // giving the destination fPtI as:
+    //     fPtI_d = f.size() - 1 - fPtI_s = 6 - 1 - 1 = 4
+    // This relationship can be verified for other points and sizes of
+    // face.
+
+    tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
+
     // Reset the face index for the next tracking operation
     if (stepFraction_ > (1.0 - SMALL))
     {
         stepFraction_ = 1.0;
-        facei_ = -1;
+        faceI_ = -1;
     }
     else
     {
-        facei_ += ppp.start();
+        faceI_ += ppp.start();
     }
 }
 
@@ -168,27 +140,55 @@ Foam::Particle<ParticleType>::Particle
 (
     const Cloud<ParticleType>& cloud,
     const vector& position,
-    const label celli
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI
 )
 :
     cloud_(cloud),
     position_(position),
-    celli_(celli),
-    facei_(-1),
+    cellI_(cellI),
+    faceI_(-1),
     stepFraction_(0.0),
+    tetFaceI_(tetFaceI),
+    tetPtI_(tetPtI),
     origProc_(Pstream::myProcNo()),
     origId_(cloud_.getNewParticleID())
 {}
 
 
+template<class ParticleType>
+Foam::Particle<ParticleType>::Particle
+(
+    const Cloud<ParticleType>& cloud,
+    const vector& position,
+    const label cellI
+)
+:
+    cloud_(cloud),
+    position_(position),
+    cellI_(cellI),
+    faceI_(-1),
+    stepFraction_(0.0),
+    tetFaceI_(-1),
+    tetPtI_(-1),
+    origProc_(Pstream::myProcNo()),
+    origId_(cloud_.getNewParticleID())
+{
+    initCellFacePt();
+}
+
+
 template<class ParticleType>
 Foam::Particle<ParticleType>::Particle(const Particle<ParticleType>& p)
 :
     cloud_(p.cloud_),
     position_(p.position_),
-    celli_(p.celli_),
-    facei_(p.facei_),
+    cellI_(p.cellI_),
+    faceI_(p.faceI_),
     stepFraction_(p.stepFraction_),
+    tetFaceI_(p.tetFaceI_),
+    tetPtI_(p.tetPtI_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -204,7 +204,7 @@ Foam::label Foam::Particle<ParticleType>::track
     TrackData& td
 )
 {
-    facei_ = -1;
+    faceI_ = -1;
 
     // Tracks to endPosition or stop on boundary
     while (!onBoundary() && stepFraction_ < 1.0 - SMALL)
@@ -212,7 +212,7 @@ Foam::label Foam::Particle<ParticleType>::track
         stepFraction_ += trackToFace(endPosition, td)*(1.0 - stepFraction_);
     }
 
-    return facei_;
+    return faceI_;
 }
 
 
@@ -224,6 +224,7 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
     return track(endPosition, dummyTd);
 }
 
+
 template<class ParticleType>
 template<class TrackData>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
@@ -234,85 +235,263 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 {
     const polyMesh& mesh = cloud_.polyMesh_;
 
-    DynamicList<label>& faces = cloud_.labels_;
-    findFaces(endPosition, faces);
+    const faceList& pFaces = mesh.faces();
+    const pointField& pPts = mesh.points();
+    const vectorField& pC = mesh.cellCentres();
+
+    faceI_ = -1;
 
-    facei_ = -1;
     scalar trackFraction = 0.0;
 
-    if (faces.empty())  // inside cell
-    {
-        trackFraction = 1.0;
-        position_ = endPosition;
-    }
-    else // hit face
-    {
-        scalar lambdaMin = GREAT;
+    // Minimum tetrahedron decomposition of each cell of the mesh into
+    // using the cell centre, base point on face, and further two
+    // points on the face.  For each face of n points, there are n - 2
+    // tets generated.
 
-        if (faces.size() == 1)
+    // The points for each tet are organised to match those used in the
+    // tetrahedron class, supplying them in the order:
+    //     Cc, basePt, pA, pB
+    // where:
+    //   + Cc is the cell centre;
+    //   + basePt is the base point on the face;
+    //   + pA and pB are the remaining points on the face, such that
+    //     the circulation, {basePt, pA, pB} produces a positive
+    //     normal by the right-hand rule.  pA and pB are chosen from
+    //     tetPtI_ do accomplish this depending if the cell owns the
+    //     face, tetPtI_ is the vertex that characterises the tet, and
+    //     is the first vertex on the tet when circulating around the
+    //     face. Therefore, the same tetPtI represents the same face
+    //     triangle for both the owner and neighbour cell.
+    //
+    // Each tet has its four triangles represented in the same order:
+    // 0) tri joining a tet to the tet across the face in next cell.
+    //    This is the triangle opposite Cc.
+    // 1) tri joining a tet to the tet that is in the same cell, but
+    //    belongs to the face that shares the edge of the current face
+    //    that doesn't contain basePt.  This is the triangle opposite
+    //    basePt.
+
+    // 2) tri joining a tet to the tet that is in the same cell, but
+    //    belongs to the face that shares the tet-edge (basePt - pB).
+    //    This may be on the same face, or a different one.  This is
+    //    the triangle opposite basePt.  This is the triangle opposite
+    //    pA.
+
+    // 4) tri joining a tet to the tet that is in the same cell, but
+    //    belongs to the face that shares the tet-edge (basePt - pA).
+    //    This may be on the same face, or a different one.  This is
+    //    the triangle opposite basePt.  This is the triangle opposite
+    //    pA.
+
+    // Which tri (0..3) of the tet has been crossed
+    label triI = -1;
+
+    // Determine which face was actually crossed.  lambdaMin < SMALL
+    // is considered a trigger for a tracking correction towards the
+    // current tet centre.
+    scalar lambdaMin = VGREAT;
+
+    DynamicList<label>& tris = cloud_.labels_;
+
+    // Tet indices that will be set by hitWallFaces if a wall face is
+    // to be hit, or are set when any wall tri of a tet is hit.
+    // Carries the description of the tet on which the cell face has
+    // been hit.  For the case of being set in hitWallFaces, this may
+    // be a different tet to the one that the particle occupies.
+    tetIndices faceHitTetIs;
+
+    do
     {
-            lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
-            facei_ = faces[0];
+        if (triI != -1)
+        {
+            // Change tet ownership because a tri face has been crossed
+            tetNeighbour(triI);
+        }
+
+        const Foam::face& f = pFaces[tetFaceI_];
+
+        bool own = (mesh.faceOwner()[tetFaceI_] == cellI_);
+
+        label tetBasePtI = mesh.tetBasePtIs()[tetFaceI_];
+
+        label basePtI = f[tetBasePtI];
+
+        label facePtI = (tetPtI_ + tetBasePtI) % f.size();
+        label otherFacePtI = f.fcIndex(facePtI);
+
+        label fPtAI = -1;
+        label fPtBI = -1;
+
+        if (own)
+        {
+            fPtAI = facePtI;
+            fPtBI = otherFacePtI;
         }
         else
         {
-            // If the particle has to cross more than one cell to reach the
-            // endPosition, we check which way to go.
-            // If one of the faces is a boundary face and the particle is
-            // outside, we choose the boundary face.
-            // The particle is outside if one of the lambda's is > 1 or < 0
-            forAll(faces, i)
+            fPtAI = otherFacePtI;
+            fPtBI = facePtI;
+        }
+
+        tetPointRef tet
+        (
+            pC[cellI_],
+            pPts[basePtI],
+            pPts[f[fPtAI]],
+            pPts[f[fPtBI]]
+        );
+
+        if (lambdaMin < SMALL)
         {
-                scalar lam =
-                    lambda(position_, endPosition, faces[i], stepFraction_);
+            // Apply tracking correction towards tet centre
+
+            position_ +=
+                Cloud<ParticleType>::trackingCorrectionTol
+               *(tet.centre() - position_);
+
+            cloud_.trackingRescue();
+
+            return trackFraction;
+        }
+
+        if (triI != -1 && mesh.moving())
+        {
+            // Mesh motion requires stepFraction to be correct for
+            // each tracking portion, so trackToFace must return after
+            // every lambda calculation.
+            return trackFraction;
+        }
+
+        FixedList<vector, 4> tetAreas;
+
+        tetAreas[0] = tet.Sa();
+        tetAreas[1] = tet.Sb();
+        tetAreas[2] = tet.Sc();
+        tetAreas[3] = tet.Sd();
+
+        FixedList<label, 4> tetPlaneBasePtIs;
+
+        tetPlaneBasePtIs[0] = basePtI;
+        tetPlaneBasePtIs[1] = f[fPtAI];
+        tetPlaneBasePtIs[2] = basePtI;
+        tetPlaneBasePtIs[3] = basePtI;
+
+        findTris(endPosition, tris, tet, tetAreas, tetPlaneBasePtIs);
+
+        // Reset variables for new track
+        triI = -1;
+        lambdaMin = VGREAT;
+
+        // Sets a value for lambdaMin and faceI_ if a wall face is hit
+        // by the track.
+        hitWallFaces(position_, endPosition, lambdaMin, faceHitTetIs);
+
+        // Did not hit any tet tri faces, and no wall face has been
+        // found to hit.
+        if (tris.empty() && faceI_ < 0)
+        {
+            position_ = endPosition;
+
+            return 1.0;
+        }
+        else
+        {
+            // Loop over all found tris and see if any of them find a
+            // lambda value smaller than that found for a wall face.
+            forAll(tris, i)
+            {
+                label tI = tris[i];
+
+                scalar lam = tetLambda
+                (
+                    position_,
+                    endPosition,
+                    triI,
+                    tetAreas[tI],
+                    tetPlaneBasePtIs[tI],
+                    cellI_,
+                    tetFaceI_,
+                    tetPtI_
+                );
 
                 if (lam < lambdaMin)
                 {
                     lambdaMin = lam;
-                    facei_ = faces[i];
+
+                    triI = tI;
                 }
             }
         }
 
-        bool internalFace = cloud_.internalFace(facei_);
+        if (triI == 0)
+        {
+            // This must be a cell face crossing
+            faceI_ = tetFaceI_;
 
-        // For warped faces the particle can be 'outside' the cell.
-        // This will yield a lambda larger than 1, or smaller than 0
-        // For values < 0, the particle travels away from the cell
-        // and we don't move the particle, only change cell.
-        // For values larger than 1, we move the particle to endPosition only.
-        if (lambdaMin > 0.0)
+            // Set the faceHitTetIs to those for the current tet in case a
+            // wall interaction is required with the cell face
+            faceHitTetIs = tetIndices
+            (
+                cellI_,
+                tetFaceI_,
+                tetBasePtI,
+                fPtAI,
+                fPtBI,
+                tetPtI_
+            );
+        }
+        else if (triI > 0)
+        {
+            // A tri was found to be crossed before a wall face was hit (if any)
+            faceI_ = -1;
+        }
+
+        // The particle can be 'outside' the tet.  This will yield a
+        // lambda larger than 1, or smaller than 0.  For values < 0,
+        // the particle travels away from the tet and we don't move
+        // the particle, only change tet/cell.  For values larger than
+        // 1, we move the particle to endPosition before the tet/cell
+        // change.
+        if (lambdaMin > SMALL)
         {
             if (lambdaMin <= 1.0)
             {
-                trackFraction = lambdaMin;
-                position_ += trackFraction*(endPosition - position_);
+                trackFraction += lambdaMin*(1 - trackFraction);
+
+                position_ += lambdaMin*(endPosition - position_);
             }
             else
             {
                 trackFraction = 1.0;
+
                 position_ = endPosition;
             }
         }
-        else if (static_cast<ParticleType&>(*this).softImpact())
+        else
         {
-            // Soft-sphere particles can travel outside the domain
-            // but we don't use lambda since this the particle
-            // is going away from face
-            trackFraction = 1.0;
-            position_ = endPosition;
+            // Set lambdaMin to zero to force a towards-tet-centre
+            // correction.
+            lambdaMin = 0.0;
         }
 
-        // change cell
-        if (internalFace) // Internal face
+    } while (faceI_ < 0);
+
+    if (cloud_.internalFace(faceI_))
     {
-            if (celli_ == mesh.faceOwner()[facei_])
+        // Change tet ownership because a tri face has been crossed,
+        // in general this is:
+        //     tetNeighbour(triI);
+        // but triI must be 0;
+        // No modifications are required for triI = 0, no call required to
+        //     tetNeighbour(0);
+
+        if (cellI_ == mesh.faceOwner()[faceI_])
         {
-                celli_ = mesh.faceNeighbour()[facei_];
+            cellI_ = mesh.faceNeighbour()[faceI_];
         }
-            else if (celli_ == mesh.faceNeighbour()[facei_])
+        else if (cellI_ == mesh.faceNeighbour()[faceI_])
         {
-                celli_ = mesh.faceOwner()[facei_];
+            cellI_ = mesh.faceOwner()[faceI_];
         }
         else
         {
@@ -327,24 +506,31 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     {
         ParticleType& p = static_cast<ParticleType&>(*this);
 
-            // Soft-sphere algorithm ignores the boundary
-            if (p.softImpact())
-            {
-                trackFraction = 1.0;
-                position_ = endPosition;
-            }
+        label origFaceI = faceI_;
+        label patchI = patch(faceI_);
 
-            label origFacei = facei_;
-            label patchi = patch(facei_);
+        // No action taken for tetPtI_ for tetFaceI_ here, handled by
+        // patch interaction call or later during processor transfer.
 
-            if (!p.hitPatch(mesh.boundaryMesh()[patchi], td, patchi))
+        if
+        (
+            !p.hitPatch
+            (
+                mesh.boundaryMesh()[patchI],
+                td,
+                patchI,
+                trackFraction,
+                faceHitTetIs
+            )
+        )
         {
             // Did patch interaction model switch patches?
-                if (facei_ != origFacei)
+            if (faceI_ != origFaceI)
             {
-                    patchi = patch(facei_);
+                patchI = patch(faceI_);
             }
-                const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+            const polyPatch& patch = mesh.boundaryMesh()[patchI];
 
             if (isA<wedgePolyPatch>(patch))
             {
@@ -378,7 +564,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             {
                 p.hitWallPatch
                 (
-                        static_cast<const wallPolyPatch&>(patch), td
+                    static_cast<const wallPolyPatch&>(patch), td, faceHitTetIs
                 );
             }
             else
@@ -387,22 +573,76 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             }
         }
     }
+
+    if (lambdaMin < SMALL)
+    {
+        // Apply tracking correction towards tet centre.
+        // Generate current tet to find centre to apply correction.
+
+        tetPointRef tet = currentTet();
+
+        position_ +=
+            Cloud<ParticleType>::trackingCorrectionTol
+           *(tet.centre() - position_);
+
+        if
+        (
+            cloud_.hasWallImpactDistance()
+            && !cloud_.internalFace(faceHitTetIs.face())
+            && cloud_.cellHasWallFaces()[faceHitTetIs.cell()]
+        )
+        {
+            const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+            label fI = faceHitTetIs.face();
+
+            label patchI = patches.patchID()[fI - mesh.nInternalFaces()];
+
+            if (isA<wallPolyPatch>(patches[patchI]))
+            {
+                // In the case of collision with a wall where there is
+                // a non-zero wallImpactDistance, it is possible for
+                // there to be a tracking correction required to bring
+                // the particle into the domain, but the position of
+                // the particle is further from the wall than the tet
+                // centre, in which case the normal correction can be
+                // counter-productive, i.e. pushes the particle
+                // further out of the domain.  In this case it is the
+                // position that hit the wall that is in need of a
+                // rescue correction.
+
+                triPointRef wallTri = faceHitTetIs.faceTri(mesh);
+
+                tetPointRef wallTet = faceHitTetIs.tet(mesh);
+
+                vector nHat = wallTri.normal();
+                nHat /= mag(nHat);
+
+                const ParticleType& p = static_cast<const ParticleType&>(*this);
+
+                scalar r = p.wallImpactDistance(nHat);
+
+                // Removing (approximately) the wallTri normal
+                // component of the existing correction, to avoid the
+                // situation where the existing correction in the wall
+                // normal direction is larger towards the wall than
+                // the new correction is away from it.
+                position_ +=
+                    Cloud<ParticleType>::trackingCorrectionTol
+                   *(
+                        (wallTet.centre() - (position_ + r*nHat))
+                      - (nHat & (tet.centre() - position_))*nHat
+                    );
+            }
         }
 
-    // If the trackFraction = 0 something went wrong.
-    // Either the particle is flipping back and forth across a face perhaps
-    // due to velocity interpolation errors or it is in a "hole" in the mesh
-    // caused by face warpage.
-    // In both cases resolve the positional ambiguity by moving the particle
-    // slightly towards the cell-centre.
-    if (trackFraction < SMALL)
-    {
-        position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
+        cloud_.trackingRescue();
     }
 
     return trackFraction;
 }
 
+
 template<class ParticleType>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
@@ -413,6 +653,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackToFace(endPosition, dummyTd);
 }
 
+
 template<class ParticleType>
 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 {
@@ -436,7 +677,9 @@ bool Foam::Particle<ParticleType>::hitPatch
 (
     const polyPatch&,
     TrackData&,
-    const label
+    const label,
+    const scalar,
+    const tetIndices&
 )
 {
     return false;
@@ -451,7 +694,17 @@ void Foam::Particle<ParticleType>::hitWedgePatch
     TrackData&
 )
 {
-    vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
+    FatalErrorIn
+    (
+        "void Foam::Particle<ParticleType>::hitWedgePatch"
+        "("
+            "const wedgePolyPatch& wpp, "
+            "TrackData&"
+        ")"
+    )   << "Hitting a wedge patch should not be possible."
+        << abort(FatalError);
+
+    vector nf = normal();
     nf /= mag(nf);
 
     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
@@ -466,7 +719,7 @@ void Foam::Particle<ParticleType>::hitSymmetryPatch
     TrackData&
 )
 {
-    vector nf = spp.faceAreas()[spp.whichFace(facei_)];
+    vector nf = normal();
     nf /= mag(nf);
 
     static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
@@ -481,11 +734,16 @@ void Foam::Particle<ParticleType>::hitCyclicPatch
     TrackData&
 )
 {
-//    label patchFacei_ = cpp.whichFace(facei_);
+    // label patchFaceI_ = cpp.whichFace(faceI_);
 
-    facei_ = cpp.transformGlobalFace(facei_);
+    faceI_ = cpp.transformGlobalFace(faceI_);
 
-    celli_ = cloud_.polyMesh_.faceOwner()[facei_];
+    cellI_ = cloud_.polyMesh_.faceOwner()[faceI_];
+
+    tetFaceI_ = faceI_;
+
+    // See note in correctAfterParallelTransfer for tetPtI_ addressing.
+    tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
 
     // Now the particle is on the receiving side
 
@@ -522,7 +780,8 @@ template<class TrackData>
 void Foam::Particle<ParticleType>::hitWallPatch
 (
     const wallPolyPatch& spp,
-    TrackData&
+    TrackData&,
+    const tetIndices&
 )
 {}
 

src/lagrangian/basic/Particle/Particle.H

@@ -38,6 +38,10 @@ Description
 #include "faceList.H"
 #include "typeInfo.H"
 #include "OFstream.H"
+#include "tetPointRef.H"
+#include "FixedList.H"
+#include "polyMeshTetDecomposition.H"
+#include "wallPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -125,14 +129,23 @@ protected:
         vector position_;
 
         //- Index of the cell it is in
-        label celli_;
+        label cellI_;
 
         //- Face index if the particle is on a face otherwise -1
-        label facei_;
+        label faceI_;
 
         //- Fraction of time-step completed
         scalar stepFraction_;
 
+        //- Index of the face that owns the decomposed tet that the
+        //  particle is in
+        label tetFaceI_;
+
+        //- Index of the point on the face that defines the decomposed
+        //  tet that the particle is in.  Relative to the face base
+        //  point.
+        label tetPtI_;
+
         //- Originating processor id
         label origProc_;
 
@@ -142,54 +155,82 @@ protected:
 
     // Private Member Functions
 
-        //- Return the 'lambda' value for the position, p, on the face,
-        // where, p = from + lamda*(to - from)
-        // for non-static meshes
-        inline scalar lambda
+        //- Find the tet tri faces between position and tet centre
+        inline void findTris
+        (
+            const vector& position,
+            DynamicList<label>& faceList,
+            const tetPointRef& tet,
+            const FixedList<vector, 4>& tetAreas,
+            const FixedList<label, 4>& tetPlaneBasePtIs
+        ) const;
+
+        //- Find the lambda value for the line to-from across the
+        //  given tri face, where p = from + lambda*(to - from)
+        inline scalar tetLambda
         (
             const vector& from,
             const vector& to,
-            const label facei,
-            const scalar stepFraction
+            const label triI,
+            const vector& tetArea,
+            const label tetPlaneBasePtI,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI
         ) const;
 
-        //- Return the 'lambda' value for the position, p, on the face,
-        // where, p = from + lamda*(to - from)
-        // for static meshes
-        inline scalar lambda
+         //- Find the lambda value for a moving tri face
+         inline scalar movingTetLambda
          (
              const vector& from,
              const vector& to,
-            const label facei
+             const label triI,
+             const vector& tetArea,
+             const label tetPlaneBasePtI,
+             const label cellI,
+             const label tetFaceI,
+             const label tetPtI
          ) const;
 
-        //- Find the faces between position and cell centre
-        void findFaces
-        (
-            const vector& position,
-            DynamicList<label>& faceList
-        ) const;
+        //- Modify the tet owner data by crossing triI
+        inline void tetNeighbour(label triI);
 
-        //- Find the faces between position and cell centre
-        void findFaces
+        //- Cross the from the given face across the given edge of the
+        //  given cell to find the resulting face and tetPtI
+        inline void crossEdgeConnectedFace
         (
-            const vector& position,
-            const label celli,
-            const scalar stepFraction,
-            DynamicList<label>& faceList
-        ) const;
+            const label& cellI,
+            label& tetFaceI,
+            label& tetPtI,
+            const edge& e
+        );
+
+        //- Hit wall faces in the current cell if the
+        //- wallImpactDistance is non-zero.  They may not be in
+        //- different tets to the current.
+        inline void hitWallFaces
+        (
+            const vector& from,
+            const vector& to,
+            scalar& lambdaMin,
+            tetIndices& closestTetIs
+        );
 
 
     // Patch interactions
 
-        //- Overridable function to handle the particle hitting a patch
-        //  Executed before other patch-hitting functions
+        //- Overridable function to handle the particle hitting a
+        //  patch.  Executed before other patch-hitting functions.
+        //  trackFraction is passed in to allow mesh motion to
+        //  interpolate in time to the correct face state.
         template<class TrackData>
         bool hitPatch
         (
             const polyPatch&,
             TrackData& td,
-            const label patchI
+            const label patchI,
+            const scalar trackFraction,
+            const tetIndices& tetIs
         );
 
         //- Overridable function to handle the particle hitting a wedgePatch
@@ -231,7 +272,8 @@ protected:
         void hitWallPatch
         (
             const wallPolyPatch&,
-            TrackData& td
+            TrackData& td,
+            const tetIndices& tetIs
         );
 
         //- Overridable function to handle the particle hitting a
@@ -264,12 +306,12 @@ protected:
         //- Convert global addressing to the processor patch
         //  local equivalents
         template<class TrackData>
-        void prepareForParallelTransfer(const label patchi, TrackData& td);
+        void prepareForParallelTransfer(const label patchI, TrackData& td);
 
         //- Convert processor patch addressing to the global equivalents
-        //  and set the celli to the face-neighbour
+        //  and set the cellI to the face-neighbour
         template<class TrackData>
-        void correctAfterParallelTransfer(const label patchi, TrackData& td);
+        void correctAfterParallelTransfer(const label patchI, TrackData& td);
 
 
 public:
@@ -293,7 +335,18 @@ public:
         (
             const Cloud<ParticleType>&,
             const vector& position,
-            const label celli
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI
+        );
+
+        //- Construct from components, tetFaceI_ and tetPtI_ are not
+        //  supplied so they will be deduced by a search
+        Particle
+        (
+            const Cloud<ParticleType>&,
+            const vector& position,
+            const label cellI
         );
 
         //- Construct from Istream
@@ -355,17 +408,6 @@ public:
 
         // Access
 
-            //- Return true if particle is in cell
-            inline bool inCell() const;
-
-            //- Return true if position is in cell i
-            inline bool inCell
-            (
-                const vector& position,
-                const label celli,
-                const scalar stepFraction
-            ) const;
-
             //- Return current particle position
             inline const vector& position() const;
 
@@ -378,6 +420,35 @@ public:
             //- Return current cell particle is in
             inline label cell() const;
 
+            //- Return current tet face particle is in
+            inline label& tetFace();
+
+            //- Return current tet face particle is in
+            inline label tetFace() const;
+
+            //- Return current tet face particle is in
+            inline label& tetPt();
+
+            //- Return current tet face particle is in
+            inline label tetPt() const;
+
+            //- Return the indices of the current tet that the
+            //  particle occupies.
+            inline tetIndices currentTetIndices() const;
+
+            //- Return the geometry of the current tet that the
+            //  particle occupies.
+            inline tetPointRef currentTet() const;
+
+            //- Return the normal of the tri on tetFaceI_ for the
+            //  current tet.
+            inline vector normal() const;
+
+            //- Return the normal of the tri on tetFaceI_ for the
+            //  current tet at the start of the timestep, i.e. based
+            //  on oldPoints
+            inline vector oldNormal() const;
+
             //- Return current face particle is on otherwise -1
             inline label& face();
 
@@ -396,17 +467,21 @@ public:
 
         // Check
 
+            //- Check the stored cell value (setting if necessary) and
+            //  initialise the tetFace and tetPt values
+            inline void initCellFacePt();
+
             //- Is the particle on the boundary/(or outside the domain)?
             inline bool onBoundary() const;
 
             //- Which patch is particle on
-            inline label patch(const label facei) const;
+            inline label patch(const label faceI) const;
 
             //- Which face of this patch is this particle on
             inline label patchFace
             (
-                const label patchi,
-                const label facei
+                const label patchI,
+                const label faceI
             ) const;
 
             //- The nearest distance to a wall that

src/lagrangian/basic/Particle/ParticleIO.C

@@ -31,7 +31,8 @@ License
 
 template<class ParticleType>
 Foam::string Foam::Particle<ParticleType>::propHeader =
-    "(Px Py Pz) cellI origProc origId";
+    "(Px Py Pz) cellI tetFaceI tetPtI origProc origId";
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -44,20 +45,24 @@ Foam::Particle<ParticleType>::Particle
 )
 :
     cloud_(cloud),
-    facei_(-1),
+    position_(),
+    cellI_(-1),
+    faceI_(-1),
     stepFraction_(0.0),
+    tetFaceI_(-1),
+    tetPtI_(-1),
     origProc_(Pstream::myProcNo()),
     origId_(-1)
 {
-
     // readFields : read additional data. Should be consistent with writeFields.
 
     if (is.format() == IOstream::ASCII)
     {
-        is >> position_ >> celli_;
+        is  >> position_ >> cellI_;
+
         if (readFields)
         {
-            is >> origProc_ >> origId_;
+            is  >> tetFaceI_ >> tetPtI_ >> origProc_ >> origId_;
         }
     }
     else
@@ -69,9 +74,11 @@ Foam::Particle<ParticleType>::Particle
             (
                 reinterpret_cast<char*>(&position_),
                 sizeof(position_)
-              + sizeof(celli_)
-              + sizeof(facei_)
+              + sizeof(cellI_)
+              + sizeof(faceI_)
               + sizeof(stepFraction_)
+              + sizeof(tetFaceI_)
+              + sizeof(tetPtI_)
               + sizeof(origProc_)
               + sizeof(origId_)
             );
@@ -82,18 +89,13 @@ Foam::Particle<ParticleType>::Particle
             (
                 reinterpret_cast<char*>(&position_),
                 sizeof(position_)
-              + sizeof(celli_)
-              + sizeof(facei_)
+              + sizeof(cellI_)
+              + sizeof(faceI_)
               + sizeof(stepFraction_)
             );
         }
     }
 
-    if (celli_ == -1)
-    {
-        celli_ = cloud_.pMesh().findCell(position_);
-    }
-
     // Check state of Istream
     is.check("Particle<ParticleType>::Particle(Istream&)");
 }
@@ -177,28 +179,32 @@ void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
         {
             // Write the additional entries
             os  << position_
-               << token::SPACE << celli_
+                << token::SPACE << cellI_
+                << token::SPACE << tetFaceI_
+                << token::SPACE << tetPtI_
                 << token::SPACE << origProc_
                 << token::SPACE << origId_;
         }
         else
         {
             os  << position_
-               << token::SPACE << celli_;
+                << token::SPACE << cellI_;
         }
     }
     else
     {
-        // In binary write both celli_ and facei_, needed for parallel transfer
+        // In binary write both cellI_ and faceI_, needed for parallel transfer
         if (writeFields)
         {
             os.write
             (
                 reinterpret_cast<const char*>(&position_),
                 sizeof(position_)
-              + sizeof(celli_)
-              + sizeof(facei_)
+              + sizeof(cellI_)
+              + sizeof(faceI_)
               + sizeof(stepFraction_)
+              + sizeof(tetFaceI_)
+              + sizeof(tetPtI_)
               + sizeof(origProc_)
               + sizeof(origId_)
             );
@@ -209,8 +215,8 @@ void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
             (
                 reinterpret_cast<const char*>(&position_),
                 sizeof(position_)
-              + sizeof(celli_)
-              + sizeof(facei_)
+              + sizeof(cellI_)
+              + sizeof(faceI_)
               + sizeof(stepFraction_)
             );
         }

src/lagrangian/basic/indexedParticle/indexedParticle.H

@@ -67,11 +67,26 @@ public:
         (
             const Cloud<indexedParticle>& c,
             const vector& position,
-            const label celli,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
             const label index = 0
         )
         :
-            Particle<indexedParticle>(c, position, celli),
+            Particle<indexedParticle>(c, position, cellI, tetFaceI, tetPtI),
+            index_(index)
+        {}
+
+        //- Construct from components, with searching for tetFace and tetPt
+        indexedParticle
+        (
+            const Cloud<indexedParticle>& c,
+            const vector& position,
+            const label cellI,
+            const label index = 0
+        )
+        :
+            Particle<indexedParticle>(c, position, cellI),
             index_(index)
         {}
 

src/lagrangian/basic/passiveParticle/passiveParticle.H

@@ -60,10 +60,23 @@ public:
         (
             const Cloud<passiveParticle>& c,
             const vector& position,
-            const label celli
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI
         )
         :
-            Particle<passiveParticle>(c, position, celli)
+            Particle<passiveParticle>(c, position, cellI, tetFaceI, tetPtI)
+        {}
+
+        //- Construct from components, with searching for tetFace and tetPt
+        passiveParticle
+        (
+            const Cloud<passiveParticle>& c,
+            const vector& position,
+            const label cellI
+        )
+        :
+            Particle<passiveParticle>(c, position, cellI)
         {}
 
         //- Construct from Istream

src/lagrangian/coalCombustion/coalParcel/coalParcel.C

@@ -31,10 +31,19 @@ Foam::coalParcel::coalParcel
 (
     ReactingMultiphaseCloud<coalParcel>& owner,
     const vector& position,
-    const label cellI
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI
 )
 :
-    ReactingMultiphaseParcel<coalParcel>(owner, position, cellI)
+    ReactingMultiphaseParcel<coalParcel>
+    (
+        owner,
+        position,
+        cellI,
+        tetFaceI,
+        tetPtI
+    )
 {}
 
 
@@ -43,6 +52,8 @@ Foam::coalParcel::coalParcel
     ReactingMultiphaseCloud<coalParcel>& owner,
     const vector& position,
     const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
     const label typeId,
     const scalar nParticle0,
     const scalar d0,
@@ -62,6 +73,8 @@ Foam::coalParcel::coalParcel
         owner,
         position,
         cellI,
+        tetFaceI,
+        tetPtI,
         typeId,
         nParticle0,
         d0,

src/lagrangian/coalCombustion/coalParcel/coalParcel.H

@@ -63,7 +63,9 @@ public:
         (
              ReactingMultiphaseCloud<coalParcel>& owner,
              const vector& position,
-             const label cellI
+             const label cellI,
+             const label tetFaceI,
+             const label tetPtI
         );
 
         //- Construct from components
@@ -72,6 +74,8 @@ public:
              ReactingMultiphaseCloud<coalParcel>& owner,
              const vector& position,
              const label cellI,
+             const label tetFaceI,
+             const label tetPtI,
              const label typeId,
              const scalar nParticle0,
              const scalar d0,

src/lagrangian/dieselSpray/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
@@ -16,6 +17,7 @@ EXE_INC = \
 
 LIB_LIBS = \
     -llagrangian \
+    -lmeshTools \
     -lfiniteVolume \
     -lcompressibleRASModels \
     -lcompressibleLESModels \

src/lagrangian/dieselSpray/parcel/boundaryTreatment.H

@@ -15,12 +15,11 @@ if (isA<wallPolyPatch>(pbMesh[patch(face())]))
 else if (isA<wedgePolyPatch>(pbMesh[patch(face())]))
 {
     // check if parcel is trying to move out of the domain
-    label patchi = patch(face());
-    label patchFacei = patchFace(patchi, face());
-    const polyPatch& patch = mesh.boundaryMesh()[patchi];
-    vector nf = patch.faceAreas()[patchFacei];
+
+    vector nf = normal();
 
     scalar Un = U() & nf;
+
     if (Un > 0)
     {
         scalar Un2 = U() & n();
@@ -30,12 +29,11 @@ else if (isA<wedgePolyPatch>(pbMesh[patch(face())]))
 else if (isA<symmetryPolyPatch>(pbMesh[patch(face())]))
 {
     // check if parcel is trying to move out of the domain
-    label patchi = patch(face());
-    label patchFacei = patchFace(patchi, face());
-    const polyPatch& patch = mesh.boundaryMesh()[patchi];
-    vector nf = patch.faceAreas()[patchFacei];
+
+    vector nf = normal();
 
     scalar Un = U() & nf;
+
     if (Un > 0)
     {
         if (sDB.twoD())

src/lagrangian/dieselSpray/parcel/parcel.C

@@ -50,6 +50,8 @@ Foam::parcel::parcel
     const Cloud<parcel>& cloud,
     const vector& position,
     const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
     const vector& n,
     const scalar d,
     const scalar T,
@@ -67,7 +69,7 @@ Foam::parcel::parcel
     const List<word>& liquidNames
 )
 :
-    Particle<parcel>(cloud, position, cellI),
+    Particle<parcel>(cloud, position, cellI, tetFaceI, tetPtI),
     liquidComponents_
     (
         liquidNames
@@ -103,12 +105,13 @@ bool Foam::parcel::move(spray& sDB)
     label Nf = fuels.components().size();
     label Ns = sDB.composition().Y().size();
 
+    tetIndices tetIs = this->currentTetIndices();
+
     // Calculate the interpolated gas properties at the position of the parcel
-    vector Up = sDB.UInterpolator().interpolate(position(), cell())
-        + Uturb();
-    scalar rhog = sDB.rhoInterpolator().interpolate(position(), cell());
-    scalar pg = sDB.pInterpolator().interpolate(position(), cell());
-    scalar Tg = sDB.TInterpolator().interpolate(position(), cell());
+    vector Up = sDB.UInterpolator().interpolate(position(), tetIs) + Uturb();
+    scalar rhog = sDB.rhoInterpolator().interpolate(position(), tetIs);
+    scalar pg = sDB.pInterpolator().interpolate(position(), tetIs);
+    scalar Tg = sDB.TInterpolator().interpolate(position(), tetIs);
 
     scalarField Yfg(Nf, 0.0);
 
@@ -119,8 +122,8 @@ bool Foam::parcel::move(spray& sDB)
         if (sDB.isLiquidFuel()[i])
         {
             label j = sDB.gasToLiquidIndex()[i];
-            scalar Yicelli = Yi[cell()];
-            Yfg[j] = Yicelli;
+            scalar YicellI = Yi[cell()];
+            Yfg[j] = YicellI;
         }
         cpMixture += Yi[cell()]*sDB.gasProperties()[i].Cp(Tg);
     }
@@ -199,8 +202,8 @@ bool Foam::parcel::move(spray& sDB)
 
         // remember which cell the parcel is in
         // since this will change if a face is hit
-        label celli = cell();
-        scalar p = sDB.p()[celli];
+        label cellI = cell();
+        scalar p = sDB.p()[cellI];
 
         // track parcel to face, or end of trajectory
         if (keepParcel)
@@ -259,7 +262,7 @@ bool Foam::parcel::move(spray& sDB)
         (
             dt,
             sDB,
-            celli,
+            cellI,
             face()
         );
 
@@ -283,11 +286,11 @@ bool Foam::parcel::move(spray& sDB)
         // Update the Spray Source Terms
         forAll(nMass, i)
         {
-            sDB.srhos()[i][celli] += oMass[i] - nMass[i];
+            sDB.srhos()[i][cellI] += oMass[i] - nMass[i];
         }
-        sDB.sms()[celli] += oMom - nMom;
+        sDB.sms()[cellI] += oMom - nMom;
 
-        sDB.shs()[celli] += oTotMass*(oH + oPE) - m()*(nH + nPE);
+        sDB.shs()[cellI] += oTotMass*(oH + oPE) - m()*(nH + nPE);
 
         // Remove evaporated mass from stripped mass
         ms() -= ms()*(oTotMass-m())/oTotMass;
@@ -300,11 +303,11 @@ bool Foam::parcel::move(spray& sDB)
             // ... and add the removed 'stuff' to the gas
             forAll(nMass, i)
             {
-                sDB.srhos()[i][celli] += nMass[i];
+                sDB.srhos()[i][cellI] += nMass[i];
             }
 
-            sDB.sms()[celli] += nMom;
-            sDB.shs()[celli] += m()*(nH + nPE);
+            sDB.sms()[cellI] += nMom;
+            sDB.shs()[cellI] += m()*(nH + nPE);
         }
 
         if (onBoundary() && keepParcel)
@@ -327,8 +330,8 @@ void Foam::parcel::updateParcelProperties
 (
     const scalar dt,
     spray& sDB,
-    const label celli,
-    const label facei
+    const label cellI,
+    const label faceI
 )
 {
     const liquidMixture& fuels = sDB.fuels();
@@ -340,34 +343,34 @@ void Foam::parcel::updateParcelProperties
     scalar W = 0.0;
     for (label i=0; i<Ns; i++)
     {
-        W += sDB.composition().Y()[i][celli]/sDB.gasProperties()[i].W();
+        W += sDB.composition().Y()[i][cellI]/sDB.gasProperties()[i].W();
 
     }
     W = 1.0/W;
 
     // Calculate the interpolated gas properties at the position of the parcel
-    vector Up = sDB.UInterpolator().interpolate(position(), celli, facei)
+    vector Up = sDB.UInterpolator().interpolate(position(), cellI, faceI)
         + Uturb();
-    scalar rhog = sDB.rhoInterpolator().interpolate(position(), celli, facei);
-    scalar pg = sDB.pInterpolator().interpolate(position(), celli, facei);
-    scalar Tg0 = sDB.TInterpolator().interpolate(position(), celli, facei);
+    scalar rhog = sDB.rhoInterpolator().interpolate(position(), cellI, faceI);
+    scalar pg = sDB.pInterpolator().interpolate(position(), cellI, faceI);
+    scalar Tg0 = sDB.TInterpolator().interpolate(position(), cellI, faceI);
 
     // correct the gaseous temperature for evaporated fuel
     scalar cpMix = 0.0;
     for (label i=0; i<Ns; i++)
     {
-        cpMix += sDB.composition().Y()[i][celli]
+        cpMix += sDB.composition().Y()[i][cellI]
                 *sDB.gasProperties()[i].Cp(T());
     }
-    scalar cellV            = sDB.mesh().V()[celli];
-    scalar rho              = sDB.rho()[celli];
+    scalar cellV            = sDB.mesh().V()[cellI];
+    scalar rho              = sDB.rho()[cellI];
     scalar cellMass         = rho*cellV;
-    scalar dh               = sDB.shs()[celli];
+    scalar dh               = sDB.shs()[cellI];
     scalarField addedMass(Nf, 0.0);
 
     forAll(addedMass, i)
     {
-        addedMass[i] += sDB.srhos()[i][celli]*cellV;
+        addedMass[i] += sDB.srhos()[i][cellI]*cellV;
     }
 
     scalar Tg  = Tg0 + dh/(cpMix*cellMass);
@@ -378,7 +381,7 @@ void Foam::parcel::updateParcelProperties
     {
         label j = sDB.liquidToGasIndex()[i];
         const volScalarField& Yj = sDB.composition().Y()[j];
-        scalar Yfg0 = Yj[celli];
+        scalar Yfg0 = Yj[cellI];
         Yfg[i] = (Yfg0*cellMass + addedMass[i])/(addedMass[i] + cellMass);
     }
 
@@ -389,7 +392,7 @@ void Foam::parcel::updateParcelProperties
 
     setRelaxationTimes
     (
-        celli,
+        cellI,
         tauMomentum,
         tauEvaporation,
         tauHeatTransfer,
@@ -492,14 +495,14 @@ void Foam::parcel::updateParcelProperties
             {
                 label j = sDB.liquidToGasIndex()[i];
                 const volScalarField& Yj = sDB.composition().Y()[j];
-                scalar Yfg0 = Yj[celli];
+                scalar Yfg0 = Yj[cellI];
                 Yfg[i] = (Yfg0*cellMass + addedMass[i] + dm)
                         /(addedMass[i] + cellMass + dm);
             }
 
             setRelaxationTimes
             (
-                celli,
+                cellI,
                 tauMomentum,
                 tauEvaporation,
                 tauHeatTransfer,
@@ -530,7 +533,7 @@ void Foam::parcel::updateParcelProperties
             }
             else
             {
-                scalar Y = sDB.composition().Y()[i][celli];
+                scalar Y = sDB.composition().Y()[i][cellI];
                 cpMix += Y*sDB.gasProperties()[i].Cp(Taverage);
             }
         }

src/lagrangian/dieselSpray/parcel/parcel.H

@@ -113,7 +113,7 @@ class parcel
         //- Set the relaxation times
         void setRelaxationTimes
         (
-            label celli,
+            label cellI,
             scalar& tauMomentum,
             scalarField& tauEvaporation,
             scalar& tauHeatTransfer,
@@ -133,8 +133,8 @@ class parcel
         (
             const scalar dt,
             spray& sprayData,
-            const label celli,
-            const label facei
+            const label cellI,
+            const label faceI
         );
 
 
@@ -150,7 +150,9 @@ public:
         (
             const Cloud<parcel>& cloud,
             const vector& position,
-            const label celli,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
             const vector& n,
             const scalar d,
             const scalar T,

src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C

@@ -36,7 +36,7 @@ License
 
 void Foam::parcel::setRelaxationTimes
 (
-    label celli,
+    label cellI,
     scalar& tauMomentum,
     scalarField& tauEvaporation,
     scalar& tauHeatTransfer,
@@ -70,7 +70,7 @@ void Foam::parcel::setRelaxationTimes
 
     for (label i=0; i<Ns; i++)
     {
-        scalar Y = sDB.composition().Y()[i][celli];
+        scalar Y = sDB.composition().Y()[i][cellI];
         W += Y/sDB.gasProperties()[i].W();
         // Using mass-fractions to average...
         kMixture += Y*sDB.gasProperties()[i].kappa(Tf);
@@ -87,7 +87,7 @@ void Foam::parcel::setRelaxationTimes
     for (label i=0; i<Nf; i++)
     {
         label j = sDB.liquidToGasIndex()[i];
-        scalar Y = sDB.composition().Y()[j][celli];
+        scalar Y = sDB.composition().Y()[j][cellI];
         scalar Wi = sDB.gasProperties()[j].W();
         Yf[i] = Y;
         Xf[i] = Y*W/Wi;
@@ -264,10 +264,10 @@ void Foam::parcel::setRelaxationTimes
                 forAll(sDB.gasProperties(), k)
                 {
                     vapourSurfaceEnthalpy +=
-                        sDB.composition().Y()[k][celli]
+                        sDB.composition().Y()[k][cellI]
                        *sDB.gasProperties()[k].H(tBoilingSurface);
                     vapourFarEnthalpy +=
-                        sDB.composition().Y()[k][celli]
+                        sDB.composition().Y()[k][cellI]
                        *sDB.gasProperties()[k].H(temperature);
                 }
 

src/lagrangian/dieselSpray/spray/findInjectorCell.H

@@ -22,7 +22,14 @@ reduce(foundCell, orOp<bool>());
 if (!foundCell)
 {
     injectionPosition = it->position(n);
-    injectorCell = mesh_.findCell(injectionPosition);
+
+    findCellFacePt
+    (
+        injectionPosition,
+        injectorCell,
+        injectorTetFaceI,
+        injectorTetPtI
+    );
 
     if (injectorCell >= 0)
     {

src/lagrangian/dieselSpray/spray/spray.C

@@ -241,16 +241,16 @@ Foam::spray::spray
     label n=0;
 
     // check for the type of boundary condition
-    forAll(bMesh, patchi)
+    forAll(bMesh, patchI)
     {
-        if (isA<symmetryPolyPatch>(bMesh[patchi]))
+        if (isA<symmetryPolyPatch>(bMesh[patchI]))
         {
             symPlaneExist = true;
         }
-        else if (isA<wedgePolyPatch>(bMesh[patchi]))
+        else if (isA<wedgePolyPatch>(bMesh[patchI]))
         {
             wedgeExist = true;
-            patches[n++] = patchi;
+            patches[n++] = patchI;
         }
     }
 

src/lagrangian/dieselSpray/spray/sprayFunctions.C

@@ -149,9 +149,9 @@ Foam::scalar Foam::spray::liquidTotalEnthalpy() const
 
     forAllConstIter(spray, *this, iter)
     {
-        label celli = iter().cell();
+        label cellI = iter().cell();
         scalar T = iter().T();
-        scalar pc = p()[celli];
+        scalar pc = p()[cellI];
         scalar rho = fuels().rho(pc, T, iter().X());
         scalar hlat = fuels().hl(pc, T, iter().X());
         scalar hg = 0.0;
@@ -351,8 +351,8 @@ Foam::scalar Foam::spray::smd() const
 
     forAllConstIter(spray, *this, iter)
     {
-        label celli = iter().cell();
-        scalar Pc = p()[celli];
+        label cellI = iter().cell();
+        scalar Pc = p()[cellI];
         scalar T = iter().T();
         scalar rho = fuels_->rho(Pc, T, iter().X());
 

src/lagrangian/dieselSpray/spray/sprayInject.C

@@ -107,7 +107,17 @@ void Foam::spray::inject()
                     scalar deviation = breakup().y0();
                     scalar ddev = breakup().yDot0();
 
-                    label injectorCell = mesh_.findCell(injectionPosition);
+                    label injectorCell = -1;
+                    label injectorTetFaceI = -1;
+                    label injectorTetPtI = -1;
+
+                    findCellFacePt
+                    (
+                        injectionPosition,
+                        injectorCell,
+                        injectorTetFaceI,
+                        injectorTetPtI
+                    );
 
 #                   include "findInjectorCell.H"
 
@@ -122,6 +132,8 @@ void Foam::spray::inject()
                             *this,
                             injectionPosition,
                             injectorCell,
+                            injectorTetFaceI,
+                            injectorTetPtI,
                             normal,
                             diameter,
                             it->T(toi),

src/lagrangian/dieselSpray/spray/sprayOps.C

@@ -88,7 +88,7 @@ void Foam::spray::breakupLoop()
         vector velocity = UInterpolator().interpolate
         (
             elmnt().position(),
-            elmnt().cell()
+            elmnt().currentTetIndices()
         );
 
         // liquidCore < 0.5 indicates discrete drops
@@ -122,7 +122,7 @@ void Foam::spray::atomizationLoop()
         vector velocity = UInterpolator().interpolate
         (
             elmnt().position(),
-            elmnt().cell()
+            elmnt().currentTetIndices()
         );
 
         // liquidCore > 0.5 indicates a liquid core

src/lagrangian/dieselSpray/spraySubModels/breakupModel/SHF/SHF.C

@@ -102,14 +102,14 @@ void Foam::SHF::breakupParcel
     const liquidMixture& fuels
 ) const
 {
-    label celli = p.cell();
+    label cellI = p.cell();
     scalar T = p.T();
-    scalar pc = spray_.p()[celli];
+    scalar pc = spray_.p()[cellI];
 
     scalar sigma = fuels.sigma(pc, T, p.X());
     scalar rhoLiquid = fuels.rho(pc, T, p.X());
     scalar muLiquid = fuels.mu(pc, T, p.X());
-    scalar rhoGas = spray_.rho()[celli];
+    scalar rhoGas = spray_.rho()[cellI];
 
     scalar weGas = p.We(vel, rhoGas, sigma);
     scalar weLiquid = p.We(vel, rhoLiquid, sigma);
@@ -237,6 +237,8 @@ void Foam::SHF::breakupParcel
                     spray_,
                     p.position(),
                     p.cell(),
+                    p.tetFace(),
+                    p.tetPt(),
                     p.n(),
                     d,
                     p.T(),

src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C

@@ -78,15 +78,15 @@ void Foam::reitzKHRT::breakupParcel
     const liquidMixture& fuels
 ) const
 {
-    label celli = p.cell();
+    label cellI = p.cell();
     scalar T = p.T();
     scalar r = 0.5*p.d();
-    scalar pc = spray_.p()[celli];
+    scalar pc = spray_.p()[cellI];
 
     scalar sigma = fuels.sigma(pc, T, p.X());
     scalar rhoLiquid = fuels.rho(pc, T, p.X());
     scalar muLiquid = fuels.mu(pc, T, p.X());
-    scalar rhoGas = spray_.rho()[celli];
+    scalar rhoGas = spray_.rho()[cellI];
     scalar Np = p.N(rhoLiquid);
     scalar semiMass = Np*pow3(p.d());
 
@@ -198,6 +198,8 @@ void Foam::reitzKHRT::breakupParcel
                         spray_,
                         p.position(),
                         p.cell(),
+                        p.tetFace(),
+                        p.tetPt(),
                         p.n(),
                         dc,
                         p.T(),

src/lagrangian/dieselSpray/spraySubModels/wallModel/reflectParcel/reflectParcel.C

@@ -70,18 +70,18 @@ Foam::reflectParcel::~reflectParcel()
 bool Foam::reflectParcel::wallTreatment
 (
     parcel& p,
-    const label globalFacei
+    const label globalFaceI
 ) const
 {
-    label patchi = p.patch(globalFacei);
-    label facei = p.patchFace(patchi, globalFacei);
+    label patchI = p.patch(globalFaceI);
+    label faceI = p.patchFace(patchI, globalFaceI);
 
     const polyMesh& mesh = spray_.mesh();
 
-    if (isA<wallPolyPatch>(mesh_.boundaryMesh()[patchi]))
+    if (isA<wallPolyPatch>(mesh_.boundaryMesh()[patchI]))
     {
         // wallNormal defined to point outwards of domain
-        vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
+        vector Sf = mesh_.Sf().boundaryField()[patchI][faceI];
         Sf /= mag(Sf);
 
         if (!mesh.moving())
@@ -97,17 +97,17 @@ bool Foam::reflectParcel::wallTreatment
         else
         {
             // moving mesh
-            vector Ub1 = U_.boundaryField()[patchi][facei];
-            vector Ub0 = U_.oldTime().boundaryField()[patchi][facei];
+            vector Ub1 = U_.boundaryField()[patchI][faceI];
+            vector Ub0 = U_.oldTime().boundaryField()[patchI][faceI];
 
             scalar dt = spray_.runTime().deltaTValue();
             const vectorField& oldPoints = mesh.oldPoints();
 
-            const vector& Cf1 = mesh.faceCentres()[globalFacei];
+            const vector& Cf1 = mesh.faceCentres()[globalFaceI];
 
-            vector Cf0 = mesh.faces()[globalFacei].centre(oldPoints);
+            vector Cf0 = mesh.faces()[globalFaceI].centre(oldPoints);
             vector Cf = Cf0 + p.stepFraction()*(Cf1 - Cf0);
-            vector Sf0 = mesh.faces()[globalFacei].normal(oldPoints);
+            vector Sf0 = mesh.faces()[globalFaceI].normal(oldPoints);
 
             // for layer addition Sf0 = vector::zero and we use Sf
             if (mag(Sf0) > SMALL)
@@ -157,10 +157,10 @@ bool Foam::reflectParcel::wallTreatment
                 {
                     Info<< "reflectParcel:: v = " << v
                         << ", Ub = " << Ub
-                        << ", facei = " << facei
-                        << ", patchi = " << patchi
-                        << ", globalFacei = " << globalFacei
-                        << ", name = " << mesh_.boundaryMesh()[patchi].name()
+                        << ", faceI = " << faceI
+                        << ", patchI = " << patchI
+                        << ", globalFaceI = " << globalFaceI
+                        << ", name = " << mesh_.boundaryMesh()[patchI].name()
                         << endl;
                 }
             */
@@ -175,7 +175,7 @@ bool Foam::reflectParcel::wallTreatment
     else
     {
         FatalErrorIn("bool reflectParcel::wallTreatment(parcel& parcel) const")
-            << " parcel has hit a boundary " << mesh_.boundary()[patchi].type()
+            << " parcel has hit a boundary " << mesh_.boundary()[patchI].type()
             << " which not yet has been implemented." << nl
             << abort(FatalError);
     }

src/lagrangian/dieselSpray/spraySubModels/wallModel/reflectParcel/reflectParcel.H

@@ -87,7 +87,7 @@ public:
 
         //- Return true if parcel is to be kept, and false if it is to be
         // removed
-        bool wallTreatment(parcel& parcel, const label facei) const;
+        bool wallTreatment(parcel& parcel, const label faceI) const;
 };
 
 

src/lagrangian/dieselSpray/spraySubModels/wallModel/removeParcel/removeParcel.C

@@ -65,7 +65,7 @@ Foam::removeParcel::~removeParcel()
 bool Foam::removeParcel::wallTreatment
 (
     parcel&,
-    const label facei
+    const label faceI
 ) const
 {
     return false;

src/lagrangian/dieselSpray/spraySubModels/wallModel/removeParcel/removeParcel.H

@@ -73,7 +73,7 @@ public:
 
         //- Return true if parcel is to be kept, and false if it is to be
         // removed
-        bool wallTreatment(parcel& parcel, const label facei) const;
+        bool wallTreatment(parcel& parcel, const label faceI) const;
 };
 
 

src/lagrangian/dieselSpray/spraySubModels/wallModel/wallModel/wallModel.H

@@ -118,7 +118,7 @@ public:
         virtual bool wallTreatment
         (
             parcel&,
-            const label facei
+            const label faceI
         ) const = 0;
 };
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C

@@ -109,57 +109,21 @@ void Foam::DsmcCloud<ParcelType>::initialise
 
     numberDensities /= nParticle_;
 
-    forAll(mesh_.cells(), cell)
+    forAll(mesh_.cells(), cellI)
     {
-        const vector& cC = mesh_.cellCentres()[cell];
-        const labelList& cellFaces = mesh_.cells()[cell];
-        const scalar cV = mesh_.cellVolumes()[cell];
+        List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
+        (
+            mesh_,
+            cellI
+        );
 
-        label nTets = 0;
-
-        // Each face is split into nEdges (or nVertices) - 2 tets.
-        forAll(cellFaces, face)
+        forAll(cellTets, tetI)
         {
-            nTets += mesh_.faces()[cellFaces[face]].size() - 2;
-        }
+            const tetIndices& cellTetIs = cellTets[tetI];
 
-        // Calculate the cumulative tet volumes circulating around the cell and
-        // record the vertex labels of each.
-        scalarList cTetVFracs(nTets, 0.0);
+            tetPointRef tet = cellTetIs.tet(mesh_);
 
-        List<labelList> tetPtIs(nTets, labelList(3,-1));
-
-        // Keep track of which tet this is.
-        label tet = 0;
-
-        forAll(cellFaces, face)
-        {
-            const labelList& facePoints = mesh_.faces()[cellFaces[face]];
-
-            label pointI = 1;
-            while ((pointI + 1) < facePoints.size())
-            {
-
-                const vector& pA = mesh_.points()[facePoints[0]];
-                const vector& pB = mesh_.points()[facePoints[pointI]];
-                const vector& pC = mesh_.points()[facePoints[pointI + 1]];
-
-                cTetVFracs[tet] =
-                    mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
-                  + cTetVFracs[max((tet - 1),0)];
-
-                tetPtIs[tet][0] = facePoints[0];
-                tetPtIs[tet][1] = facePoints[pointI];
-                tetPtIs[tet][2] = facePoints[pointI + 1];
-
-                pointI++;
-                tet++;
-            }
-        }
-
-        // Force the last volume fraction value to 1.0 to avoid any
-        // rounding/non-flat face errors giving a value < 1.0
-        cTetVFracs[nTets - 1] = 1.0;
+            scalar tetVolume = tet.mag();
 
             forAll(molecules, i)
             {
@@ -180,67 +144,25 @@ void Foam::DsmcCloud<ParcelType>::initialise
                 scalar numberDensity = numberDensities[i];
 
                 // Calculate the number of particles required
-            scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
+                scalar particlesRequired = numberDensity*tetVolume;
 
                 // Only integer numbers of particles can be inserted
                 label nParticlesToInsert = label(particlesRequired);
 
                 // Add another particle with a probability proportional to the
                 // remainder of taking the integer part of particlesRequired
-            if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
+                if
+                (
+                    (particlesRequired - nParticlesToInsert)
+                  > rndGen_.scalar01()
+                )
                 {
                     nParticlesToInsert++;
                 }
 
                 for (label pI = 0; pI < nParticlesToInsert; pI++)
                 {
-                // Choose a random point in a generic tetrahedron
-
-                scalar s = rndGen_.scalar01();
-                scalar t = rndGen_.scalar01();
-                scalar u = rndGen_.scalar01();
-
-                if (s + t > 1.0)
-                {
-                    s = 1.0 - s;
-                    t = 1.0 - t;
-                }
-
-                if (t + u > 1.0)
-                {
-                    scalar tmp = u;
-                    u = 1.0 - s - t;
-                    t = 1.0 - tmp;
-                }
-                else if (s + t + u > 1.0)
-                {
-                    scalar tmp = u;
-                    u = s + t + u - 1.0;
-                    s = 1.0 - t - tmp;
-                }
-
-                // Choose a tetrahedron to insert in, based on their relative
-                // volumes
-                scalar tetSelection = rndGen_.scalar01();
-
-                // Selected tetrahedron
-                label sTet = -1;
-
-                forAll(cTetVFracs, tet)
-                {
-                    sTet = tet;
-
-                    if (cTetVFracs[tet] >= tetSelection)
-                    {
-                        break;
-                    }
-                }
-
-                vector p =
-                    (1 - s - t - u)*cC
-                  + s*mesh_.points()[tetPtIs[sTet][0]]
-                  + t*mesh_.points()[tetPtIs[sTet][1]]
-                  + u*mesh_.points()[tetPtIs[sTet][2]];
+                    point p = tet.randomPoint(rndGen_);
 
                     vector U = equipartitionLinearVelocity
                     (
@@ -261,12 +183,15 @@ void Foam::DsmcCloud<ParcelType>::initialise
                         p,
                         U,
                         Ei,
-                    cell,
+                        cellI,
+                        cellTetIs.face(),
+                        cellTetIs.tetPt(),
                         typeId
                     );
                 }
             }
         }
+    }
 
     // Initialise the sigmaTcRMax_ field to the product of the cross section of
     // the most abundant species and the most probable thermal speed (Bird,
@@ -306,9 +231,9 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
     label collisions = 0;
 
-    forAll(cellOccupancy_, celli)
+    forAll(cellOccupancy_, cellI)
     {
-        const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
+        const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[cellI]);
 
         label nC(cellParcels.size());
 
@@ -327,13 +252,13 @@ void Foam::DsmcCloud<ParcelType>::collisions()
             // Inverse addressing specifying which subCell a parcel is in
             List<label> whichSubCell(cellParcels.size());
 
-            const point& cC = mesh_.cellCentres()[celli];
+            const point& cC = mesh_.cellCentres()[cellI];
 
             forAll(cellParcels, i)
             {
-                ParcelType* p = cellParcels[i];
+                const ParcelType& p = *cellParcels[i];
 
-                vector relPos = p->position() - cC;
+                vector relPos = p.position() - cC;
 
                 label subCell =
                     pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
@@ -345,16 +270,16 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-            scalar sigmaTcRMax = sigmaTcRMax_[celli];
+            scalar sigmaTcRMax = sigmaTcRMax_[cellI];
 
             scalar selectedPairs =
-                collisionSelectionRemainder_[celli]
+                collisionSelectionRemainder_[cellI]
               + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
-               /mesh_.cellVolumes()[celli];
+               /mesh_.cellVolumes()[cellI];
 
             label nCandidates(selectedPairs);
 
-            collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
+            collisionSelectionRemainder_[cellI] = selectedPairs - nCandidates;
 
             collisionCandidates += nCandidates;
 
@@ -415,39 +340,30 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-                ParcelType* parcelP = cellParcels[candidateP];
-                ParcelType* parcelQ = cellParcels[candidateQ];
-
-                label typeIdP = parcelP->typeId();
-                label typeIdQ = parcelQ->typeId();
+                ParcelType& parcelP = *cellParcels[candidateP];
+                ParcelType& parcelQ = *cellParcels[candidateQ];
 
                 scalar sigmaTcR = binaryCollision().sigmaTcR
                 (
-                    typeIdP,
-                    typeIdQ,
-                    parcelP->U(),
-                    parcelQ->U()
+                    parcelP,
+                    parcelQ
                 );
 
                 // Update the maximum value of sigmaTcR stored, but use the
                 // initial value in the acceptance-rejection criteria because
                 // the number of collision candidates selected was based on this
 
-                if (sigmaTcR > sigmaTcRMax_[celli])
+                if (sigmaTcR > sigmaTcRMax_[cellI])
                 {
-                    sigmaTcRMax_[celli] = sigmaTcR;
+                    sigmaTcRMax_[cellI] = sigmaTcR;
                 }
 
                 if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
                 {
                     binaryCollision().collide
                     (
-                        typeIdP,
-                        typeIdQ,
-                        parcelP->U(),
-                        parcelQ->U(),
-                        parcelP->Ei(),
-                        parcelQ->Ei()
+                        parcelP,
+                        parcelQ
                     );
 
                     collisions++;
@@ -577,7 +493,9 @@ void Foam::DsmcCloud<ParcelType>::addNewParcel
     const vector& position,
     const vector& U,
     const scalar Ei,
-    const label cellId,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
     const label typeId
 )
 {
@@ -587,7 +505,9 @@ void Foam::DsmcCloud<ParcelType>::addNewParcel
         position,
         U,
         Ei,
-        cellId,
+        cellI,
+        tetFaceI,
+        tetPtI,
         typeId
     );
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H

@@ -214,6 +214,10 @@ public:
     virtual ~DsmcCloud();
 
 
+    //- Type of parcel the cloud was instantiated for
+    typedef ParcelType parcelType;
+
+
     // Member Functions
 
         // Access
@@ -452,7 +456,9 @@ public:
                 const vector& position,
                 const vector& U,
                 const scalar Ei,
-                const label cellId,
+                const label cellI,
+                const label tetFaceI,
+                const label tetPtI,
                 const label typeId
             );
 

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C

@@ -92,7 +92,9 @@ bool Foam::DsmcParcel<ParcelType>::hitPatch
 (
     const polyPatch&,
     TrackData& td,
-    const label patchI
+    const label,
+    const scalar,
+    const tetIndices&
 )
 {
     return false;
@@ -125,7 +127,8 @@ template<class TrackData>
 void Foam::DsmcParcel<ParcelType>::hitWallPatch
 (
     const wallPolyPatch& wpp,
-    TrackData& td
+    TrackData& td,
+    const tetIndices& tetIs
 )
 {
     label wppIndex = wpp.index();
@@ -171,11 +174,8 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
 
     td.cloud().wallInteraction().correct
     (
-        wpp,
-        this->face(),
-        U_,
-        Ei_,
-        typeId_
+        static_cast<ParcelType&>(*this),
+        wpp
     );
 
     U_dot_nw = U_ & nw;
@@ -219,7 +219,8 @@ template<class ParcelType>
 void Foam::DsmcParcel<ParcelType>::hitWallPatch
 (
     const wallPolyPatch&,
-    int&
+    int&,
+    const tetIndices&
 )
 {}
 

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H

@@ -185,7 +185,9 @@ public:
             const vector& position,
             const vector& U,
             const scalar Ei,
-            const label celli,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
             const label typeId
         );
 
@@ -244,7 +246,9 @@ public:
             (
                 const polyPatch&,
                 TrackData& td,
-                const label patchI
+                const label patchI,
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 
             //- Overridable function to handle the particle hitting a
@@ -269,7 +273,8 @@ public:
             void hitWallPatch
             (
                 const wallPolyPatch&,
-                TrackData& td
+                TrackData& td,
+                const tetIndices&
             );
 
             //- Overridable function to handle the particle hitting a wallPatch
@@ -277,7 +282,8 @@ public:
             void hitWallPatch
             (
                 const wallPolyPatch&,
-                int&
+                int&,
+                const tetIndices&
             );
 
             //- Overridable function to handle the particle hitting a polyPatch
@@ -352,4 +358,3 @@ public:
 #endif
 
 // ************************************************************************* //
-

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H

@@ -69,11 +69,13 @@ inline Foam::DsmcParcel<ParcelType>::DsmcParcel
     const vector& position,
     const vector& U,
     const scalar Ei,
-    const label celli,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
     const label typeId
 )
 :
-    Particle<ParcelType>(owner, position, celli),
+    Particle<ParcelType>(owner, position, cellI, tetFaceI, tetPtI),
     U_(U),
     Ei_(Ei),
     typeId_(typeId)

src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C

@@ -43,7 +43,9 @@ Foam::dsmcParcel::dsmcParcel
     const vector& position,
     const vector& U,
     const scalar Ei,
-    const label celli,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
     const label typeId
 )
 :
@@ -53,7 +55,9 @@ Foam::dsmcParcel::dsmcParcel
         position,
         U,
         Ei,
-        celli,
+        cellI,
+        tetFaceI,
+        tetPtI,
         typeId
     )
 {}

src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H

@@ -66,7 +66,9 @@ public:
             const vector& position,
             const vector& U,
             const scalar Ei,
-            const label celli,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
             const label typeId
         );