ENH: Reinstated local tet-based intersection

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

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::tetPoints
+
+Description
+    Tet storage. Null constructable (unfortunately tetrahedron<point, point>
+    is not)
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetPoints_H
+#define tetPoints_H
+
+#include "tetrahedron.H"
+#include "FixedList.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class tetPoints Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetPoints
+:
+    public FixedList<point, 4>
+{
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline tetPoints()
+        {}
+
+        //- Construct from four points
+        inline tetPoints
+        (
+            const point& a,
+            const point& b,
+            const point& c,
+            const point& d
+        )
+        {
+            operator[](0) = a;
+            operator[](1) = b;
+            operator[](2) = c;
+            operator[](3) = d;
+        }
+
+    // Member Functions
+
+        //- Return the tetrahedron
+        inline tetPointRef tet() const
+        {
+            return tetPointRef
+            (
+                operator[](0),
+                operator[](1),
+                operator[](2),
+                operator[](3)
+            );
+        }
+
+        //- Calculate the bounding box
+        inline treeBoundBox bounds() const
+        {
+            treeBoundBox bb(operator[](0));
+            for (label i = 1; i < size(); ++i)
+            {
+                bb.add(operator[](i));
+            }
+            return bb;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,8 @@ namespace Foam
 
 class Istream;
 class Ostream;
+class tetPoints;
+class plane;
 
 // Forward declaration of friend functions and operators
 
@@ -75,6 +77,8 @@ inline Ostream& operator<<
     const tetrahedron<Point, PointRef>&
 );
 
+typedef tetrahedron<point, const point&> tetPointRef;
+
 /*---------------------------------------------------------------------------*\
                            class tetrahedron Declaration
 \*---------------------------------------------------------------------------*/
@@ -82,12 +86,78 @@ inline Ostream& operator<<
 template<class Point, class PointRef>
 class tetrahedron
 {
+public:
+
+    // Public typedefs
+
+        //- Storage type for tets originating from intersecting tets.
+        //  (can possibly be smaller than 200)
+        typedef FixedList<tetPoints, 200> tetIntersectionList;
+
+
+        // Classes for use in sliceWithPlane. What to do with decomposition
+        // of tet.
+
+            //- Dummy
+            class dummyOp
+            {
+            public:
+                inline void operator()(const tetPoints&);
+            };
+
+            //- Sum resulting volumes
+            class sumVolOp
+            {
+            public:
+                scalar vol_;
+
+                inline sumVolOp();
+
+                inline void operator()(const tetPoints&);
+            };
+
+            //- Store resulting tets
+            class storeOp
+            {
+                tetIntersectionList& tets_;
+                label& nTets_;
+
+            public:
+                inline storeOp(tetIntersectionList&, label&);
+
+                inline void operator()(const tetPoints&);
+            };
+
 private:
 
     // Private data
 
         PointRef a_, b_, c_, d_;
 
+        inline static point planeIntersection
+        (
+            const FixedList<scalar, 4>&,
+            const tetPoints&,
+            const label,
+            const label
+        );
+
+        template<class TetOp>
+        inline static void decomposePrism
+        (
+            const FixedList<point, 6>& points,
+            TetOp& op
+        );
+
+        template<class AboveTetOp, class BelowTetOp>
+        inline static void tetSliceWithPlane
+        (
+            const plane& pl,
+            const tetPoints& tet,
+            AboveTetOp& aboveOp,
+            BelowTetOp& belowOp
+        );
+
 
 public:
 
@@ -170,10 +240,6 @@ public:
             //  uniform distribution
             inline Point randomPoint(Random& rndGen) const;
 
-            //- Return a random point in the tetrahedron from a
-            //  uniform distribution
-            inline Point randomPoint(cachedRandom& rndGen) const;
-
             //- Calculate the point from the given barycentric coordinates.
             inline Point barycentricToPoint(const barycentric& bary) const;
 
@@ -195,6 +261,26 @@ public:
             //- Return true if point is inside tetrahedron
             inline bool inside(const point& pt) const;
 
+            //- Decompose tet into tets above and below plane
+            template<class AboveTetOp, class BelowTetOp>
+            inline void sliceWithPlane
+            (
+                const plane& pl,
+                AboveTetOp& aboveOp,
+                BelowTetOp& belowOp
+            ) const;
+
+            //- Decompose tet into tets inside and outside other tet
+            inline void tetOverlap
+            (
+                const tetrahedron<Point, PointRef>& tetB,
+                tetIntersectionList& insideTets,
+                label& nInside,
+                tetIntersectionList& outsideTets,
+                label& nOutside
+            ) 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

@@ -25,6 +25,7 @@ License
 
 #include "triangle.H"
 #include "IOstreams.H"
+#include "tetPoints.H"
 #include "plane.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -485,6 +486,548 @@ bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
 }
 
 
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::dummyOp::operator()
+(
+    const tetPoints&
+)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::sumVolOp::sumVolOp()
+:
+    vol_(0.0)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
+(
+    const tetPoints& tet
+)
+{
+    vol_ += tet.tet().mag();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
+(
+    tetIntersectionList& tets,
+    label& nTets
+)
+:
+    tets_(tets),
+    nTets_(nTets)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::storeOp::operator()
+(
+    const tetPoints& tet
+)
+{
+    tets_[nTets_++] = tet;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::point Foam::tetrahedron<Point, PointRef>::planeIntersection
+(
+    const FixedList<scalar, 4>& d,
+    const tetPoints& t,
+    const label negI,
+    const label posI
+)
+{
+    return
+        (d[posI]*t[negI] - d[negI]*t[posI])
+      / (-d[negI]+d[posI]);
+}
+
+
+template<class Point, class PointRef>
+template<class TetOp>
+inline void Foam::tetrahedron<Point, PointRef>::decomposePrism
+(
+    const FixedList<point, 6>& points,
+    TetOp& op
+)
+{
+    op(tetPoints(points[1], points[3], points[2], points[0]));
+    op(tetPoints(points[1], points[2], points[3], points[4]));
+    op(tetPoints(points[4], points[2], points[3], points[5]));
+}
+
+
+template<class Point, class PointRef>
+template<class AboveTetOp, class BelowTetOp>
+inline void Foam::tetrahedron<Point, PointRef>::tetSliceWithPlane
+(
+    const plane& pl,
+    const tetPoints& tet,
+    AboveTetOp& aboveOp,
+    BelowTetOp& belowOp
+)
+{
+    // Distance to plane
+    FixedList<scalar, 4> d;
+    label nPos = 0;
+    forAll(tet, i)
+    {
+        d[i] = ((tet[i] - pl.refPoint()) & pl.normal());
+        if (d[i] > 0)
+        {
+            nPos++;
+        }
+    }
+
+    if (nPos == 4)
+    {
+        aboveOp(tet);
+    }
+    else if (nPos == 3)
+    {
+        // Sliced into below tet and above prism. Prism gets split into
+        // two tets.
+
+        // Find the below tet
+        label i0 = -1;
+        forAll(d, i)
+        {
+            if (d[i] <= 0)
+            {
+                i0 = i;
+                break;
+            }
+        }
+
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+        label i3 = d.fcIndex(i2);
+
+        point p01(planeIntersection(d, tet, i0, i1));
+        point p02(planeIntersection(d, tet, i0, i2));
+        point p03(planeIntersection(d, tet, i0, i3));
+
+        // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
+        //          ,,         1 :     ,,     inwards pointing triad
+        //          ,,         2 :     ,,     outwards pointing triad
+        //          ,,         3 :     ,,     inwards pointing triad
+
+        //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
+
+        if (i0 == 0 || i0 == 2)
+        {
+            tetPoints t(tet[i0], p01, p02, p03);
+            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 3, belowTet i0==0 or 2");
+            belowOp(t);
+
+            // Prism
+            FixedList<point, 6> p
+            (
+                {
+                    tet[i1],
+                    tet[i3],
+                    tet[i2],
+                    p01,
+                    p03,
+                    p02
+                }
+            );
+
+            //Pout<< "    aboveprism:" << p << endl;
+            decomposePrism(p, aboveOp);
+        }
+        else
+        {
+            tetPoints t(p01, p02, p03, tet[i0]);
+            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 3, belowTet i0==1 or 3");
+            belowOp(t);
+
+            // Prism
+            FixedList<point, 6> p
+            (
+                {
+                    tet[i3],
+                    tet[i1],
+                    tet[i2],
+                    p03,
+                    p01,
+                    p02
+                }
+            );
+            //Pout<< "    aboveprism:" << p << endl;
+            decomposePrism(p, aboveOp);
+        }
+    }
+    else if (nPos == 2)
+    {
+        // Tet cut into two prisms. Determine the positive one.
+        label pos0 = -1;
+        label pos1 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                if (pos0 == -1)
+                {
+                    pos0 = i;
+                }
+                else
+                {
+                    pos1 = i;
+                }
+            }
+        }
+
+        //Pout<< "Split 2pos tet " << tet << " d:" << d
+        //    << " around pos0:" << pos0 << " pos1:" << pos1
+        //    << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
+
+        const edge posEdge(pos0, pos1);
+
+        if (posEdge == edge(0, 1))
+        {
+            point p02(planeIntersection(d, tet, 0, 2));
+            point p03(planeIntersection(d, tet, 0, 3));
+            point p12(planeIntersection(d, tet, 1, 2));
+            point p13(planeIntersection(d, tet, 1, 3));
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[0],
+                        p02,
+                        p03,
+                        tet[1],
+                        p12,
+                        p13
+                    }
+                );
+
+                //Pout<< "    01 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[2],
+                        p02,
+                        p12,
+                        tet[3],
+                        p03,
+                        p13
+                    }
+                );
+
+                //Pout<< "    01 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(1, 2))
+        {
+            point p01(planeIntersection(d, tet, 0, 1));
+            point p13(planeIntersection(d, tet, 1, 3));
+            point p02(planeIntersection(d, tet, 0, 2));
+            point p23(planeIntersection(d, tet, 2, 3));
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[1],
+                        p01,
+                        p13,
+                        tet[2],
+                        p02,
+                        p23
+                    }
+                );
+
+                //Pout<< "    12 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[3],
+                        p23,
+                        p13,
+                        tet[0],
+                        p02,
+                        p01
+                    }
+                );
+
+                //Pout<< "    12 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(2, 0))
+        {
+            point p01(planeIntersection(d, tet, 0, 1));
+            point p03(planeIntersection(d, tet, 0, 3));
+            point p12(planeIntersection(d, tet, 1, 2));
+            point p23(planeIntersection(d, tet, 2, 3));
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[2],
+                        p12,
+                        p23,
+                        tet[0],
+                        p01,
+                        p03
+                    }
+                );
+
+                //Pout<< "    20 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[1],
+                        p12,
+                        p01,
+                        tet[3],
+                        p23,
+                        p03
+                    }
+                );
+
+                //Pout<< "    20 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(0, 3))
+        {
+            point p01(planeIntersection(d, tet, 0, 1));
+            point p02(planeIntersection(d, tet, 0, 2));
+            point p13(planeIntersection(d, tet, 1, 3));
+            point p23(planeIntersection(d, tet, 2, 3));
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[3],
+                        p23,
+                        p13,
+                        tet[0],
+                        p02,
+                        p01
+                    }
+                );
+
+                //Pout<< "    03 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[2],
+                        p23,
+                        p02,
+                        tet[1],
+                        p13,
+                        p01
+                    }
+                );
+
+                //Pout<< "    03 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(1, 3))
+        {
+            point p01(planeIntersection(d, tet, 0, 1));
+            point p12(planeIntersection(d, tet, 1, 2));
+            point p03(planeIntersection(d, tet, 0, 3));
+            point p23(planeIntersection(d, tet, 2, 3));
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[1],
+                        p12,
+                        p01,
+                        tet[3],
+                        p23,
+                        p03
+                    }
+                );
+
+                //Pout<< "    13 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[2],
+                        p12,
+                        p23,
+                        tet[0],
+                        p01,
+                        p03
+                    }
+                );
+
+                //Pout<< "    13 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(2, 3))
+        {
+            point p02(planeIntersection(d, tet, 0, 2));
+            point p12(planeIntersection(d, tet, 1, 2));
+            point p03(planeIntersection(d, tet, 0, 3));
+            point p13(planeIntersection(d, tet, 1, 3));
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[2],
+                        p02,
+                        p12,
+                        tet[3],
+                        p03,
+                        p13
+                    }
+                );
+
+                //Pout<< "    23 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p
+                (
+                    {
+                        tet[0],
+                        p02,
+                        p03,
+                        tet[1],
+                        p12,
+                        p13
+                    }
+                );
+
+                //Pout<< "    23 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Missed edge:" << posEdge
+                << abort(FatalError);
+        }
+    }
+    else if (nPos == 1)
+    {
+        // Find the positive tet
+        label i0 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                i0 = i;
+                break;
+            }
+        }
+
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+        label i3 = d.fcIndex(i2);
+
+        point p01(planeIntersection(d, tet, i0, i1));
+        point p02(planeIntersection(d, tet, i0, i2));
+        point p03(planeIntersection(d, tet, i0, i3));
+
+        //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
+
+        if (i0 == 0 || i0 == 2)
+        {
+            tetPoints t(tet[i0], p01, p02, p03);
+            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
+            aboveOp(t);
+
+            // Prism
+            FixedList<point, 6> p
+            (
+                {
+                    tet[i1],
+                    tet[i3],
+                    tet[i2],
+                    p01,
+                    p03,
+                    p02
+                }
+            );
+
+            //Pout<< "    belowprism:" << p << endl;
+            decomposePrism(p, belowOp);
+        }
+        else
+        {
+            tetPoints t(p01, p02, p03, tet[i0]);
+            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
+            aboveOp(t);
+
+            // Prism
+            FixedList<point, 6> p
+            (
+                {
+                    tet[i3],
+                    tet[i1],
+                    tet[i2],
+                    p03,
+                    p01,
+                    p02
+                }
+            );
+
+            //Pout<< "    belowprism:" << p << endl;
+            decomposePrism(p, belowOp);
+        }
+    }
+    else    // nPos == 0
+    {
+        belowOp(tet);
+    }
+}
+
+
+template<class Point, class PointRef>
+template<class AboveTetOp, class BelowTetOp>
+inline void Foam::tetrahedron<Point, PointRef>::sliceWithPlane
+(
+    const plane& pl,
+    AboveTetOp& aboveOp,
+    BelowTetOp& belowOp
+) const
+{
+    tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp);
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>

src/meshTools/tetOverlapVolume/tetOverlapVolume.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -24,12 +24,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tetOverlapVolume.H"
+#include "tetrahedron.H"
+#include "tetPoints.H"
 #include "polyMesh.H"
 #include "OFstream.H"
 #include "treeBoundBox.H"
 #include "indexedOctree.H"
 #include "treeDataCell.H"
-#include "cut.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -47,67 +48,6 @@ Foam::tetOverlapVolume::tetOverlapVolume()
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
-(
-    const tetPointRef& tetA,
-    const tetPointRef& tetB
-) const
-{
-    // A maximum of three cuts are made (the tets that result from the final cut
-    // are not stored), and each cut can create at most three tets. The
-    // temporary storage must therefore extend to 3^3 = 27 tets.
-    typedef cutTetList<27> tetListType;
-    static tetListType cutTetList1, cutTetList2;
-
-    // face 0
-    const plane pl0(tetB.b(), tetB.d(), tetB.c());
-    const FixedList<point, 4> t({tetA.a(), tetA.b(), tetA.c(), tetA.d()});
-    cutTetList1.clear();
-    tetCut(t, pl0, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
-    if (cutTetList1.size() == 0)
-    {
-        return 0;
-    }
-
-    // face 1
-    const plane pl1(tetB.a(), tetB.c(), tetB.d());
-    cutTetList2.clear();
-    for (label i = 0; i < cutTetList1.size(); i++)
-    {
-        const FixedList<point, 4>& t = cutTetList1[i];
-        tetCut(t, pl1, cut::appendOp<tetListType>(cutTetList2), cut::noOp());
-    }
-    if (cutTetList2.size() == 0)
-    {
-        return 0;
-    }
-
-    // face 2
-    const plane pl2(tetB.a(), tetB.d(), tetB.b());
-    cutTetList1.clear();
-    for (label i = 0; i < cutTetList2.size(); i++)
-    {
-        const FixedList<point, 4>& t = cutTetList2[i];
-        tetCut(t, pl2, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
-    }
-    if (cutTetList1.size() == 0)
-    {
-        return 0;
-    }
-
-    // face 3
-    const plane pl3(tetB.a(), tetB.b(), tetB.c());
-    scalar v = 0;
-    for (label i = 0; i < cutTetList1.size(); i++)
-    {
-        const FixedList<point, 4>& t = cutTetList1[i];
-        v += tetCut(t, pl3, cut::volumeOp(), cut::noOp());
-    }
-
-    return v;
-}
-
-
 Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
 (
     const pointField& points,
@@ -134,122 +74,43 @@ bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp
     const scalar threshold
 ) const
 {
-    const cell& cFacesA = meshA.cells()[cellAI];
+    hasOverlapOp overlapCheckOp(threshold);
-    const point& ccA = meshA.cellCentres()[cellAI];
+    cellCellOverlapMinDecomp<hasOverlapOp>
-
-    const cell& cFacesB = meshB.cells()[cellBI];
-    const point& ccB = meshB.cellCentres()[cellBI];
-
-    scalar vol = 0.0;
-
-    forAll(cFacesA, cFA)
-    {
-        label faceAI = cFacesA[cFA];
-
-        const face& fA = meshA.faces()[faceAI];
-        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
-        if (!pyrA.overlaps(cellBbB))
-        {
-            continue;
-        }
-
-        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
-
-        label tetBasePtAI = 0;
-
-        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
-
-        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
-        {
-            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
-            label otherFacePtAI = fA.fcIndex(facePtAI);
-
-            label pt0I = -1;
-            label pt1I = -1;
-
-            if (ownA)
-            {
-                pt0I = fA[facePtAI];
-                pt1I = fA[otherFacePtAI];
-            }
-            else
-            {
-                pt0I = fA[otherFacePtAI];
-                pt1I = fA[facePtAI];
-            }
-
-            const tetPointRef tetA
     (
-                ccA,
+        meshA,
-                tetBasePtA,
+        cellAI,
-                meshA.points()[pt0I],
+        meshB,
-                meshA.points()[pt1I]
+        cellBI,
-            );
+        cellBbB,
-            const treeBoundBox tetABb(tetA.bounds());
+        overlapCheckOp
-
-
-            // Loop over tets of cellB
-            forAll(cFacesB, cFB)
-            {
-                label faceBI = cFacesB[cFB];
-
-                const face& fB = meshB.faces()[faceBI];
-                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
-                if (!pyrB.overlaps(pyrA))
-                {
-                    continue;
-                }
-
-                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
-
-                label tetBasePtBI = 0;
-
-                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
-
-                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
-                {
-                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
-                    label otherFacePtBI = fB.fcIndex(facePtBI);
-
-                    label pt0I = -1;
-                    label pt1I = -1;
-
-                    if (ownB)
-                    {
-                        pt0I = fB[facePtBI];
-                        pt1I = fB[otherFacePtBI];
-                    }
-                    else
-                    {
-                        pt0I = fB[otherFacePtBI];
-                        pt1I = fB[facePtBI];
-                    }
-
-                    const tetPointRef tetB
-                    (
-                        ccB,
-                        tetBasePtB,
-                        meshB.points()[pt0I],
-                        meshB.points()[pt1I]
     );
 
-                    if (!tetB.bounds().overlaps(tetABb))
+    return overlapCheckOp.ok_;
+}
+
+
+Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
+
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB
+) const
 {
-                        continue;
+    sumOverlapOp overlapSumOp;
-                    }
+    cellCellOverlapMinDecomp<sumOverlapOp>
+    (
+        meshA,
+        cellAI,
+        meshB,
+        cellBI,
+        cellBbB,
+        overlapSumOp
+    );
 
-                    vol += tetTetOverlapVol(tetA, tetB);
+    return overlapSumOp.iop_.vol_;
-
-                    if (vol > threshold)
-                    {
-                        return true;
-                    }
-                }
-            }
-        }
-    }
-
-    return false;
 }
 
 
@@ -264,116 +125,18 @@ Foam::tetOverlapVolume::cellCellOverlapMomentMinDecomp
     const treeBoundBox& cellBbB
 ) const
 {
-    const cell& cFacesA = meshA.cells()[cellAI];
+    sumOverlapMomentOp overlapSumOp;
-    const point& ccA = meshA.cellCentres()[cellAI];
+    cellCellOverlapMinDecomp<sumOverlapMomentOp>
-
-    const cell& cFacesB = meshB.cells()[cellBI];
-    const point& ccB = meshB.cellCentres()[cellBI];
-
-    scalar vol = 0.0;
-
-    forAll(cFacesA, cFA)
-    {
-        label faceAI = cFacesA[cFA];
-
-        const face& fA = meshA.faces()[faceAI];
-        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
-        if (!pyrA.overlaps(cellBbB))
-        {
-            continue;
-        }
-
-        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
-
-        label tetBasePtAI = 0;
-
-        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
-
-        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
-        {
-            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
-            label otherFacePtAI = fA.fcIndex(facePtAI);
-
-            label pt0I = -1;
-            label pt1I = -1;
-
-            if (ownA)
-            {
-                pt0I = fA[facePtAI];
-                pt1I = fA[otherFacePtAI];
-            }
-            else
-            {
-                pt0I = fA[otherFacePtAI];
-                pt1I = fA[facePtAI];
-            }
-
-            const tetPointRef tetA
     (
-                ccA,
+        meshA,
-                tetBasePtA,
+        cellAI,
-                meshA.points()[pt0I],
+        meshB,
-                meshA.points()[pt1I]
+        cellBI,
+        cellBbB,
+        overlapSumOp
     );
-            const treeBoundBox tetABb(tetA.bounds());
 
-
+    return overlapSumOp.iop_.vol_;
-            // Loop over tets of cellB
-            forAll(cFacesB, cFB)
-            {
-                label faceBI = cFacesB[cFB];
-
-                const face& fB = meshB.faces()[faceBI];
-                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
-                if (!pyrB.overlaps(pyrA))
-                {
-                    continue;
-                }
-
-                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
-
-                label tetBasePtBI = 0;
-
-                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
-
-                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
-                {
-                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
-                    label otherFacePtBI = fB.fcIndex(facePtBI);
-
-                    label pt0I = -1;
-                    label pt1I = -1;
-
-                    if (ownB)
-                    {
-                        pt0I = fB[facePtBI];
-                        pt1I = fB[otherFacePtBI];
-                    }
-                    else
-                    {
-                        pt0I = fB[otherFacePtBI];
-                        pt1I = fB[facePtBI];
-                    }
-
-                    const tetPointRef tetB
-                    (
-                        ccB,
-                        tetBasePtB,
-                        meshB.points()[pt0I],
-                        meshB.points()[pt1I]
-                    );
-                    if (!tetB.bounds().overlaps(tetABb))
-                    {
-                        continue;
-                    }
-
-                    vol += tetTetOverlapVol(tetA, tetB);
-                }
-            }
-        }
-    }
-
-    return vol;
 }
 
 

src/meshTools/tetOverlapVolume/tetOverlapVolume.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ Description
 
 SourceFiles
     tetOverlapVolume.C
+    tetOverlapVolumeTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,7 +39,8 @@ SourceFiles
 #include "FixedList.H"
 #include "labelList.H"
 #include "treeBoundBox.H"
-#include "tetPointRef.H"
+#include "Tuple2.H"
+#include "tetrahedron.H"
 
 namespace Foam
 {
@@ -52,67 +54,124 @@ class polyMesh;
 
 class tetOverlapVolume
 {
+    // Private classes
+
+        //- tetPoints handling : sum resulting volumes
+        class sumMomentOp
+        {
+        public:
+            Tuple2<scalar, point> vol_;
+
+            inline sumMomentOp()
+            :
+                vol_(0.0, Zero)
+            {}
+
+            inline void operator()(const tetPoints& tet)
+            {
+                const tetPointRef t(tet.tet());
+                scalar tetVol = t.mag();
+                vol_.first() += tetVol;
+                vol_.second() += (tetVol*t.centre());
+            }
+        };
+
+        //- tetPoints combining : check for overlap
+        class hasOverlapOp
+        {
+        public:
+
+            const scalar threshold_;
+            tetPointRef::sumVolOp iop_;
+            bool ok_;
+
+            inline hasOverlapOp(const scalar threshold)
+            :
+                threshold_(threshold),
+                iop_(),
+                ok_(false)
+            {}
+
+            //- Overlap two tets
+            inline bool operator()(const tetPoints& A, const tetPoints& B)
+            {
+                tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
+                ok_ = (iop_.vol_ > threshold_);
+                return ok_;
+            }
+        };
+
+        //- tetPoints combining : sum overlap volume
+        class sumOverlapOp
+        {
+        public:
+
+            tetPointRef::sumVolOp iop_;
+
+            inline sumOverlapOp()
+            :
+                iop_()
+            {}
+
+            //- Overlap two tets
+            inline bool operator()(const tetPoints& A, const tetPoints& B)
+            {
+                tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
+                return false;
+            }
+        };
+
+        //- tetPoints combining : sum overlap volume
+        class sumOverlapMomentOp
+        {
+        public:
+
+            sumMomentOp iop_;
+
+            inline sumOverlapMomentOp()
+            :
+                iop_()
+            {}
+
+            //- Overlap two tets
+            inline bool operator()(const tetPoints& A, const tetPoints& B)
+            {
+                tetTetOverlap<sumMomentOp>(A, B, iop_);
+                return false;
+            }
+        };
+
+
     // Private member functions
 
-        //- Tet overlap volume
+        //- Tet overlap calculation
-        scalar tetTetOverlapVol
+        template<class tetPointsOp>
+        static void tetTetOverlap
         (
-            const tetPointRef& tetA,
+            const tetPoints& tetA,
-            const tetPointRef& tetB
+            const tetPoints& tetB,
-        ) const;
+            tetPointsOp& insideOp
+        );
+
+        //- Cell overlap calculation
+        template<class tetsOp>
+        static void cellCellOverlapMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB,
+            tetsOp& combineTetsOp
+        );
 
         //- Return a const treeBoundBox
-        treeBoundBox pyrBb
+        static treeBoundBox pyrBb
         (
             const pointField& points,
             const face& f,
             const point& fc
-        ) const;
+        );
-
-
-    // Private classes
-
-        //- A fixed list of tets which simulates a dynamic list by incrementing
-        //  a counter whenever its append method is called. This is used as an
-        //  optimisation for the tetTetOverlapVol method.
-        template<unsigned Size>
-        class cutTetList
-        :
-            public FixedList<FixedList<point, 4>, Size>
-        {
-        private:
-
-            //- The number of stored elements
-            label n_;
-
-
-        public:
-
-            //- Construct null
-            cutTetList()
-            :
-                n_(0)
-            {}
-
-            //- Clear the array
-            void clear()
-            {
-                n_ = 0;
-            }
-
-            //- Get the current size
-            label size() const
-            {
-                return n_;
-            }
-
-            //- Add a new tet to the end of the array
-            void append(const FixedList<point, 4>& t)
-            {
-                this->operator[](n_) = t;
-                ++ n_;
-            }
-        };
 
 
 public:
@@ -159,6 +218,17 @@ public:
             const label cellBI,
             const treeBoundBox& cellBbB
         ) const;
+
+        //- Calculates the overlap volume and moment
+        Tuple2<scalar, point> cellCellOverlapMomentMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB
+        ) const;
 };
 
 
@@ -168,6 +238,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "tetOverlapVolumeTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C

@@ -0,0 +1,259 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "tetOverlapVolume.H"
+#include "primitiveMesh.H"
+
+// * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+template<class tetPointsOp>
+void Foam::tetOverlapVolume::tetTetOverlap
+(
+    const tetPoints& tetA,
+    const tetPoints& tetB,
+    tetPointsOp& insideOp
+)
+{
+    static tetPointRef::tetIntersectionList insideTets;
+    label nInside = 0;
+    static tetPointRef::tetIntersectionList cutInsideTets;
+    label nCutInside = 0;
+
+    tetPointRef::storeOp inside(insideTets, nInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::dummyOp outside;
+
+    // Precompute the tet face areas and exit early if any face area is
+    // too small
+    static FixedList<vector, 4> tetAFaceAreas;
+    static FixedList<scalar, 4> tetAMag2FaceAreas;
+    tetPointRef tetATet = tetA.tet();
+    for (label facei = 0; facei < 4; ++facei)
+    {
+        tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
+        tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
+        if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
+    }
+
+    static FixedList<vector, 4> tetBFaceAreas;
+    static FixedList<scalar, 4> tetBMag2FaceAreas;
+    tetPointRef tetBTet = tetB.tet();
+    for (label facei = 0; facei < 4; ++facei)
+    {
+        tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
+        tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
+        if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
+    }
+
+
+    // Face 0
+    {
+        vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
+        plane pl0(tetBTet.tri(0).a(), n, false);
+
+        tetA.tet().sliceWithPlane(pl0, cutInside, outside);
+        if (nCutInside == 0)
+        {
+            return;
+        }
+    }
+
+    // Face 1
+    {
+        vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
+        plane pl1(tetBTet.tri(1).a(), n, false);
+
+        nInside = 0;
+        for (label i = 0; i < nCutInside; i++)
+        {
+            const tetPointRef t = cutInsideTets[i].tet();
+            t.sliceWithPlane(pl1, inside, outside);
+        }
+        if (nInside == 0)
+        {
+            return;
+        }
+    }
+
+    // Face 2
+    {
+        vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
+        plane pl2(tetBTet.tri(2).a(), n, false);
+
+        nCutInside = 0;
+        for (label i = 0; i < nInside; i++)
+        {
+            const tetPointRef t = insideTets[i].tet();
+            t.sliceWithPlane(pl2, cutInside, outside);
+        }
+        if (nCutInside == 0)
+        {
+            return;
+        }
+    }
+
+    // Face 3
+    {
+        vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
+        plane pl3(tetBTet.tri(3).a(), n, false);
+        for (label i = 0; i < nCutInside; i++)
+        {
+            const tetPointRef t = cutInsideTets[i].tet();
+            t.sliceWithPlane(pl3, insideOp, outside);
+        }
+    }
+}
+
+
+template<class tetsOp>
+void Foam::tetOverlapVolume::cellCellOverlapMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
+
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB,
+    tetsOp& combineTetsOp
+)
+{
+    const cell& cFacesA = meshA.cells()[cellAI];
+    const point& ccA = meshA.cellCentres()[cellAI];
+
+    const cell& cFacesB = meshB.cells()[cellBI];
+    const point& ccB = meshB.cellCentres()[cellBI];
+
+    forAll(cFacesA, cFA)
+    {
+        label faceAI = cFacesA[cFA];
+
+        const face& fA = meshA.faces()[faceAI];
+        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
+        if (!pyrA.overlaps(cellBbB))
+        {
+            continue;
+        }
+
+        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
+
+        label tetBasePtAI = 0;
+
+        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
+
+        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
+        {
+            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
+            label otherFacePtAI = fA.fcIndex(facePtAI);
+
+            label pt0I = -1;
+            label pt1I = -1;
+
+            if (ownA)
+            {
+                pt0I = fA[facePtAI];
+                pt1I = fA[otherFacePtAI];
+            }
+            else
+            {
+                pt0I = fA[otherFacePtAI];
+                pt1I = fA[facePtAI];
+            }
+
+            const tetPoints tetA
+            (
+                ccA,
+                tetBasePtA,
+                meshA.points()[pt0I],
+                meshA.points()[pt1I]
+            );
+            const treeBoundBox tetABb(tetA.bounds());
+
+            // Loop over tets of cellB
+            forAll(cFacesB, cFB)
+            {
+                label faceBI = cFacesB[cFB];
+
+                const face& fB = meshB.faces()[faceBI];
+                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
+                if (!pyrB.overlaps(pyrA))
+                {
+                    continue;
+                }
+
+                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
+
+                label tetBasePtBI = 0;
+
+                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
+
+                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
+                {
+                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
+                    label otherFacePtBI = fB.fcIndex(facePtBI);
+
+                    label pt0I = -1;
+                    label pt1I = -1;
+
+                    if (ownB)
+                    {
+                        pt0I = fB[facePtBI];
+                        pt1I = fB[otherFacePtBI];
+                    }
+                    else
+                    {
+                        pt0I = fB[otherFacePtBI];
+                        pt1I = fB[facePtBI];
+                    }
+
+                    const tetPoints tetB
+                    (
+                        ccB,
+                        tetBasePtB,
+                        meshB.points()[pt0I],
+                        meshB.points()[pt1I]
+                    );
+                    if (!tetB.bounds().overlaps(tetABb))
+                    {
+                        continue;
+                    }
+
+                    if (combineTetsOp(tetA, tetB))
+                    {
+                        return;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //