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 .
+
+Class
+ Foam::tetPoints
+
+Description
+ Tet storage. Null constructable (unfortunately tetrahedron
+ 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
+{
+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&
);
+typedef tetrahedron tetPointRef;
+
/*---------------------------------------------------------------------------*\
class tetrahedron Declaration
\*---------------------------------------------------------------------------*/
@@ -82,12 +86,78 @@ inline Ostream& operator<<
template
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.
+
+ //- 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&,
+ const tetPoints&,
+ const label,
+ const label
+ );
+
+ template
+ inline static void decomposePrism
+ (
+ const FixedList& points,
+ TetOp& op
+ );
+
+ template
+ 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
+ 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& 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::inside(const point& pt) const
}
+template
+inline void Foam::tetrahedron::dummyOp::operator()
+(
+ const tetPoints&
+)
+{}
+
+
+template
+inline Foam::tetrahedron::sumVolOp::sumVolOp()
+:
+ vol_(0.0)
+{}
+
+
+template
+inline void Foam::tetrahedron::sumVolOp::operator()
+(
+ const tetPoints& tet
+)
+{
+ vol_ += tet.tet().mag();
+}
+
+
+template
+inline Foam::tetrahedron::storeOp::storeOp
+(
+ tetIntersectionList& tets,
+ label& nTets
+)
+:
+ tets_(tets),
+ nTets_(nTets)
+{}
+
+
+template
+inline void Foam::tetrahedron::storeOp::operator()
+(
+ const tetPoints& tet
+)
+{
+ tets_[nTets_++] = tet;
+}
+
+
+template
+inline Foam::point Foam::tetrahedron::planeIntersection
+(
+ const FixedList& d,
+ const tetPoints& t,
+ const label negI,
+ const label posI
+)
+{
+ return
+ (d[posI]*t[negI] - d[negI]*t[posI])
+ / (-d[negI]+d[posI]);
+}
+
+
+template
+template
+inline void Foam::tetrahedron::decomposePrism
+(
+ const FixedList& 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
+template
+inline void Foam::tetrahedron::tetSliceWithPlane
+(
+ const plane& pl,
+ const tetPoints& tet,
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+)
+{
+ // Distance to plane
+ FixedList 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 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 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 p
+ (
+ {
+ tet[0],
+ p02,
+ p03,
+ tet[1],
+ p12,
+ p13
+ }
+ );
+
+ //Pout<< " 01 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList 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 p
+ (
+ {
+ tet[1],
+ p01,
+ p13,
+ tet[2],
+ p02,
+ p23
+ }
+ );
+
+ //Pout<< " 12 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList 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 p
+ (
+ {
+ tet[2],
+ p12,
+ p23,
+ tet[0],
+ p01,
+ p03
+ }
+ );
+
+ //Pout<< " 20 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList 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 p
+ (
+ {
+ tet[3],
+ p23,
+ p13,
+ tet[0],
+ p02,
+ p01
+ }
+ );
+
+ //Pout<< " 03 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList 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 p
+ (
+ {
+ tet[1],
+ p12,
+ p01,
+ tet[3],
+ p23,
+ p03
+ }
+ );
+
+ //Pout<< " 13 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList 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 p
+ (
+ {
+ tet[2],
+ p02,
+ p12,
+ tet[3],
+ p03,
+ p13
+ }
+ );
+
+ //Pout<< " 23 aboveprism:" << p << endl;
+ decomposePrism(p, aboveOp);
+ }
+ {
+ FixedList 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 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 p
+ (
+ {
+ tet[i3],
+ tet[i1],
+ tet[i2],
+ p03,
+ p01,
+ p02
+ }
+ );
+
+ //Pout<< " belowprism:" << p << endl;
+ decomposePrism(p, belowOp);
+ }
+ }
+ else // nPos == 0
+ {
+ belowOp(tet);
+ }
+}
+
+
+template
+template
+inline void Foam::tetrahedron::sliceWithPlane
+(
+ const plane& pl,
+ AboveTetOp& aboveOp,
+ BelowTetOp& belowOp
+) const
+{
+ tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp);
+}
+
+
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template
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 t({tetA.a(), tetA.b(), tetA.c(), tetA.d()});
- cutTetList1.clear();
- tetCut(t, pl0, cut::appendOp(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& t = cutTetList1[i];
- tetCut(t, pl1, cut::appendOp(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& t = cutTetList2[i];
- tetCut(t, pl2, cut::appendOp(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& 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];
+ hasOverlapOp overlapCheckOp(threshold);
+ cellCellOverlapMinDecomp
+ (
+ meshA,
+ cellAI,
+ meshB,
+ cellBI,
+ cellBbB,
+ overlapCheckOp
+ );
- 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());
+ return overlapCheckOp.ok_;
+}
- // Loop over tets of cellB
- forAll(cFacesB, cFB)
- {
- label faceBI = cFacesB[cFB];
+Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
+(
+ const primitiveMesh& meshA,
+ const label cellAI,
- const face& fB = meshB.faces()[faceBI];
- const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
- if (!pyrB.overlaps(pyrA))
- {
- continue;
- }
+ const primitiveMesh& meshB,
+ const label cellBI,
+ const treeBoundBox& cellBbB
+) const
+{
+ sumOverlapOp overlapSumOp;
+ cellCellOverlapMinDecomp
+ (
+ meshA,
+ cellAI,
+ meshB,
+ cellBI,
+ cellBbB,
+ overlapSumOp
+ );
- 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;
+ 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];
+ sumOverlapMomentOp overlapSumOp;
+ cellCellOverlapMinDecomp
+ (
+ meshA,
+ cellAI,
+ meshB,
+ cellBI,
+ cellBbB,
+ overlapSumOp
+ );
- 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;
+ 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,67 +54,124 @@ class polyMesh;
class tetOverlapVolume
{
+ // Private classes
+
+ //- tetPoints handling : sum resulting volumes
+ class sumMomentOp
+ {
+ public:
+ Tuple2 vol_;
+
+ inline sumMomentOp()
+ :
+ vol_(0.0, Zero)
+ {}
+
+ inline void operator()(const tetPoints& tet)
+ {
+ const tetPointRef t(tet.tet());
+ scalar tetVol = t.mag();
+ vol_.first() += tetVol;
+ vol_.second() += (tetVol*t.centre());
+ }
+ };
+
+ //- tetPoints combining : check for overlap
+ class hasOverlapOp
+ {
+ public:
+
+ const scalar threshold_;
+ tetPointRef::sumVolOp iop_;
+ bool ok_;
+
+ inline hasOverlapOp(const scalar threshold)
+ :
+ threshold_(threshold),
+ iop_(),
+ ok_(false)
+ {}
+
+ //- Overlap two tets
+ inline bool operator()(const tetPoints& A, const tetPoints& B)
+ {
+ tetTetOverlap(A, B, iop_);
+ ok_ = (iop_.vol_ > threshold_);
+ return ok_;
+ }
+ };
+
+ //- tetPoints combining : sum overlap volume
+ class sumOverlapOp
+ {
+ public:
+
+ tetPointRef::sumVolOp iop_;
+
+ inline sumOverlapOp()
+ :
+ iop_()
+ {}
+
+ //- Overlap two tets
+ inline bool operator()(const tetPoints& A, const tetPoints& B)
+ {
+ tetTetOverlap(A, B, iop_);
+ return false;
+ }
+ };
+
+ //- tetPoints combining : sum overlap volume
+ class sumOverlapMomentOp
+ {
+ public:
+
+ sumMomentOp iop_;
+
+ inline sumOverlapMomentOp()
+ :
+ iop_()
+ {}
+
+ //- Overlap two tets
+ inline bool operator()(const tetPoints& A, const tetPoints& B)
+ {
+ tetTetOverlap(A, B, iop_);
+ return false;
+ }
+ };
+
+
// Private member functions
- //- Tet overlap volume
- scalar tetTetOverlapVol
+ //- Tet overlap calculation
+ template
+ static void tetTetOverlap
(
- const tetPointRef& tetA,
- const tetPointRef& tetB
- ) const;
+ const tetPoints& tetA,
+ const tetPoints& tetB,
+ tetPointsOp& insideOp
+ );
+
+ //- Cell overlap calculation
+ template
+ static void cellCellOverlapMinDecomp
+ (
+ const primitiveMesh& meshA,
+ const label cellAI,
+ const primitiveMesh& meshB,
+ const label cellBI,
+ const treeBoundBox& cellBbB,
+ tetsOp& combineTetsOp
+ );
//- Return a const treeBoundBox
- treeBoundBox pyrBb
+ static treeBoundBox pyrBb
(
const pointField& points,
const face& f,
const point& fc
- ) const;
-
-
- // Private classes
-
- //- A fixed list of tets which simulates a dynamic list by incrementing
- // a counter whenever its append method is called. This is used as an
- // optimisation for the tetTetOverlapVol method.
- template
- class cutTetList
- :
- public FixedList, Size>
- {
- private:
-
- //- The number of stored elements
- label n_;
-
-
- public:
-
- //- Construct null
- cutTetList()
- :
- n_(0)
- {}
-
- //- Clear the array
- void clear()
- {
- n_ = 0;
- }
-
- //- Get the current size
- label size() const
- {
- return n_;
- }
-
- //- Add a new tet to the end of the array
- void append(const FixedList& t)
- {
- this->operator[](n_) = t;
- ++ n_;
- }
- };
+ );
public:
@@ -159,6 +218,17 @@ public:
const label cellBI,
const treeBoundBox& cellBbB
) const;
+
+ //- Calculates the overlap volume and moment
+ Tuple2 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetOverlapVolume.H"
+#include "primitiveMesh.H"
+
+// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+template
+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 tetAFaceAreas;
+ static FixedList 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 tetBFaceAreas;
+ static FixedList 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
+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;
+ }
+ }
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //