From 2b9b4c1aee532d8d7487d0328e286e432b267e1a Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 5 Sep 2017 09:32:21 +0100
Subject: [PATCH] ENH: Reinstated local tet-based intersection
---
.../primitiveShapes/tetrahedron/tetPoints.H | 113 ++++
.../primitiveShapes/tetrahedron/tetrahedron.H | 96 +++-
.../tetrahedron/tetrahedronI.H | 543 ++++++++++++++++++
.../tetOverlapVolume/tetOverlapVolume.C | 341 ++---------
.../tetOverlapVolume/tetOverlapVolume.H | 162 ++++--
.../tetOverlapVolumeTemplates.C | 259 +++++++++
6 files changed, 1177 insertions(+), 337 deletions(-)
create mode 100644 src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
create mode 100644 src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
new file mode 100644
index 0000000000..8ac906ca62
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
index 6114f4fc41..1006382339 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
+++ b/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
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
index 0cbee519e3..2ed07ccf7d 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
+++ b/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>
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index dcae8ae427..73a5b67fbd 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
+++ b/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];
- const point& ccA = meshA.cellCentres()[cellAI];
-
- 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,
- 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 tetPointRef tetB
- (
- ccB,
- tetBasePtB,
- meshB.points()[pt0I],
- meshB.points()[pt1I]
- );
-
- if (!tetB.bounds().overlaps(tetABb))
- {
- continue;
- }
-
- vol += tetTetOverlapVol(tetA, tetB);
-
- if (vol > threshold)
- {
- return true;
- }
- }
- }
- }
- }
-
- return false;
+ hasOverlapOp overlapCheckOp(threshold);
+ cellCellOverlapMinDecomp<hasOverlapOp>
+ (
+ meshA,
+ cellAI,
+ meshB,
+ cellBI,
+ cellBbB,
+ overlapCheckOp
+ );
+
+ return overlapCheckOp.ok_;
+}
+
+
+Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
+(
+ const primitiveMesh& meshA,
+ const label cellAI,
+
+ const primitiveMesh& meshB,
+ const label cellBI,
+ const treeBoundBox& cellBbB
+) const
+{
+ sumOverlapOp overlapSumOp;
+ cellCellOverlapMinDecomp<sumOverlapOp>
+ (
+ meshA,
+ cellAI,
+ meshB,
+ cellBI,
+ cellBbB,
+ overlapSumOp
+ );
+
+ return overlapSumOp.iop_.vol_;
}
@@ -264,116 +125,18 @@ Foam::tetOverlapVolume::cellCellOverlapMomentMinDecomp
const treeBoundBox& cellBbB
) const
{
- const cell& cFacesA = meshA.cells()[cellAI];
- const point& ccA = meshA.cellCentres()[cellAI];
-
- 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,
- 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 tetPointRef tetB
- (
- ccB,
- tetBasePtB,
- meshB.points()[pt0I],
- meshB.points()[pt1I]
- );
- if (!tetB.bounds().overlaps(tetABb))
- {
- continue;
- }
-
- vol += tetTetOverlapVol(tetA, tetB);
- }
- }
- }
- }
-
- return vol;
+ sumOverlapMomentOp overlapSumOp;
+ cellCellOverlapMinDecomp<sumOverlapMomentOp>
+ (
+ meshA,
+ cellAI,
+ meshB,
+ cellBI,
+ cellBbB,
+ overlapSumOp
+ );
+
+ return overlapSumOp.iop_.vol_;
}
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index 06b315d47f..a77e608c6e 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
+++ b/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,69 +54,126 @@ class polyMesh;
class tetOverlapVolume
{
- // Private member functions
-
- //- Tet overlap volume
- scalar tetTetOverlapVol
- (
- const tetPointRef& tetA,
- const tetPointRef& tetB
- ) const;
-
- //- Return a const treeBoundBox
- 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>
+ //- tetPoints handling : sum resulting volumes
+ class sumMomentOp
{
- private:
+ public:
+ Tuple2<scalar, point> vol_;
- //- The number of stored elements
- label n_;
+ 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:
- //- Construct null
- cutTetList()
+ const scalar threshold_;
+ tetPointRef::sumVolOp iop_;
+ bool ok_;
+
+ inline hasOverlapOp(const scalar threshold)
:
- n_(0)
+ threshold_(threshold),
+ iop_(),
+ ok_(false)
{}
- //- Clear the array
- void clear()
+ //- Overlap two tets
+ inline bool operator()(const tetPoints& A, const tetPoints& B)
{
- n_ = 0;
+ tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
+ ok_ = (iop_.vol_ > threshold_);
+ return ok_;
}
+ };
+
+ //- tetPoints combining : sum overlap volume
+ class sumOverlapOp
+ {
+ public:
+
+ tetPointRef::sumVolOp iop_;
- //- Get the current size
- label size() const
+ inline sumOverlapOp()
+ :
+ iop_()
+ {}
+
+ //- Overlap two tets
+ inline bool operator()(const tetPoints& A, const tetPoints& B)
{
- return n_;
+ tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
+ return false;
}
+ };
+
+ //- tetPoints combining : sum overlap volume
+ class sumOverlapMomentOp
+ {
+ public:
+
+ sumMomentOp iop_;
- //- Add a new tet to the end of the array
- void append(const FixedList<point, 4>& t)
+ inline sumOverlapMomentOp()
+ :
+ iop_()
+ {}
+
+ //- Overlap two tets
+ inline bool operator()(const tetPoints& A, const tetPoints& B)
{
- this->operator[](n_) = t;
- ++ n_;
+ tetTetOverlap<sumMomentOp>(A, B, iop_);
+ return false;
}
};
+ // Private member functions
+
+ //- Tet overlap calculation
+ template<class tetPointsOp>
+ static void tetTetOverlap
+ (
+ const tetPoints& tetA,
+ const tetPoints& tetB,
+ 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
+ static treeBoundBox pyrBb
+ (
+ const pointField& points,
+ const face& f,
+ const point& fc
+ );
+
+
public:
//- Runtime type information
@@ -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
// ************************************************************************* //
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
new file mode 100644
index 0000000000..34eded1ff3
--- /dev/null
+++ b/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;
+ }
+ }
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
--
GitLab