Merge branch 'master' of ssh://graham@hunt//home/noisy3/OpenFOAM/OpenFOAM-dev

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

@@ -1,10 +1,265 @@
 #include "checkGeometry.H"
 #include "polyMesh.H"
-#include "globalMeshData.H"
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
 #include "EdgeMap.H"
+#include "wedgePolyPatch.H"
+#include "unitConversion.H"
+
+
+// Find wedge with opposite orientation. Note: does not actually check that
+// it is opposite, only that it has opposite normal and same axis
+Foam::label Foam::findOppositeWedge
+(
+    const polyMesh& mesh,
+    const wedgePolyPatch& wpp
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
+
+    forAll(patches, patchI)
+    {
+        if
+        (
+            patchI != wpp.index()
+         && patches[patchI].size()
+         && isA<wedgePolyPatch>(patches[patchI])
+        )
+        {
+            const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
+            (
+                patches[patchI]
+            );
+
+            // Calculate (cos of) angle to wpp (not pp!) centre normal
+            scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
+
+            if
+            (
+                pp.size() == wpp.size()
+             && mag(pp.axis() & wpp.axis()) >= (1-1E-3)
+             && mag(ppCosAngle - wppCosAngle) >= 1E-3
+            )
+            {
+                return patchI;
+            }
+        }
+    }
+    return -1;
+}
+
+
+bool Foam::checkWedges
+(
+    const polyMesh& mesh,
+    const bool report,
+    const Vector<label>& directions,
+    labelHashSet* setPtr
+)
+{
+    // To mark edges without calculating edge addressing
+    EdgeMap<label> edgesInError;
+
+    const pointField& p = mesh.points();
+    const faceList& fcs = mesh.faces();
+
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    forAll(patches, patchI)
+    {
+        if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
+        {
+            const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
+            (
+                patches[patchI]
+            );
+
+            scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
+
+            if (report)
+            {
+                Info<< "    Wedge " << pp.name() << " with angle "
+                    << radToDeg(wedgeAngle) << " degrees"
+                    << endl;                
+            }
+
+            // Find opposite
+            label oppositePatchI = findOppositeWedge(mesh, pp);
+
+            if (oppositePatchI == -1)
+            {
+                if (report)
+                {
+                    Info<< " ***Cannot find opposite wedge for wedge "
+                        << pp.name() << endl;
+                }
+                return true;
+            }
+
+            const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
+            (
+                patches[oppositePatchI]
+            );
+
+
+            if (mag(opp.axis() & pp.axis()) < (1-1E-3))
+            {
+                if (report)
+                {
+                    Info<< " ***Wedges do not have the same axis."
+                        << " Encountered " << pp.axis()
+                        << " on patch " << pp.name()
+                        << " which differs from " << opp.axis()
+                        << " on opposite wedge patch" << opp.axis()
+                        << endl;
+                }
+                return true;
+            }
+
+
+
+            // Mark edges on wedgePatches
+            forAll(pp, i)
+            {
+                const face& f = pp[i];
+                forAll(f, fp)
+                {
+                    label p0 = f[fp];
+                    label p1 = f.nextLabel(fp);
+                    edgesInError.insert(edge(p0, p1), -1);  // non-error value
+                }
+            }
+
+
+            // Check that wedge patch is flat
+            const point& p0 = p[pp.meshPoints()[0]];
+            forAll(pp.meshPoints(), i)
+            {
+                const point& pt = p[pp.meshPoints()[i]];
+                scalar d = mag((pt-p0) & pp.patchNormal());
+
+                if (d > sqrt(SMALL))
+                {
+                    if (report)
+                    {
+                        Info<< " ***Wedge patch " << pp.name() << " not planar."
+                            << " Point " << pt << " is not in patch plane by "
+                            << d << " meter."
+                            << endl;
+                    }
+                    return true;
+                }
+            }
+        }
+    }
+
+
+
+    // Check all non-wedge faces
+    label nEdgesInError = 0;
+
+    forAll(fcs, faceI)
+    {
+        const face& f = fcs[faceI];
+
+        forAll(f, fp)
+        {
+            label p0 = f[fp];
+            label p1 = f.nextLabel(fp);
+            if (p0 < p1)
+            {
+                vector d(p[p1]-p[p0]);
+                scalar magD = mag(d);
+
+                if (magD > ROOTVSMALL)
+                {
+                    d /= magD;
+
+                    // Check how many empty directions are used by the edge.
+                    label nEmptyDirs = 0;
+                    label nNonEmptyDirs = 0;
+                    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+                    {
+                        if (mag(d[cmpt]) > 1e-6)
+                        {
+                            if (directions[cmpt] == 0)
+                            {
+                                nEmptyDirs++;
+                            }
+                            else
+                            {
+                                nNonEmptyDirs++;
+                            }
+                        }
+                    }
+
+                    if (nEmptyDirs == 0)
+                    {
+                        // Purely in ok directions.
+                    }
+                    else if (nEmptyDirs == 1)
+                    {
+                        // Ok if purely in empty directions.
+                        if (nNonEmptyDirs > 0)
+                        {
+                            if (edgesInError.insert(edge(p0, p1), faceI))
+                            {
+                                nEdgesInError++;
+                            }
+                        }
+                    }
+                    else if (nEmptyDirs > 1)
+                    {
+                        // Always an error
+                        if (edgesInError.insert(edge(p0, p1), faceI))
+                        {
+                            nEdgesInError++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());
+
+    if (nErrorEdges > 0)
+    {
+        if (report)
+        {
+            Info<< " ***Number of edges not aligned with or perpendicular to "
+                << "non-empty directions: " << nErrorEdges << endl;
+        }
+
+        if (setPtr)
+        {
+            setPtr->resize(2*nEdgesInError);
+            forAllConstIter(EdgeMap<label>, edgesInError, iter)
+            {
+                if (iter() >= 0)
+                {
+                    setPtr->insert(iter.key()[0]);
+                    setPtr->insert(iter.key()[1]);
+                }
+            }
+        }
+
+        return true;
+    }
+    else
+    {
+        if (report)
+        {
+            Info<< "    All edges aligned with or perpendicular to "
+                << "non-empty directions." << endl;
+        }
+        return false;
+    }
+}
+
 
 Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
 {
@@ -33,7 +288,17 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     {
         pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
 
-        if (mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints))
+        if
+        (
+            (
+                validDirs != solDirs
+             && checkWedges(mesh, true, validDirs, &nonAlignedPoints)
+            )
+         || (
+                validDirs == solDirs
+             && mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
+            )
+        )
         {
             noFailedChecks++;
             label nNonAligned = returnReduce
@@ -58,7 +323,16 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     {
         cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
         cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
-        if (mesh.checkClosedCells(true, &cells, &aspectCells))
+        if
+        (
+            mesh.checkClosedCells
+            (
+                true,
+                &cells,
+                &aspectCells,
+                mesh.geometricD()
+            )
+        )
         {
             noFailedChecks++;
 
@@ -86,7 +360,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     }
 
     {
-        faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100 + 1);
+        faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100+1);
         if (mesh.checkFaceAreas(true, &faces))
         {
             noFailedChecks++;
@@ -104,7 +378,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     }
 
     {
-        cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100 + 1);
+        cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100+1);
         if (mesh.checkCellVolumes(true, &cells))
         {
             noFailedChecks++;
@@ -122,7 +396,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     }
 
     {
-        faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100 + 1);
+        faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100+1);
         if (mesh.checkFaceOrthogonality(true, &faces))
         {
             noFailedChecks++;
@@ -160,7 +434,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     }
 
     {
-        faceSet faces(mesh, "wrongOrientedTriangleFaces", mesh.nFaces()/100 + 1);
+        faceSet faces(mesh, "wrongOrientedTriangleFaces", mesh.nFaces()/100+1);
         if (mesh.checkFaceTets(true, 0, &faces))
         {
             noFailedChecks++;
@@ -179,7 +453,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     }
 
     {
-        faceSet faces(mesh, "skewFaces", mesh.nFaces()/100 + 1);
+        faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
         if (mesh.checkFaceSkewness(true, &faces))
         {
             noFailedChecks++;
@@ -278,7 +552,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     if (allGeometry)
     {
         cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
-        if (mesh.checkCellDeterminant(true, &cells))
+        if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
         {
             noFailedChecks++;
 

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

@@ -1,8 +1,21 @@
 #include "label.H"
+#include "HashSet.H"
+#include "labelVector.H"
 
 namespace Foam
 {
     class polyMesh;
+    class wedgePolyPatch;
+
+    label findOppositeWedge(const polyMesh&, const wedgePolyPatch&);
+
+    bool checkWedges
+    (
+        const polyMesh&,
+        const bool report,
+        const Vector<label>&,
+        labelHashSet*
+    );
 
     label checkGeometry(const polyMesh& mesh, const bool allGeometry);
 }

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -544,7 +544,8 @@ public:
                 (
                     const bool report = false,
                     labelHashSet* setPtr = NULL,
-                    labelHashSet* highAspectSetPtr = NULL
+                    labelHashSet* highAspectSetPtr = NULL,
+                    const Vector<label>& solutionD = Vector<label>::one
                 ) const;
 
                 //- Check for negative face areas
@@ -645,7 +646,8 @@ public:
                 bool checkCellDeterminant
                 (
                     const bool report = false,
-                    labelHashSet* setPtr = NULL
+                    labelHashSet* setPtr = NULL,
+                    const Vector<label>& solutionD = Vector<label>::one
                 ) const;
 
                 //- Check for concave cells by the planes of faces

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,13 +98,15 @@ bool Foam::primitiveMesh::checkClosedCells
 (
     const bool report,
     labelHashSet* setPtr,
-    labelHashSet* aspectSetPtr
+    labelHashSet* aspectSetPtr,
+    const Vector<label>& meshD
 ) const
 {
     if (debug)
     {
         Info<< "bool primitiveMesh::checkClosedCells("
-            << "const bool, labelHashSet*, labelHashSet*) const: "
+            << "const bool, labelHashSet*, labelHashSet*"
+            << ", const Vector<label>&) const: "
             << "checking whether cells are closed" << endl;
     }
 
@@ -171,6 +173,16 @@ bool Foam::primitiveMesh::checkClosedCells
 
     const scalarField& vols = cellVolumes();
 
+    label nDims = 0;
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (meshD[dir] == 1)
+        {
+            nDims++;
+        }
+    }
+
+
     // Check the sums
     forAll(sumClosed, cellI)
     {
@@ -200,12 +212,26 @@ bool Foam::primitiveMesh::checkClosedCells
 
         // Calculate the aspect ration as the maximum of Cartesian component
         // aspect ratio to the total area hydraulic area aspect ratio
-        scalar aspectRatio = max
+        scalar minCmpt = VGREAT;
-        (
+        scalar maxCmpt = -VGREAT;
-            cmptMax(sumMagClosed[cellI])
+        for (direction dir = 0; dir < vector::nComponents; dir++)
-           /(cmptMin(sumMagClosed[cellI]) + VSMALL),
+        {
-            1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
+            if (meshD[dir] == 1)
-        );
+            {
+                minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
+                maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
+            }
+        }
+
+        scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
+        if (nDims == 3)
+        {
+            aspectRatio = max
+            (
+                aspectRatio,
+                1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
+            );
+        }
 
         maxAspectRatio = max(maxAspectRatio, aspectRatio);
 
@@ -1977,7 +2003,8 @@ bool Foam::primitiveMesh::checkFaceFaces
 bool Foam::primitiveMesh::checkCellDeterminant
 (
     const bool report,    // report,
-    labelHashSet* setPtr  // setPtr
+    labelHashSet* setPtr, // setPtr
+    const Vector<label>& meshD
 ) const
 {
     if (debug)
@@ -1987,6 +2014,22 @@ bool Foam::primitiveMesh::checkCellDeterminant
             << "checking for under-determined cells" << endl;
     }
 
+    // Determine number of dimensions and (for 2D) missing dimension
+    label nDims = 0;
+    label twoD = -1;
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (meshD[dir] == 1)
+        {
+            nDims++;
+        }
+        else
+        {
+            twoD = dir;
+        }
+    }
+
+
     const cellList& c = cells();
 
     label nErrorCells = 0;
@@ -1995,55 +2038,34 @@ bool Foam::primitiveMesh::checkCellDeterminant
     scalar sumDet = 0;
     label nSummed = 0;
 
-    forAll(c, cellI)
+    if (nDims == 1)
     {
-        const labelList& curFaces = c[cellI];
+        minDet = 1;
-
+        sumDet = c.size()*minDet;
-        // Calculate local normalization factor
+        nSummed = c.size();
-        scalar avgArea = 0;
+    }
-
+    else
-        label nInternalFaces = 0;
+    {
-
+        forAll (c, cellI)
-        forAll(curFaces, i)
         {
-            if (isInternalFace(curFaces[i]))
+            const labelList& curFaces = c[cellI];
-            {
-                avgArea += mag(faceAreas()[curFaces[i]]);
 
-                nInternalFaces++;
+            // Calculate local normalization factor
-            }
+            scalar avgArea = 0;
-        }
 
-        if (nInternalFaces == 0)
+            label nInternalFaces = 0;
-        {
-            if (setPtr)
-            {
-                setPtr->insert(cellI);
-            }
-
-            nErrorCells++;
-        }
-        else
-        {
-            avgArea /= nInternalFaces;
-
-            symmTensor areaTensor(symmTensor::zero);
 
             forAll(curFaces, i)
             {
                 if (isInternalFace(curFaces[i]))
                 {
-                    areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
+                    avgArea += mag(faceAreas()[curFaces[i]]);
+
+                    nInternalFaces++;
                 }
             }
 
-            scalar determinant = mag(det(areaTensor));
+            if (nInternalFaces == 0)
-
-            minDet = min(determinant, minDet);
-            sumDet += determinant;
-            nSummed++;
-
-            if (determinant < 1e-3)
             {
                 if (setPtr)
                 {
@@ -2052,6 +2074,54 @@ bool Foam::primitiveMesh::checkCellDeterminant
 
                 nErrorCells++;
             }
+            else
+            {
+                avgArea /= nInternalFaces;
+
+                symmTensor areaTensor(symmTensor::zero);
+
+                forAll(curFaces, i)
+                {
+                    if (isInternalFace(curFaces[i]))
+                    {
+                        areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
+                    }
+                }
+
+                if (nDims == 2)
+                {
+                    // Add the missing eigenvector (such that it does not
+                    // affect the determinant)
+                    if (twoD == 0)
+                    {
+                        areaTensor.xx() = 1;
+                    }
+                    else if (twoD == 1)
+                    {
+                        areaTensor.yy() = 1;
+                    }
+                    else
+                    {
+                        areaTensor.zz() = 1;
+                    }
+                }
+
+                scalar determinant = mag(det(areaTensor));
+
+                minDet = min(determinant, minDet);
+                sumDet += determinant;
+                nSummed++;
+
+                if (determinant < 1e-3)
+                {
+                    if (setPtr)
+                    {
+                        setPtr->insert(cellI);
+                    }
+
+                    nErrorCells++;
+                }
+            }
         }
     }
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,13 +30,14 @@ License
 void Foam::setRefCell
 (
     const volScalarField& field,
+    const volScalarField& fieldRef,
     const dictionary& dict,
     label& refCelli,
     scalar& refValue,
     const bool forceReference
 )
 {
-    if (field.needReference() || forceReference)
+    if (fieldRef.needReference() || forceReference)
     {
         word refCellName = field.name() + "RefCell";
         word refPointName = field.name() + "RefPoint";
@@ -56,6 +57,7 @@ void Foam::setRefCell
                         "void Foam::setRefCell\n"
                          "(\n"
                          "    const volScalarField&,\n"
+                         "    const volScalarField&,\n"
                          "    const dictionary&,\n"
                          "    label& scalar&,\n"
                          "    bool\n"
@@ -76,7 +78,7 @@ void Foam::setRefCell
             point refPointi(dict.lookup(refPointName));
             refCelli = field.mesh().findCell(refPointi);
             label hasRef = (refCelli >= 0 ? 1 : 0);
-            label sumHasRef = returnReduce(hasRef, sumOp<label>());
+            label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
             if (sumHasRef != 1)
             {
                 FatalIOErrorIn
@@ -84,6 +86,7 @@ void Foam::setRefCell
                     "void Foam::setRefCell\n"
                      "(\n"
                      "    const volScalarField&,\n"
+                     "    const volScalarField&,\n"
                      "    const dictionary&,\n"
                      "    label& scalar&,\n"
                      "    bool\n"
@@ -103,6 +106,7 @@ void Foam::setRefCell
                 "void Foam::setRefCell\n"
                  "(\n"
                  "    const volScalarField&,\n"
+                 "    const volScalarField&,\n"
                  "    const dictionary&,\n"
                  "    label& scalar&,\n"
                  "    bool\n"
@@ -119,6 +123,19 @@ void Foam::setRefCell
 }
 
 
+void Foam::setRefCell
+(
+    const volScalarField& field,
+    const dictionary& dict,
+    label& refCelli,
+    scalar& refValue,
+    const bool forceReference
+)
+{
+    setRefCell(field, field, dict, refCelli, refValue, forceReference);
+}
+
+
 Foam::scalar Foam::getRefCellValue
 (
     const volScalarField& field,

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,8 +44,22 @@ SourceFiles
 namespace Foam
 {
 
-//- Find the reference cell nearest (in index) to the given cell,
+//- If the field fieldRef needs referencing find the reference cell nearest
-//  but which is not on a cyclic, symmetry or processor patch.
+//  (in index) to the given cell looked-up for field, but which is not on a
+//  cyclic, symmetry or processor patch.
+void setRefCell
+(
+    const volScalarField& field,
+    const volScalarField& fieldRef,
+    const dictionary& dict,
+    label& refCelli,
+    scalar& refValue,
+    const bool forceReference = false
+);
+
+//- If the field needs referencing find the reference cell nearest
+//  (in index) to the given cell looked-up for field, but which is not on a
+//  cyclic, symmetry or processor patch.
 void setRefCell
 (
     const volScalarField& field,

src/lagrangian/basic/Particle/Particle.C

@@ -334,11 +334,18 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 position_ = endPosition;
             }
 
+            label origFacei = facei_;
             label patchi = patch(facei_);
-            const polyPatch& patch = mesh.boundaryMesh()[patchi];
 
-            if (!p.hitPatch(patch, td, patchi))
+            if (!p.hitPatch(mesh.boundaryMesh()[patchi], td, patchi))
             {
+                // Did patch interaction model switch patches?
+                if (facei_ != origFacei)
+                {
+                    patchi = patch(facei_);
+                }
+                const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
                 if (isA<wedgePolyPatch>(patch))
                 {
                     p.hitWedgePatch

src/lagrangian/basic/Particle/Particle.H

@@ -378,6 +378,9 @@ public:
             //- Return current cell particle is in
             inline label cell() const;
 
+            //- Return current face particle is on otherwise -1
+            inline label& face();
+
             //- Return current face particle is on otherwise -1
             inline label face() const;
 

src/lagrangian/basic/Particle/ParticleI.H

@@ -325,6 +325,13 @@ inline Foam::label Foam::Particle<ParticleType>::face() const
 }
 
 
+template<class ParticleType>
+inline Foam::label& Foam::Particle<ParticleType>::face()
+{
+    return facei_;
+}
+
+
 template<class ParticleType>
 inline bool Foam::Particle<ParticleType>::onBoundary() const
 {

src/lagrangian/intermediate/Make/options

@@ -31,5 +31,8 @@ LIB_LIBS = \
     -lspecie \
     -lbasicThermophysicalModels \
     -lreactionThermophysicalModels \
+    -lchemistryModel \
+    -lradiation \
+    -lODE \
     -lcompressibleRASModels \
     -lcompressibleLESModels