diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C index 6a16d0df73..9180eb0986 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C +++ b/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 . -Description - Calculation of shape function product for a tetrahedron - \*---------------------------------------------------------------------------*/ #include "tetrahedron.H" @@ -32,6 +29,96 @@ Description // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +void Foam::tetrahedron::tetOverlap +( + const tetrahedron& 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 diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index 5ff25795b2..803cc7599e 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H +++ b/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 tetIntersectionList; + + // Classes for use in sliceWithPlane. What to do with decomposition // of tet. @@ -111,11 +118,11 @@ public: //- Store resulting tets class storeOp { - FixedList& tets_; + tetIntersectionList& tets_; label& nTets_; public: - inline storeOp(FixedList&, 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& 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: diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 5a6787be83..7831e2614b 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -556,7 +556,7 @@ inline void Foam::tetrahedron::sumVolOp::operator() template inline Foam::tetrahedron::storeOp::storeOp ( - FixedList& tets, + tetIntersectionList& tets, label& nTets ) : diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C index 16d311815c..a16a386f8a 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C +++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C @@ -49,14 +49,14 @@ void Foam::tetOverlapVolume::tetTetOverlap ( const tetPoints& tetA, const tetPoints& tetB, - FixedList& insideTets, + tetIntersectionList& insideTets, label& nInside, - FixedList& outsideTets, + tetIntersectionList& outsideTets, label& nOutside ) const { // Work storage - FixedList 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 insideTets; + tetIntersectionList insideTets; label nInside = 0; - FixedList 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]; diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H index 388f645f3f..c898f427d1 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H +++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H @@ -60,9 +60,9 @@ class tetOverlapVolume ( const tetPoints& tetA, const tetPoints& tetB, - FixedList& insideTets, + tetIntersectionList& insideTets, label& nInside, - FixedList& outsideTets, + tetIntersectionList& outsideTets, label& nOutside ) const; @@ -115,7 +115,7 @@ public: const primitiveMesh& meshB, const label cellBI, const treeBoundBox& cellBbB - ); + ) const; };