ENH: tetrahedron: move slicing with plane. Removed tetPointRef.H

applications/test/momentOfInertia/Test-momentOfInertia.C

@@ -35,7 +35,7 @@ Description
 #include "polyMesh.H"
 #include "ListOps.H"
 #include "face.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "triFaceList.H"
 #include "OFstream.H"
 #include "meshTools.H"

applications/utilities/mesh/manipulation/createBaffles/createBaffles.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H

@@ -43,7 +43,7 @@ SourceFiles
 #include "triFace.H"
 #include "edge.H"
 #include "pointField.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -42,7 +42,7 @@ SourceFiles
 #include "polyMesh.H"
 #include "coupledPolyPatch.H"
 #include "syncTools.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "tetIndices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -38,7 +38,7 @@ SourceFiles
 #define tetIndices_H
 
 #include "label.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "triPointRef.H"
 #include "polyMesh.H"
 #include "triFace.H"

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

@@ -25,7 +25,7 @@ License
 
 #include "primitiveMesh.H"
 #include "pyramidPointFaceRef.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "ListOps.H"
 #include "unitConversion.H"
 #include "SortableList.H"

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

@@ -1,52 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 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/>.
-
-Typedef
-    Foam::tetPointRef
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef tetPointRef_H
-#define tetPointRef_H
-
-#include "point.H"
-#include "tetrahedron.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-typedef tetrahedron<point, const point&> tetPointRef;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-#endif
-
-// ************************************************************************* //

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

@@ -35,9 +35,9 @@ SourceFiles
 #ifndef tetPoints_H
 #define tetPoints_H
 
+#include "tetrahedron.H"
 #include "FixedList.H"
 #include "treeBoundBox.H"
-#include "tetPointRef.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -54,6 +54,8 @@ namespace Foam
 
 class Istream;
 class Ostream;
+class tetPoints;
+class plane;
 
 // Forward declaration of friend functions and operators
 
@@ -73,6 +75,7 @@ inline Ostream& operator<<
     const tetrahedron<Point, PointRef>&
 );
 
+typedef tetrahedron<point, const point&> tetPointRef;
 
 /*---------------------------------------------------------------------------*\
                            class tetrahedron Declaration
@@ -81,10 +84,71 @@ inline Ostream& operator<<
 template<class Point, class PointRef>
 class tetrahedron
 {
+public:
+
+        // 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
+            {
+                FixedList<tetPoints, 200>& tets_;
+                label& nTets_;
+
+            public:
+                inline storeOp(FixedList<tetPoints, 200>&, 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:
 
@@ -184,6 +248,16 @@ 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;
+
+
             //- 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,7 +25,8 @@ License
 
 #include "triangle.H"
 #include "IOstreams.H"
-#include "triPointRef.H"
+#include "tetPoints.H"
+#include "plane.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -492,6 +493,483 @@ 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
+(
+    FixedList<tetPoints, 200>& 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;
+            p[0] = tet[i1];
+            p[1] = tet[i3];
+            p[2] = tet[i2];
+            p[3] = p01;
+            p[4] = p03;
+            p[5] = 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;
+            p[0] = tet[i3];
+            p[1] = tet[i1];
+            p[2] = tet[i2];
+            p[3] = p03;
+            p[4] = p01;
+            p[5] = 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;
+        label neg0 = -1;
+        label neg1 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                if (pos0 == -1)
+                {
+                    pos0 = i;
+                }
+                else
+                {
+                    pos1 = i;
+                }
+            }
+            else
+            {
+                if (neg0 == -1)
+                {
+                    neg0 = i;
+                }
+                else
+                {
+                    neg1 = 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;
+                p[0] = tet[0];
+                p[1] = p02;
+                p[2] = p03;
+                p[3] = tet[1];
+                p[4] = p12;
+                p[5] = p13;
+                //Pout<< "    01 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p02;
+                p[2] = p12;
+                p[3] = tet[3];
+                p[4] = p03;
+                p[5] = 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;
+                p[0] = tet[1];
+                p[1] = p01;
+                p[2] = p13;
+                p[3] = tet[2];
+                p[4] = p02;
+                p[5] = p23;
+                //Pout<< "    12 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[3];
+                p[1] = p23;
+                p[2] = p13;
+                p[3] = tet[0];
+                p[4] = p02;
+                p[5] = 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;
+                p[0] = tet[2];
+                p[1] = p12;
+                p[2] = p23;
+                p[3] = tet[0];
+                p[4] = p01;
+                p[5] = p03;
+                //Pout<< "    20 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p12;
+                p[2] = p01;
+                p[3] = tet[3];
+                p[4] = p23;
+                p[5] = 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;
+                p[0] = tet[3];
+                p[1] = p23;
+                p[2] = p13;
+                p[3] = tet[0];
+                p[4] = p02;
+                p[5] = p01;
+                //Pout<< "    03 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p23;
+                p[2] = p02;
+                p[3] = tet[1];
+                p[4] = p13;
+                p[5] = 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;
+                p[0] = tet[1];
+                p[1] = p12;
+                p[2] = p01;
+                p[3] = tet[3];
+                p[4] = p23;
+                p[5] = p03;
+                //Pout<< "    13 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p12;
+                p[2] = p23;
+                p[3] = tet[0];
+                p[4] = p01;
+                p[5] = 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;
+                p[0] = tet[2];
+                p[1] = p02;
+                p[2] = p12;
+                p[3] = tet[3];
+                p[4] = p03;
+                p[5] = p13;
+                //Pout<< "    23 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[0];
+                p[1] = p02;
+                p[2] = p03;
+                p[3] = tet[1];
+                p[4] = p12;
+                p[5] = p13;
+                //Pout<< "    23 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else
+        {
+            FatalErrorIn("tetSliceWithPlane(..)")
+                << "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;
+            p[0] = tet[i1];
+            p[1] = tet[i3];
+            p[2] = tet[i2];
+            p[3] = p01;
+            p[4] = p03;
+            p[5] = 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;
+            p[0] = tet[i3];
+            p[1] = tet[i1];
+            p[2] = tet[i2];
+            p[3] = p03;
+            p[4] = p01;
+            p[5] = 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/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H

@@ -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) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -26,7 +26,7 @@ License
 #include "polyMeshGeometry.H"
 #include "polyMeshTetDecomposition.H"
 #include "pyramidPointFaceRef.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "syncTools.H"
 #include "unitConversion.H"
 

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

@@ -25,7 +25,6 @@ License
 
 #include "cellPointWeight.H"
 #include "polyMesh.H"
-#include "tetPointRef.H"
 #include "polyMeshTetDecomposition.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

src/lagrangian/basic/Cloud/Cloud.H

@@ -41,9 +41,6 @@ SourceFiles
 #include "IOField.H"
 #include "CompactIOField.H"
 #include "polyMesh.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
-#include "tetPointRef.H"
 #include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/basic/particle/particle.H

@@ -35,11 +35,10 @@ Description
 #include "vector.H"
 #include "Cloud.H"
 #include "IDLList.H"
-#include "labelList.H"
 #include "pointField.H"
 #include "faceList.H"
 #include "OFstream.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "FixedList.H"
 #include "polyMeshTetDecomposition.H"
 

src/meshTools/momentOfInertia/momentOfInertia.H

@@ -37,7 +37,6 @@ SourceFiles
 #ifndef momentOfInertia_H
 #define momentOfInertia_H
 
-#include "tetPointRef.H"
 #include "triFaceList.H"
 #include "triSurface.H"
 #include "polyMesh.H"

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C

@@ -24,436 +24,27 @@ 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"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::tetOverlapVolume, 0);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::tetOverlapVolume::tetOverlapVolume()
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::tetOverlapVolume::~tetOverlapVolume()
-{}
-
-
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-Foam::point Foam::tetOverlapVolume::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 TetOp>
-inline void Foam::tetOverlapVolume::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 AboveTetOp, class BelowTetOp>
-inline void Foam::tetOverlapVolume::tetSliceWithPlane
-(
-    const tetPoints& tet,
-    const plane& pl,
-
-    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;
-            p[0] = tet[i1];
-            p[1] = tet[i3];
-            p[2] = tet[i2];
-            p[3] = p01;
-            p[4] = p03;
-            p[5] = 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;
-            p[0] = tet[i3];
-            p[1] = tet[i1];
-            p[2] = tet[i2];
-            p[3] = p03;
-            p[4] = p01;
-            p[5] = 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;
-        label neg0 = -1;
-        label neg1 = -1;
-        forAll(d, i)
-        {
-            if (d[i] > 0)
-            {
-                if (pos0 == -1)
-                {
-                    pos0 = i;
-                }
-                else
-                {
-                    pos1 = i;
-                }
-            }
-            else
-            {
-                if (neg0 == -1)
-                {
-                    neg0 = i;
-                }
-                else
-                {
-                    neg1 = 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;
-                p[0] = tet[0];
-                p[1] = p02;
-                p[2] = p03;
-                p[3] = tet[1];
-                p[4] = p12;
-                p[5] = p13;
-                //Pout<< "    01 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p02;
-                p[2] = p12;
-                p[3] = tet[3];
-                p[4] = p03;
-                p[5] = 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;
-                p[0] = tet[1];
-                p[1] = p01;
-                p[2] = p13;
-                p[3] = tet[2];
-                p[4] = p02;
-                p[5] = p23;
-                //Pout<< "    12 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[3];
-                p[1] = p23;
-                p[2] = p13;
-                p[3] = tet[0];
-                p[4] = p02;
-                p[5] = 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;
-                p[0] = tet[2];
-                p[1] = p12;
-                p[2] = p23;
-                p[3] = tet[0];
-                p[4] = p01;
-                p[5] = p03;
-                //Pout<< "    20 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[1];
-                p[1] = p12;
-                p[2] = p01;
-                p[3] = tet[3];
-                p[4] = p23;
-                p[5] = 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;
-                p[0] = tet[3];
-                p[1] = p23;
-                p[2] = p13;
-                p[3] = tet[0];
-                p[4] = p02;
-                p[5] = p01;
-                //Pout<< "    03 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p23;
-                p[2] = p02;
-                p[3] = tet[1];
-                p[4] = p13;
-                p[5] = 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;
-                p[0] = tet[1];
-                p[1] = p12;
-                p[2] = p01;
-                p[3] = tet[3];
-                p[4] = p23;
-                p[5] = p03;
-                //Pout<< "    13 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p12;
-                p[2] = p23;
-                p[3] = tet[0];
-                p[4] = p01;
-                p[5] = 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;
-                p[0] = tet[2];
-                p[1] = p02;
-                p[2] = p12;
-                p[3] = tet[3];
-                p[4] = p03;
-                p[5] = p13;
-                //Pout<< "    23 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[0];
-                p[1] = p02;
-                p[2] = p03;
-                p[3] = tet[1];
-                p[4] = p12;
-                p[5] = p13;
-                //Pout<< "    23 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else
-        {
-            FatalErrorIn("tetSliceWithPlane(..)") << "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;
-            p[0] = tet[i1];
-            p[1] = tet[i3];
-            p[2] = tet[i2];
-            p[3] = p01;
-            p[4] = p03;
-            p[5] = 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;
-            p[0] = tet[i3];
-            p[1] = tet[i1];
-            p[2] = tet[i2];
-            p[3] = p03;
-            p[4] = p01;
-            p[5] = p02;
-            //Pout<< "    belowprism:" << p << endl;
-            decomposePrism(p, belowOp);
-        }
-    }
-    else    // nPos == 0
-    {
-        belowOp(tet);
-    }
-}
-
-
 void Foam::tetOverlapVolume::tetTetOverlap
 (
     const tetPoints& tetA,
@@ -462,16 +53,15 @@ void Foam::tetOverlapVolume::tetTetOverlap
     label& nInside,
     FixedList<tetPoints, 200>& outsideTets,
     label& nOutside
-)
+) const
 {
     // Work storage
     FixedList<tetPoints, 200> cutInsideTets;
     label nCutInside = 0;
 
-    storeTetOp inside(insideTets, nInside);
+    tetPointRef::storeOp inside(insideTets, nInside);
-    storeTetOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
-    dummyTetOp outside;
+    tetPointRef::dummyOp outside;
-
 
 
     // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
@@ -484,7 +74,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
         // Cut and insert subtets into cutInsideTets (either by getting
         // an index from freeSlots or by appending to insideTets) or
         // insert into outsideTets
-        tetSliceWithPlane(tetA, pl0, cutInside, outside);
+        tetA.tet().sliceWithPlane(pl0, cutInside, outside);
     }
 
     if (nCutInside == 0)
@@ -501,7 +91,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
 
         for (label i = 0; i < nCutInside; i++)
         {
-            tetSliceWithPlane(cutInsideTets[i], pl1, inside, outside);
+            cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
         }
 
         if (nInside == 0)
@@ -518,7 +108,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
 
         for (label i = 0; i < nInside; i++)
         {
-            tetSliceWithPlane(insideTets[i], pl2, cutInside, outside);
+            insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
         }
 
         if (nCutInside == 0)
@@ -536,28 +126,27 @@ void Foam::tetOverlapVolume::tetTetOverlap
 
         for (label i = 0; i < nCutInside; i++)
         {
-            tetSliceWithPlane(cutInsideTets[i], pl3, inside, outside);
+            cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
         }
     }
 }
 
 
-inline Foam::scalar
+Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
-Foam::tetOverlapVolume::tetTetOverlapVol
 (
     const tetPoints& tetA,
     const tetPoints& tetB
-)
+) const
 {
     FixedList<tetPoints, 200> insideTets;
     label nInside = 0;
     FixedList<tetPoints, 200> cutInsideTets;
     label nCutInside = 0;
 
-    storeTetOp inside(insideTets, nInside);
+    tetPointRef::storeOp inside(insideTets, nInside);
-    storeTetOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
-    sumTetVolOp volInside;
+    tetPointRef::sumVolOp volInside;
-    dummyTetOp outside;
+    tetPointRef::dummyOp outside;
 
     // face0
     plane pl0(tetB[1], tetB[3], tetB[2]);
@@ -605,12 +194,12 @@ Foam::tetOverlapVolume::tetTetOverlapVol
 }
 
 
-inline const Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
+Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
 (
     const pointField& points,
     const face& f,
     const point& fc
-)
+) const
 {
     treeBoundBox bb(fc, fc);
     forAll(f, fp)
@@ -750,8 +339,8 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
 
 Foam::labelList Foam::tetOverlapVolume::overlappingCells
 (
-    const fvMesh& fromMesh,
+    const polyMesh& fromMesh,
-    const fvMesh& toMesh,
+    const polyMesh& toMesh,
     const label iTo
 ) const
 {
@@ -766,30 +355,4 @@ Foam::labelList Foam::tetOverlapVolume::overlappingCells
 }
 
 
-/*
-
-        forAll(cellsA, i)
-        {
-            label cellAI = cellsA[i];
-            treeBoundBox bbA
-            (
-                pointField(meshA.points(), meshA.cellPoints()[cellAI])
-            );
-
-            scalar v = cellCellOverlapVolumeMinDecomp
-            (
-                meshA,
-                cellAI,
-                bbA,
-                meshB,
-                cellBI,
-                bbB
-            );
-
-            overlapVol += v;
-            nOverlapTests++;
-        }
-
-*/
-
 // ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H

@@ -36,21 +36,17 @@ SourceFiles
 #ifndef tetOverlapVolume_H
 #define tetOverlapVolume_H
 
-#include "tetrahedron.H"
+#include "FixedList.H"
-#include "fvMesh.H"
+#include "labelList.H"
-#include "plane.H"
+#include "treeBoundBox.H"
-#include "tetPointRef.H"
-#include "OFstream.H"
-#include "meshTools.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
-#include "tetPoints.H"
-#include "tetCell.H"
-#include "EdgeMap.H"
 
 namespace Foam
 {
 
+class primitiveMesh;
+class polyMesh;
+class tetPoints;
+
 /*---------------------------------------------------------------------------*\
                          Class tetOverlapVolume Declaration
 \*---------------------------------------------------------------------------*/
@@ -59,82 +55,6 @@ class tetOverlapVolume
 {
     // Private member functions
 
-        //- Plane intersection
-        inline point planeIntersection
-        (
-            const FixedList<scalar, 4>& d,
-            const tetPoints& t,
-            const label negI,
-            const label posI
-        );
-
-
-        //- Decompose prism
-        template <class TetOp> inline void decomposePrism
-        (
-            const FixedList<point, 6>& points,
-            TetOp& op
-        );
-
-
-        //- Helping cľasses
-        class dummyTetOp
-        {
-         public:
-
-            inline void operator()(const tetPoints&){}
-        };
-
-
-        class sumTetVolOp
-        {
-         public:
-            scalar vol_;
-
-            inline sumTetVolOp()
-            :
-                vol_(0.0)
-            {}
-
-            inline void operator()(const tetPoints& tet)
-            {
-                vol_ += tet.tet().mag();
-            }
-        };
-
-
-        class storeTetOp
-        {
-            FixedList<tetPoints, 200>& tets_;
-            label& nTets_;
-
-         public:
-
-            inline storeTetOp(FixedList<tetPoints, 200>& tets, label& nTets)
-            :
-                tets_(tets),
-                nTets_(nTets)
-            {}
-
-            inline void operator()(const tetPoints& tet)
-            {
-                tets_[nTets_++] = tet;
-            }
-        };
-
-
-        //- Slice. Split tet into subtets above and below plane
-        template <class AboveTetOp, class BelowTetOp>
-        inline void tetSliceWithPlane
-        (
-            const tetPoints& tet,
-            const plane& pl,
-
-            AboveTetOp& aboveOp,
-            BelowTetOp& belowOp
-        );
-
-
         //- Tet overlap
         void tetTetOverlap
         (
@@ -144,31 +64,28 @@ class tetOverlapVolume
             label& nInside,
             FixedList<tetPoints, 200>& outsideTets,
             label& nOutside
-        );
+        ) const;
-
 
         //- Tet Overlap Vol
-        inline scalar tetTetOverlapVol
+        scalar tetTetOverlapVol
         (
             const tetPoints& tetA,
             const tetPoints& tetB
-        );
+        ) const;
-
 
         //- Return a const treeBoundBox
-        inline const treeBoundBox pyrBb
+        treeBoundBox pyrBb
         (
             const pointField& points,
             const face& f,
             const point& fc
-        );
+        ) const;
-
 
 public:
 
 
     //- Runtime type information
-    TypeName("tetOverlapVolume");
+    ClassName("tetOverlapVolume");
 
 
     // Constructors
@@ -177,18 +94,14 @@ public:
         tetOverlapVolume();
 
 
-    //- Destructor
-     virtual ~tetOverlapVolume();
-
-
     // Public members
 
         //- Return a list of cells in meshA which overlaps with cellBI in
         // meshB
         labelList overlappingCells
         (
-            const fvMesh& meshA,
+            const polyMesh& meshA,
-            const fvMesh& meshB,
+            const polyMesh& meshB,
             const label cellBI
         ) const;
 

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C

@@ -1,341 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Description
-    Calculation of shape function product for a tetrahedron
-
-\*---------------------------------------------------------------------------*/
-
-#include "tetrahedron.H"
-#include "triPointRef.H"
-#include "scalarField.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// (Probably very inefficient) minimum containment sphere calculation.
-// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
-// Sphere ctr is smallest one of
-// - tet circumcentre
-// - triangle circumcentre
-// - edge mids
-template<class Point, class PointRef>
-Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
-(
-    const scalar tol
-) const
-{
-    const scalar fac = 1 + tol;
-
-    // Halve order of tolerance for comparisons of sqr.
-    const scalar facSqr = Foam::sqrt(fac);
-
-
-    // 1. Circumcentre itself.
-
-    pointHit pHit(circumCentre());
-    pHit.setHit();
-    scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
-
-
-    // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
-    // check if 4th point is inside.
-
-    {
-        point ctr = triPointRef(a_, b_, c_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - a_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-    {
-        point ctr = triPointRef(a_, b_, d_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - a_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-    {
-        point ctr = triPointRef(a_, c_, d_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - a_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-    {
-        point ctr = triPointRef(b_, c_, d_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - b_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-
-    // 3. Try midpoints of edges
-
-    // mid of edge A-B
-    {
-        point ctr = 0.5*(a_ + b_);
-        scalar radiusSqr = magSqr(a_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(c_-ctr) <= testRadSrq
-         && magSqr(d_-ctr) <= testRadSrq)
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge A-C
-    {
-        point ctr = 0.5*(a_ + c_);
-        scalar radiusSqr = magSqr(a_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(b_-ctr) <= testRadSrq
-         && magSqr(d_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge A-D
-    {
-        point ctr = 0.5*(a_ + d_);
-        scalar radiusSqr = magSqr(a_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(b_-ctr) <= testRadSrq
-         && magSqr(c_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge B-C
-    {
-        point ctr = 0.5*(b_ + c_);
-        scalar radiusSqr = magSqr(b_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(a_-ctr) <= testRadSrq
-         && magSqr(d_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge B-D
-    {
-        point ctr = 0.5*(b_ + d_);
-        scalar radiusSqr = magSqr(b_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(a_-ctr) <= testRadSrq
-         && magSqr(c_-ctr) <= testRadSrq)
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge C-D
-    {
-        point ctr = 0.5*(c_ + d_);
-        scalar radiusSqr = magSqr(c_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(a_-ctr) <= testRadSrq
-         && magSqr(b_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-
-    pHit.setDistance(sqrt(minRadiusSqr));
-
-    return pHit;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiSquared
-(
-    scalarField& buffer
-) const
-{
-    // Change of sign between face area vector and gradient
-    // does not matter because of square
-
-    // Warning: Added a mag to produce positive coefficients even if
-    // the tetrahedron is twisted inside out.  This is pretty
-    // dangerous, but essential for mesh motion.
-    scalar magVol = Foam::mag(mag());
-
-    buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
-    buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
-    buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
-    buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiDotGradNj
-(
-    scalarField& buffer
-) const
-{
-    // Warning. Ordering of edges needs to be the same for a tetrahedron
-    // class, a tetrahedron cell shape model and a tetCell
-
-    // Warning: Added a mag to produce positive coefficients even if
-    // the tetrahedron is twisted inside out.  This is pretty
-    // dangerous, but essential for mesh motion.
-
-    // Double change of sign between face area vector and gradient
-
-    scalar magVol = Foam::mag(mag());
-    vector sa = Sa();
-    vector sb = Sb();
-    vector sc = Sc();
-    vector sd = Sd();
-
-    buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
-    buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
-    buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
-    buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
-    buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
-    buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiGradNi
-(
-    tensorField& buffer
-) const
-{
-    // Change of sign between face area vector and gradient
-    // does not matter because of square
-
-    scalar magVol = Foam::mag(mag());
-
-    buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
-    buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
-    buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
-    buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
-(
-    tensorField& buffer
-) const
-{
-    // Warning. Ordering of edges needs to be the same for a tetrahedron
-    // class, a tetrahedron cell shape model and a tetCell
-
-    // Double change of sign between face area vector and gradient
-
-    scalar magVol = Foam::mag(mag());
-    vector sa = Sa();
-    vector sb = Sb();
-    vector sc = Sc();
-    vector sd = Sd();
-
-    buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
-    buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
-    buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
-    buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
-    buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
-    buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
-}
-
-
-// ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H

@@ -1,315 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::tetrahedron
-
-Description
-    A tetrahedron primitive.
-
-    Ordering of edges needs to be the same for a tetrahedron
-    class, a tetrahedron cell shape model and a tetCell.
-
-SourceFiles
-    tetrahedronI.H
-    tetrahedron.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef tetrahedron_H
-#define tetrahedron_H
-
-#include "point.H"
-#include "primitiveFieldsFwd.H"
-#include "pointHit.H"
-#include "cachedRandom.H"
-#include "Random.H"
-#include "FixedList.H"
-#include "UList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class Istream;
-class Ostream;
-class tetPoints;
-class plane;
-
-// Forward declaration of friend functions and operators
-
-template<class Point, class PointRef> class tetrahedron;
-
-template<class Point, class PointRef>
-inline Istream& operator>>
-(
-    Istream&,
-    tetrahedron<Point, PointRef>&
-);
-
-template<class Point, class PointRef>
-inline Ostream& operator<<
-(
-    Ostream&,
-    const tetrahedron<Point, PointRef>&
-);
-
-
-/*---------------------------------------------------------------------------*\
-                           class tetrahedron Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Point, class PointRef>
-class tetrahedron
-{
-public:
-
-        // 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
-            {
-                FixedList<tetPoints, 200>& tets_;
-                label& nTets_;
-
-            public:
-                inline storeOp(FixedList<tetPoints, 200>&, 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:
-
-    // Member constants
-
-        enum
-        {
-            nVertices = 4,  // Number of vertices in tetrahedron
-            nEdges = 6      // Number of edges in tetrahedron
-        };
-
-
-    // Constructors
-
-        //- Construct from points
-        inline tetrahedron
-        (
-            const Point& a,
-            const Point& b,
-            const Point& c,
-            const Point& d
-        );
-
-        //- Construct from four points in the list of points
-        inline tetrahedron
-        (
-            const UList<Point>&,
-            const FixedList<label, 4>& indices
-        );
-
-        //- Construct from Istream
-        inline tetrahedron(Istream&);
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return vertices
-            inline const Point& a() const;
-
-            inline const Point& b() const;
-
-            inline const Point& c() const;
-
-            inline const Point& d() const;
-
-
-        // Properties
-
-            //- Return face normal
-            inline vector Sa() const;
-
-            inline vector Sb() const;
-
-            inline vector Sc() const;
-
-            inline vector Sd() const;
-
-            //- Return centre (centroid)
-            inline Point centre() const;
-
-            //- Return volume
-            inline scalar mag() const;
-
-            //- Return circum-centre
-            inline Point circumCentre() const;
-
-            //- Return circum-radius
-            inline scalar circumRadius() const;
-
-            //- Return quality: Ratio of tetrahedron and circum-sphere
-            //  volume, scaled so that a regular tetrahedron has a
-            //  quality of 1
-            inline scalar quality() const;
-
-            //- Return a random point in the tetrahedron from a
-            //  uniform distribution
-            inline Point randomPoint(Random& rndGen) const;
-
-            //- Return a random point in the tetrahedron from a
-            //  uniform distribution
-            inline Point randomPoint(cachedRandom& rndGen) const;
-
-            //- Calculate the barycentric coordinates of the given
-            //  point, in the same order as a, b, c, d.  Returns the
-            //  determinant of the solution.
-            inline scalar barycentric
-            (
-                const point& pt,
-                List<scalar>& bary
-            ) const;
-
-            //- Return nearest point to p on tetrahedron
-            inline pointHit nearestPoint(const point& p) const;
-
-            //- Return 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;
-
-
-            //- Return (min)containment sphere, i.e. the smallest sphere with
-            //  all points inside. Returns pointHit with:
-            //  - hit         : if sphere is equal to circumsphere
-            //                  (biggest sphere)
-            //  - point       : centre of sphere
-            //  - distance    : radius of sphere
-            //  - eligiblemiss: false
-            // Tol (small compared to 1, e.g. 1E-9) is used to determine
-            // whether point is inside: mag(pt - ctr) < (1+tol)*radius.
-            pointHit containmentSphere(const scalar tol) const;
-
-            //- Fill buffer with shape function products
-            void gradNiSquared(scalarField& buffer) const;
-
-            void gradNiDotGradNj(scalarField& buffer) const;
-
-            void gradNiGradNi(tensorField& buffer) const;
-
-            void gradNiGradNj(tensorField& buffer) const;
-
-
-    // IOstream operators
-
-        friend Istream& operator>> <Point, PointRef>
-        (
-            Istream&,
-            tetrahedron&
-        );
-
-        friend Ostream& operator<< <Point, PointRef>
-        (
-            Ostream&,
-            const tetrahedron&
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "tetrahedronI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "tetrahedron.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //