From 6061cf971c87ce9f207d471cfcf074be79cd71b1 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 26 Apr 2012 17:14:10 +0100
Subject: [PATCH 1/5] ENH: polyTopoChange: override readOpt flag in makeMesh;
 does not make sense to read when constructing new mesh

---
 .../polyTopoChange/polyTopoChange/polyTopoChange.C            | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index 9000380aa1..b9b7bde14b 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -3394,11 +3394,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     // Create the mesh
     // ~~~~~~~~~~~~~~~
 
+    IOobject noReadIO(io);
+    noReadIO.readOpt() = IOobject::NO_READ;
     newMeshPtr.reset
     (
         new fvMesh
         (
-            io,
+            noReadIO,
             xferMove(newPoints),
             faces_.xfer(),
             faceOwner_.xfer(),

From b1d84337c7080784019fcd7cbaecd464deef7145 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 27 Apr 2012 11:38:40 +0100
Subject: [PATCH 2/5] ENH: tetDecomposer: class to split cells into tets

---
 src/dynamicMesh/Make/files                    |   1 +
 .../polyTopoChange/tetDecomposer.C            | 713 ++++++++++++++++++
 .../polyTopoChange/tetDecomposer.H            | 198 +++++
 3 files changed, 912 insertions(+)
 create mode 100644 src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
 create mode 100644 src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H

diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index a7f72ef9fa..d7310029a7 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -37,6 +37,7 @@ polyTopoChange/polyTopoChange/removePoints.C
 polyTopoChange/polyTopoChange/combineFaces.C
 polyTopoChange/polyTopoChange/localPointRegion.C
 polyTopoChange/polyTopoChange/duplicatePoints.C
+polyTopoChange/polyTopoChange/tetDecomposer.C
 
 slidingInterface/slidingInterface.C
 slidingInterface/slidingInterfaceProjectPoints.C
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
new file mode 100644
index 0000000000..e8948c4fad
--- /dev/null
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
@@ -0,0 +1,713 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetDecomposer.H"
+#include "meshTools.H"
+#include "polyMesh.H"
+#include "polyTopoChange.H"
+#include "mapPolyMesh.H"
+#include "OFstream.H"
+#include "EdgeMap.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::tetDecomposer, 0);
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum<tetDecomposer::decompositionType, 2>::names[] =
+    {
+        "faceCentre",
+        "faceDiagonal"
+    };
+
+    const NamedEnum<tetDecomposer::decompositionType, 2>
+        tetDecomposer::decompositionTypeNames;
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::tetDecomposer::modifyFace
+(
+    polyTopoChange& meshMod,
+    const face& f,
+    const label faceI,
+    const label own,
+    const label nei,
+    const label patchI,
+    const label zoneI,
+    const bool zoneFlip
+) const
+{
+    // First usage of face. Modify.
+    if (nei == -1 || own < nei)
+    {
+        meshMod.modifyFace
+        (
+            f,                          // modified face
+            faceI,                      // label of face
+            own,                        // owner
+            nei,                        // neighbour
+            false,                      // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            zoneFlip                    // face flip in zone
+        );
+    }
+    else
+    {
+        meshMod.modifyFace
+        (
+            f.reverseFace(),            // modified face
+            faceI,                      // label of face
+            nei,                        // owner
+            own,                        // neighbour
+            true,                       // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            !zoneFlip                   // face flip in zone
+        );
+    }
+}
+
+
+void Foam::tetDecomposer::addFace
+(
+    polyTopoChange& meshMod,
+    const face& f,
+    const label own,
+    const label nei,
+    const label masterPointID,
+    const label masterEdgeID,
+    const label masterFaceID,
+    const label patchI,
+    const label zoneI,
+    const bool zoneFlip
+) const
+{
+    // Second or more usage of face. Add.
+    if (nei == -1 || own < nei)
+    {
+        meshMod.addFace
+        (
+            f,                          // modified face
+            own,                        // owner
+            nei,                        // neighbour
+            masterPointID,              // master point
+            masterEdgeID,               // master edge
+            masterFaceID,               // master face
+            false,                      // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            zoneFlip                    // face flip in zone
+        );
+    }
+    else
+    {
+        meshMod.addFace
+        (
+            f.reverseFace(),            // modified face
+            nei,                        // owner
+            own,                        // neighbour
+            masterPointID,              // master point
+            masterEdgeID,               // master edge
+            masterFaceID,               // master face
+            true,                       // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            !zoneFlip                   // face flip in zone
+        );
+    }
+}
+
+
+// Work out triangle index given the starting vertex in the face
+Foam::label Foam::tetDecomposer::triIndex(const label faceI, const label fp)
+const
+{
+    const face& f = mesh_.faces()[faceI];
+    const label fp0 = mesh_.tetBasePtIs()[faceI];
+
+    // Work out triangle index on this face
+    label thisTriI;
+    if (fp == fp0)
+    {
+        thisTriI = 0;
+    }
+    else if (fp == f.rcIndex(fp0))
+    {
+        thisTriI = f.size()-3;
+    }
+    else
+    {
+        thisTriI = (fp-fp0-1) % (f.size()-2);
+    }
+    return thisTriI;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::tetDecomposer::tetDecomposer(const polyMesh& mesh)
+:
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::tetDecomposer::setRefinement
+(
+    const decompositionType decomposeType,
+    polyTopoChange& meshMod
+)
+{
+    cellToPoint_.setSize(mesh_.nCells());
+    forAll(mesh_.cellCentres(), cellI)
+    {
+        // Any point on the cell
+        label masterPointI = mesh_.faces()[mesh_.cells()[cellI][0]][0];
+
+        cellToPoint_[cellI] = meshMod.addPoint
+        (
+            mesh_.cellCentres()[cellI],
+            masterPointI,
+            -1,
+            true
+        );
+    }
+
+
+    // Add face centre points
+    if (decomposeType == FACECENTRETETS)
+    {
+        faceToPoint_.setSize(mesh_.nFaces());
+        forAll(mesh_.faceCentres(), faceI)
+        {
+            // Any point on the face
+            const label masterPointI = mesh_.faces()[faceI][0];
+
+            faceToPoint_[faceI] = meshMod.addPoint
+            (
+                mesh_.faceCentres()[faceI],
+                masterPointI,
+                -1,
+                true
+            );
+        }
+    }
+
+
+    // Per face, per point (faceCentre) or triangle (faceDiag) the added cell
+    faceOwnerCells_.setSize(mesh_.nFaces());
+    faceNeighbourCells_.setSize(mesh_.nFaces());
+
+    if (decomposeType == FACECENTRETETS)
+    {
+        forAll(faceOwnerCells_, faceI)
+        {
+            const face& f = mesh_.faces()[faceI];
+            faceOwnerCells_[faceI].setSize(f.size(), -1);
+            faceNeighbourCells_[faceI].setSize(f.size(), -1);
+        }
+    }
+    else
+    {
+        // Force construction of diagonal decomposition
+        (void)mesh_.tetBasePtIs();
+
+        forAll(faceOwnerCells_, faceI)
+        {
+            const face& f = mesh_.faces()[faceI];
+            faceOwnerCells_[faceI].setSize(f.size()-2, -1);
+            faceNeighbourCells_[faceI].setSize(f.size()-2, -1);
+        }
+    }
+
+
+    forAll(mesh_.cells(), cellI)
+    {
+        const cell& cFaces = mesh_.cells()[cellI];
+
+        EdgeMap<label> edgeToFace(8*cFaces.size());
+
+        forAll(cFaces, cFaceI)
+        {
+            label faceI = cFaces[cFaceI];
+            const face& f = mesh_.faces()[faceI];
+
+            // Get reference to either owner or neighbour
+            labelList& added =
+            (
+                (mesh_.faceOwner()[faceI] == cellI)
+              ? faceOwnerCells_[faceI]
+              : faceNeighbourCells_[faceI]
+            );
+
+            if (decomposeType == FACECENTRETETS)
+            {
+                forAll(f, fp)
+                {
+                    if (cFaceI == 0 && fp == 0)
+                    {
+                        // Reuse cell itself
+                        added[fp] = cellI;
+                    }
+                    else
+                    {
+                        added[fp] = meshMod.addCell
+                        (
+                            -1,     // masterPoint
+                            -1,     // masterEdge
+                            -1,     // masterFace
+                            cellI,  // masterCell
+                            mesh_.cellZones().whichZone(cellI)
+                        );
+                    }
+                }
+            }
+            else
+            {
+                for (label triI = 0; triI < f.size()-2; triI++)
+                {
+                    if (cFaceI == 0 && triI == 0)
+                    {
+                        // Reuse cell itself
+                        added[triI] = cellI;
+                    }
+                    else
+                    {
+                        added[triI] = meshMod.addCell
+                        (
+                            -1,     // masterPoint
+                            -1,     // masterEdge
+                            -1,     // masterFace
+                            cellI,  // masterCell
+                            mesh_.cellZones().whichZone(cellI)
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // Add triangle faces
+    face triangle(3);
+
+    forAll(mesh_.faces(), faceI)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        const labelList& addedOwn = faceOwnerCells_[faceI];
+        const labelList& addedNei = faceNeighbourCells_[faceI];
+        const face& f = mesh_.faces()[faceI];
+
+        label patchI = -1;
+        if (faceI >= mesh_.nInternalFaces())
+        {
+            patchI = mesh_.boundaryMesh().whichPatch(faceI);
+        }
+
+        label zoneI = mesh_.faceZones().whichZone(faceI);
+        bool zoneFlip = false;
+        if (zoneI != -1)
+        {
+            const faceZone& fz = mesh_.faceZones()[zoneI];
+            zoneFlip = fz.flipMap()[fz.whichFace(faceI)];
+        }
+
+
+        if (decomposeType == FACECENTRETETS)
+        {
+            forAll(f, fp)
+            {
+                // 1. Front triangle (decomposition of face itself)
+                //    (between owner and neighbour cell)
+                {
+                    triangle[0] = f[fp];
+                    triangle[1] = f[f.fcIndex(fp)];
+                    triangle[2] = faceToPoint_[faceI];
+
+                    if (fp == 0)
+                    {
+                        modifyFace
+                        (
+                            meshMod,
+                            triangle,
+                            faceI,
+                            addedOwn[fp],
+                            addedNei[fp],
+                            patchI,
+                            zoneI,
+                            zoneFlip
+                        );
+                    }
+                    else
+                    {
+                        addFace
+                        (
+                            meshMod,
+                            triangle,
+                            addedOwn[fp],
+                            addedNei[fp],
+                            -1,                 //point
+                            -1,                 //edge
+                            faceI,              //face
+                            patchI,
+                            zoneI,
+                            zoneFlip
+                        );
+                    }
+                }
+
+
+                // 2. Within owner cell - to cell centre
+                {
+                    label newOwn = addedOwn[f.rcIndex(fp)];
+                    label newNei = addedOwn[fp];
+
+                    triangle[0] = f[fp];
+                    triangle[1] = cellToPoint_[own];
+                    triangle[2] = faceToPoint_[faceI];
+
+                    addFace
+                    (
+                        meshMod,
+                        triangle,
+                        newOwn,
+                        newNei,
+                        f[fp],      //point
+                        -1,         //edge
+                        -1,         //face
+                        -1,         //patchI
+                        zoneI,
+                        zoneFlip
+                    );
+                }
+                // 2b. Within neighbour cell - to cell centre
+                if (faceI < mesh_.nInternalFaces())
+                {
+                    label newOwn = addedNei[f.rcIndex(fp)];
+                    label newNei = addedNei[fp];
+
+                    triangle[0] = f[fp];
+                    triangle[1] = faceToPoint_[faceI];
+                    triangle[2] = cellToPoint_[mesh_.faceNeighbour()[faceI]];
+
+                    addFace
+                    (
+                        meshMod,
+                        triangle,
+                        newOwn,
+                        newNei,
+                        f[fp],      //point
+                        -1,         //edge
+                        -1,         //face
+                        -1,         //patchI
+                        zoneI,
+                        zoneFlip
+                    );
+                }
+            }
+        }
+        else
+        {
+            label fp0 = mesh_.tetBasePtIs()[faceI];
+            label fp = f.fcIndex(fp0);
+
+            for (label triI = 0; triI < f.size()-2; triI++)
+            {
+                label nextTri = triI+1;
+                if (nextTri >= f.size()-2)
+                {
+                    nextTri -= f.size()-2;
+                }
+                label nextFp = f.fcIndex(fp);
+
+
+                // Triangle triI consisiting of f[fp0], f[fp], f[nextFp]
+
+
+                // 1. Front triangle (decomposition of face itself)
+                //    (between owner and neighbour cell)
+                {
+                    triangle[0] = f[fp0];
+                    triangle[1] = f[fp];
+                    triangle[2] = f[nextFp];
+
+                    if (triI == 0)
+                    {
+                        modifyFace
+                        (
+                            meshMod,
+                            triangle,
+                            faceI,
+                            addedOwn[triI],
+                            addedNei[triI],
+                            patchI,
+                            zoneI,
+                            zoneFlip
+                        );
+                    }
+                    else
+                    {
+                        addFace
+                        (
+                            meshMod,
+                            triangle,
+                            addedOwn[triI],
+                            addedNei[triI],
+                            -1,                 //point
+                            -1,                 //edge
+                            faceI,              //face
+                            patchI,
+                            zoneI,
+                            zoneFlip
+                        );
+                    }
+                }
+
+
+                // 2. Within owner cell - diagonal to cell centre
+                if (triI < f.size()-3)
+                {
+                    label newOwn = addedOwn[triI];
+                    label newNei = addedOwn[nextTri];
+
+                    triangle[0] = f[fp0];
+                    triangle[1] = f[nextFp];
+                    triangle[2] = cellToPoint_[own];
+
+                    addFace
+                    (
+                        meshMod,
+                        triangle,
+                        newOwn,
+                        newNei,
+                        f[fp],      //point
+                        -1,         //edge
+                        -1,         //face
+                        -1,         //patchI
+                        zoneI,
+                        zoneFlip
+                    );
+
+                    // 2b. Within neighbour cell - to cell centre
+                    if (faceI < mesh_.nInternalFaces())
+                    {
+                        label newOwn = addedNei[triI];
+                        label newNei = addedNei[nextTri];
+
+                        triangle[0] = f[nextFp];
+                        triangle[1] = f[fp0];
+                        triangle[2] =
+                            cellToPoint_[mesh_.faceNeighbour()[faceI]];
+
+                        addFace
+                        (
+                            meshMod,
+                            triangle,
+                            newOwn,
+                            newNei,
+                            f[fp],      //point
+                            -1,         //edge
+                            -1,         //face
+                            -1,         //patchI
+                            zoneI,
+                            zoneFlip
+                        );
+                    }
+                }
+
+
+                fp = nextFp;
+            }
+        }
+    }
+
+
+
+    // Add triangles for all edges.
+    EdgeMap<label> edgeToFace;
+
+    forAll(mesh_.cells(), cellI)
+    {
+        const cell& cFaces = mesh_.cells()[cellI];
+
+        edgeToFace.clear();
+
+        forAll(cFaces, cFaceI)
+        {
+            label faceI = cFaces[cFaceI];
+
+            label zoneI = mesh_.faceZones().whichZone(faceI);
+            bool zoneFlip = false;
+            if (zoneI != -1)
+            {
+                const faceZone& fz = mesh_.faceZones()[zoneI];
+                zoneFlip = fz.flipMap()[fz.whichFace(faceI)];
+            }
+
+            const face& f = mesh_.faces()[faceI];
+            //const labelList& fEdges = mesh_.faceEdges()[faceI];
+            forAll(f, fp)
+            {
+                label p0 = f[fp];
+                label p1 = f[f.fcIndex(fp)];
+                const edge e(p0, p1);
+
+                EdgeMap<label>::const_iterator edgeFnd = edgeToFace.find(e);
+                if (edgeFnd == edgeToFace.end())
+                {
+                    edgeToFace.insert(e, faceI);
+                }
+                else
+                {
+                    // Found the other face on the edge.
+                    label otherFaceI = edgeFnd();
+                    const face& otherF = mesh_.faces()[otherFaceI];
+
+                    // Found the other face on the edge. Note that since
+                    // we are looping in the same order the tets added for
+                    // otherFaceI will be before those of faceI
+
+                    label otherFp = findIndex(otherF, p0);
+                    if (otherF.nextLabel(otherFp) == p1)
+                    {
+                        // ok. otherFp is first vertex of edge.
+                    }
+                    else if (otherF.prevLabel(otherFp) == p1)
+                    {
+                        otherFp = otherF.rcIndex(otherFp);
+                    }
+                    else
+                    {
+                        FatalErrorIn("tetDecomposer::setRefinement(..)")
+                            << "problem." << abort(FatalError);
+                    }
+
+
+                    // Triangle from edge to cell centre
+                    if (mesh_.faceOwner()[faceI] == cellI)
+                    {
+                        triangle[0] = p0;
+                        triangle[1] = p1;
+                        triangle[2] = cellToPoint_[cellI];
+                    }
+                    else
+                    {
+                        triangle[0] = p1;
+                        triangle[1] = p0;
+                        triangle[2] = cellToPoint_[cellI];
+                    }
+
+                    // Determine tets on either side
+                    label thisTet, otherTet;
+
+                    if (decomposeType == FACECENTRETETS)
+                    {
+                        if (mesh_.faceOwner()[faceI] == cellI)
+                        {
+                            thisTet = faceOwnerCells_[faceI][fp];
+                        }
+                        else
+                        {
+                            thisTet = faceNeighbourCells_[faceI][fp];
+                        }
+
+                        if (mesh_.faceOwner()[otherFaceI] == cellI)
+                        {
+                            otherTet = faceOwnerCells_[otherFaceI][otherFp];
+                        }
+                        else
+                        {
+                            otherTet =
+                                faceNeighbourCells_[otherFaceI][otherFp];
+                        }
+                    }
+                    else
+                    {
+                        label thisTriI = triIndex(faceI, fp);
+                        if (mesh_.faceOwner()[faceI] == cellI)
+                        {
+                            thisTet = faceOwnerCells_[faceI][thisTriI];
+                        }
+                        else
+                        {
+                            thisTet = faceNeighbourCells_[faceI][thisTriI];
+                        }
+
+                        label otherTriI = triIndex(otherFaceI, otherFp);
+                        if (mesh_.faceOwner()[otherFaceI] == cellI)
+                        {
+                            otherTet = faceOwnerCells_[otherFaceI][otherTriI];
+                        }
+                        else
+                        {
+                            otherTet =
+                                faceNeighbourCells_[otherFaceI][otherTriI];
+                        }
+                    }
+
+
+                    addFace
+                    (
+                        meshMod,
+                        triangle,
+                        otherTet,
+                        thisTet,
+                        -1,         //masterPoint
+                        -1,         //fEdges[fp], //masterEdge
+                        faceI,      //masterFace
+                        -1,         //patchI
+                        zoneI,
+                        zoneFlip
+                    );
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::tetDecomposer::updateMesh(const mapPolyMesh& map)
+{
+    inplaceRenumber(map.pointMap(), cellToPoint_);
+    inplaceRenumber(map.pointMap(), faceToPoint_);
+
+    forAll(faceOwnerCells_, faceI)
+    {
+        inplaceRenumber(map.cellMap(), faceOwnerCells_[faceI]);
+    }
+    forAll(faceNeighbourCells_, faceI)
+    {
+        inplaceRenumber(map.cellMap(), faceNeighbourCells_[faceI]);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H
new file mode 100644
index 0000000000..8485f3cd8d
--- /dev/null
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H
@@ -0,0 +1,198 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    tetDecomposer
+
+Description
+    Decomposes polyMesh into tets.
+
+SourceFiles
+    tetDecomposer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetDecomposer_H
+#define tetDecomposer_H
+
+#include "DynamicList.H"
+#include "PackedBoolList.H"
+#include "boolList.H"
+#include "typeInfo.H"
+#include "NamedEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class polyTopoChange;
+class face;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class tetDecomposer Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetDecomposer
+{
+public:
+
+    // Public data types
+
+        enum decompositionType
+        {
+            FACECENTRETETS, // tet decomposition using facectr and cellctr
+            FACEDIAGTETS    // tet decomposition using face diag and cellctr
+        };
+        static const NamedEnum<decompositionType, 2> decompositionTypeNames;
+
+private:
+
+    // Private data
+
+        const polyMesh& mesh_;
+
+        //- From cell to tet point
+        labelList cellToPoint_;
+
+        //- From face to tet point
+        labelList faceToPoint_;
+
+
+        // Per face, per point (faceCentre) or triangle (faceDiag)
+        // the added tet on the owner side
+        labelListList faceOwnerCells_;
+
+        // Per face, per point (faceCentre) or triangle (faceDiag)
+        // the added tet on the neighbour side
+        labelListList faceNeighbourCells_;
+
+    // Private Member Functions
+
+        //- Modify a face
+        void modifyFace
+        (
+            polyTopoChange& meshMod,
+            const face& f,
+            const label faceI,
+            const label own,
+            const label nei,
+            const label patchI,
+            const label zoneI,
+            const bool zoneFlip
+        ) const;
+
+        //- Add a face
+        void addFace
+        (
+            polyTopoChange& meshMod,
+            const face& f,
+            const label own,
+            const label nei,
+            const label masterPointID,
+            const label masterEdgeID,
+            const label masterFaceID,
+            const label patchI,
+            const label zoneI,
+            const bool zoneFlip
+        ) const;
+
+        //- Work out triangle index given the starting vertex in the face
+        label triIndex(const label faceI, const label fp) const;
+
+        //- Disallow default bitwise copy construct
+        tetDecomposer(const tetDecomposer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const tetDecomposer&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("tetDecomposer");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        tetDecomposer(const polyMesh&);
+
+
+    // Member Functions
+
+        // Access
+
+            //- From cell to tet point
+            const labelList& cellToPoint() const
+            {
+                return cellToPoint_;
+            }
+
+            //- From face to tet point
+            const labelList& faceToPoint() const
+            {
+                return faceToPoint_;
+            }
+
+
+            //- Per face, per point (faceCentre) or triangle (faceDiag)
+            //  the added tet on the owner side
+            const labelListList& faceOwnerCells() const
+            {
+                return faceOwnerCells_;
+            }
+
+            //- Per face, per point (faceCentre) or triangle (faceDiag)
+            //  the added tet on the neighbour side
+            const labelListList& faceNeighbourCells() const
+            {
+                return faceNeighbourCells_;
+            }
+
+
+        // Edit
+
+            //- Insert all changes into meshMod to convert the polyMesh into
+            //  tets.
+            void setRefinement
+            (
+                const decompositionType decomposeType,
+                polyTopoChange& meshMod
+            );
+
+            //- Force recalculation of locally stored data on topological change
+            void updateMesh(const mapPolyMesh&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 3ba9c53475d21c85503dd72c16c5772a08466faf Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 27 Apr 2012 11:39:24 +0100
Subject: [PATCH 3/5] ENH: motionSmoother: stabilise edge weight calculation

---
 src/dynamicMesh/motionSmoother/motionSmoother.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C
index cf686b7dd8..fa9c2f1baa 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.C
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.C
@@ -283,7 +283,7 @@ Foam::tmp<Foam::scalarField> Foam::motionSmoother::calcEdgeWeights
 
     forAll(edges, edgeI)
     {
-        wght[edgeI] = min(GREAT, 1.0/edges[edgeI].mag(points));
+        wght[edgeI] = 1.0/(edges[edgeI].mag(points)+SMALL);
     }
     return twght;
 }

From cba9a30d3d88ea10c24e507a1a28c323898f1fe5 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 27 Apr 2012 11:41:36 +0100
Subject: [PATCH 4/5] ENH: autoSnapDriverFeature: make bounding box 3D

---
 .../autoHexMeshDriver/autoSnapDriverFeature.C   | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 2ebf0b140d..6347d5dd16 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1202,10 +1202,23 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     // Reverse lookup : go through all edgeAttractors and find the
     // nearest point on pp
 
+    // Get search domain and extend it a bit
+    treeBoundBox bb(pp.localPoints());
+    {
+        // Random number generator. Bit dodgy since not exactly random ;-)
+        Random rndGen(65431);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        bb = bb.extend(rndGen, 1e-4);
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    }
+
     indexedOctree<treeDataPoint> ppTree
     (
         treeDataPoint(pp.localPoints()),
-        treeBoundBox(pp.localPoints()), // overall search domain
+        bb,                             // overall search domain
         8,                              // maxLevel
         10,                             // leafsize
         3.0                             // duplicity

From 156839b16a85997f69f31066c9a16973bee98799 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 27 Apr 2012 12:18:54 +0100
Subject: [PATCH 5/5] BUG: tetDecomposer: updated tet numbering

---
 .../polyTopoChange/polyTopoChange/tetDecomposer.C         | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
index e8948c4fad..bdbd851527 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
@@ -696,16 +696,16 @@ void Foam::tetDecomposer::setRefinement
 
 void Foam::tetDecomposer::updateMesh(const mapPolyMesh& map)
 {
-    inplaceRenumber(map.pointMap(), cellToPoint_);
-    inplaceRenumber(map.pointMap(), faceToPoint_);
+    inplaceRenumber(map.reversePointMap(), cellToPoint_);
+    inplaceRenumber(map.reversePointMap(), faceToPoint_);
 
     forAll(faceOwnerCells_, faceI)
     {
-        inplaceRenumber(map.cellMap(), faceOwnerCells_[faceI]);
+        inplaceRenumber(map.reverseCellMap(), faceOwnerCells_[faceI]);
     }
     forAll(faceNeighbourCells_, faceI)
     {
-        inplaceRenumber(map.cellMap(), faceNeighbourCells_[faceI]);
+        inplaceRenumber(map.reverseCellMap(), faceNeighbourCells_[faceI]);
     }
 }