ENH: tetrahedron: add tet-tet overlap

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     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"
@@ -32,6 +29,96 @@ Description
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::tetOverlap
+(
+    const tetrahedron<Point, PointRef>& tetB,
+    tetIntersectionList& insideTets,
+    label& nInside,
+    tetIntersectionList& outsideTets,
+    label& nOutside
+) const
+{
+    // Work storage
+    tetIntersectionList cutInsideTets;
+    label nCutInside = 0;
+
+    nInside = 0;
+    storeOp inside(insideTets, nInside);
+    storeOp cutInside(cutInsideTets, nCutInside);
+
+    nOutside = 0;
+    storeOp outside(outsideTets, nOutside);
+
+
+    // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
+    // in aboveTets are inside tetB.
+
+    {
+        // face0
+        plane pl0(tetB.b_, tetB.d_, tetB.c_);
+
+        // Cut and insert subtets into cutInsideTets (either by getting
+        // an index from freeSlots or by appending to insideTets) or
+        // insert into outsideTets
+        sliceWithPlane(pl0, cutInside, outside);
+    }
+
+    if (nCutInside == 0)
+    {
+        nInside = nCutInside;
+        return;
+    }
+
+    {
+        // face1
+        plane pl1(tetB.a_, tetB.c_, tetB.d_);
+
+        nInside = 0;
+
+        for (label i = 0; i < nCutInside; i++)
+        {
+            cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
+        }
+
+        if (nInside == 0)
+        {
+            return;
+        }
+    }
+
+    {
+        // face2
+        plane pl2(tetB.a_, tetB.d_, tetB.b_);
+
+        nCutInside = 0;
+
+        for (label i = 0; i < nInside; i++)
+        {
+            insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
+        }
+
+        if (nCutInside == 0)
+        {
+            nInside = nCutInside;
+            return;
+        }
+    }
+
+    {
+        // face3
+        plane pl3(tetB.a_, tetB.b_, tetB.c_);
+
+        nInside = 0;
+
+        for (label i = 0; i < nCutInside; i++)
+        {
+            cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
+        }
+    }
+}
+
+
 // (Probably very inefficient) minimum containment sphere calculation.
 // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
 // Sphere ctr is smallest one of

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

@@ -87,6 +87,13 @@ 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.
 
@@ -111,11 +118,11 @@ public:
             //- Store resulting tets
             class storeOp
             {
-                FixedList<tetPoints, 200>& tets_;
+                tetIntersectionList& tets_;
                 label& nTets_;
 
             public:
-                inline storeOp(FixedList<tetPoints, 200>&, label&);
+                inline storeOp(tetIntersectionList&, label&);
 
                 inline void operator()(const tetPoints&);
             };
@@ -261,6 +268,16 @@ public:
                 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:

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

@@ -556,7 +556,7 @@ inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
 template<class Point, class PointRef>
 inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
 (
-    FixedList<tetPoints, 200>& tets,
+    tetIntersectionList& tets,
     label& nTets
 )
 :

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C

@@ -49,14 +49,14 @@ void Foam::tetOverlapVolume::tetTetOverlap
 (
     const tetPoints& tetA,
     const tetPoints& tetB,
-    FixedList<tetPoints, 200>& insideTets,
+    tetIntersectionList& insideTets,
     label& nInside,
-    FixedList<tetPoints, 200>& outsideTets,
+    tetIntersectionList& outsideTets,
     label& nOutside
 ) const
 {
     // Work storage
-    FixedList<tetPoints, 200> cutInsideTets;
+    tetIntersectionList cutInsideTets;
     label nCutInside = 0;
 
     tetPointRef::storeOp inside(insideTets, nInside);
@@ -138,9 +138,9 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
     const tetPoints& tetB
 ) const
 {
-    FixedList<tetPoints, 200> insideTets;
+    tetIntersectionList insideTets;
     label nInside = 0;
-    FixedList<tetPoints, 200> cutInsideTets;
+    tetIntersectionList cutInsideTets;
     label nCutInside = 0;
 
     tetPointRef::storeOp inside(insideTets, nInside);
@@ -222,7 +222,7 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
     const primitiveMesh& meshB,
     const label cellBI,
     const treeBoundBox& cellBbB
-)
+) const
 {
     const cell& cFacesA = meshA.cells()[cellAI];
     const point& ccA = meshA.cellCentres()[cellAI];

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H

@@ -60,9 +60,9 @@ class tetOverlapVolume
         (
             const tetPoints& tetA,
             const tetPoints& tetB,
-            FixedList<tetPoints, 200>& insideTets,
+            tetIntersectionList& insideTets,
             label& nInside,
-            FixedList<tetPoints, 200>& outsideTets,
+            tetIntersectionList& outsideTets,
             label& nOutside
         ) const;
 
@@ -115,7 +115,7 @@ public:
             const primitiveMesh& meshB,
             const label cellBI,
             const treeBoundBox& cellBbB
-        );
+        ) const;
 };