diff --git a/applications/test/momentOfInertia/Test-momentOfInertia.C b/applications/test/momentOfInertia/Test-momentOfInertia.C
index 5abba3d382a4e58894c60905a654f8e27c9488bd..3961caaffc1ade2f72099029e7079f23361cfa14 100644
--- a/applications/test/momentOfInertia/Test-momentOfInertia.C
+++ b/applications/test/momentOfInertia/Test-momentOfInertia.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -35,7 +35,7 @@ Description
#include "polyMesh.H"
#include "ListOps.H"
#include "face.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "triFaceList.H"
#include "OFstream.H"
#include "meshTools.H"
diff --git a/applications/test/tetTetOverlap/Test-tetTetOverlap.C b/applications/test/tetTetOverlap/Test-tetTetOverlap.C
index e31493161484898ea134c47de7bd7581e67b2b11..b85db9c681b79ad8c6b468352e78be0fedc8ac68 100644
--- a/applications/test/tetTetOverlap/Test-tetTetOverlap.C
+++ b/applications/test/tetTetOverlap/Test-tetTetOverlap.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -29,9 +29,10 @@ Description
\*---------------------------------------------------------------------------*/
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "OFstream.H"
#include "meshTools.H"
+#include "cut.H"
using namespace Foam;
@@ -41,7 +42,7 @@ void writeOBJ
(
Ostream& os,
label& vertI,
- const tetPoints& tet
+ const FixedList<point, 4>& tet
)
{
forAll(tet, fp)
@@ -58,105 +59,126 @@ void writeOBJ
}
+tetPointRef makeTetPointRef(const FixedList<point, 4>& p)
+{
+ return tetPointRef(p[0], p[1], p[2], p[3]);
+}
+
+
int main(int argc, char *argv[])
{
- tetPoints A
- (
+ // Tets to test
+ FixedList<point, 4> tetA
+ ({
point(0, 0, 0),
point(1, 0, 0),
point(1, 1, 0),
point(1, 1, 1)
- );
- const tetPointRef tetA = A.tet();
-
- tetPoints B
- (
+ });
+ FixedList<point, 4> tetB
+ ({
point(0.1, 0.1, 0.1),
point(1.1, 0.1, 0.1),
point(1.1, 1.1, 0.1),
point(1.1, 1.1, 1.1)
- );
- const tetPointRef tetB = B.tet();
+ });
- tetPointRef::tetIntersectionList insideTets;
- label nInside = 0;
- tetPointRef::tetIntersectionList outsideTets;
- label nOutside = 0;
+ // Do intersection
+ typedef DynamicList<FixedList<point, 4>> tetList;
+ tetList tetsIn1, tetsIn2, tetsOut;
+ cut::appendOp<tetList> tetOpIn1(tetsIn1);
+ cut::appendOp<tetList> tetOpIn2(tetsIn2);
+ cut::appendOp<tetList> tetOpOut(tetsOut);
- tetA.tetOverlap
- (
- tetB,
- insideTets,
- nInside,
- outsideTets,
- nOutside
- );
+ const plane p0(tetB[1], tetB[3], tetB[2]);
+ tetsIn1.clear();
+ tetCut(tetA, p0, tetOpIn1, tetOpOut);
+ const plane p1(tetB[0], tetB[2], tetB[3]);
+ tetsIn2.clear();
+ forAll(tetsIn1, i)
+ {
+ tetCut(tetsIn1[i], p1, tetOpIn2, tetOpOut);
+ }
- // Dump to file
- // ~~~~~~~~~~~~
+ const plane p2(tetB[0], tetB[3], tetB[1]);
+ tetsIn1.clear();
+ forAll(tetsIn2, i)
+ {
+ tetCut(tetsIn2[i], p2, tetOpIn1, tetOpOut);
+ }
+ const plane p3(tetB[0], tetB[1], tetB[2]);
+ tetsIn2.clear();
+ forAll(tetsIn1, i)
{
- OFstream str("tetA.obj");
+ tetCut(tetsIn1[i], p3, tetOpIn2, tetOpOut);
+ }
+
+ const tetList& tetsIn = tetsIn2;
+
+
+ // Dump to file
+ {
+ OFstream str("A.obj");
Info<< "Writing A to " << str.name() << endl;
label vertI = 0;
- writeOBJ(str, vertI, A);
+ writeOBJ(str, vertI, tetA);
}
{
- OFstream str("tetB.obj");
+ OFstream str("B.obj");
Info<< "Writing B to " << str.name() << endl;
label vertI = 0;
- writeOBJ(str, vertI, B);
+ writeOBJ(str, vertI, tetB);
}
{
- OFstream str("inside.obj");
- Info<< "Writing parts of A inside B to " << str.name() << endl;
+ OFstream str("AInB.obj");
+ Info<< "Writing parts of B inside A to " << str.name() << endl;
label vertI = 0;
- for (label i = 0; i < nInside; ++i)
+ forAll(tetsIn, i)
{
- writeOBJ(str, vertI, insideTets[i]);
+ writeOBJ(str, vertI, tetsIn[i]);
}
}
{
- OFstream str("outside.obj");
- Info<< "Writing parts of A outside B to " << str.name() << endl;
+ OFstream str("AOutB.obj");
+ Info<< "Writing parts of B inside A to " << str.name() << endl;
label vertI = 0;
- for (label i = 0; i < nOutside; ++i)
+ forAll(tetsOut, i)
{
- writeOBJ(str, vertI, outsideTets[i]);
+ writeOBJ(str, vertI, tetsOut[i]);
}
}
- // Check
- // ~~~~~
+ // Check the volumes
+ Info<< "Vol A: " << makeTetPointRef(tetA).mag() << endl;
- Info<< "Vol A:" << tetA.mag() << endl;
-
- scalar volInside = 0;
- for (label i = 0; i < nInside; ++i)
+ scalar volIn = 0;
+ forAll(tetsIn, i)
{
- volInside += insideTets[i].tet().mag();
+ volIn += makeTetPointRef(tetsIn[i]).mag();
}
- Info<< "Vol A inside B:" << volInside << endl;
+ Info<< "Vol A inside B: " << volIn << endl;
- scalar volOutside = 0;
- for (label i = 0; i < nOutside; ++i)
+ scalar volOut = 0;
+ forAll(tetsOut, i)
{
- volOutside += outsideTets[i].tet().mag();
+ volOut += makeTetPointRef(tetsOut[i]).mag();
}
- Info<< "Vol A outside B:" << volOutside << endl;
+ Info<< "Vol A outside B: " << volOut << endl;
- Info<< "Sum inside and outside:" << volInside+volOutside << endl;
+ Info<< "Sum inside and outside: " << volIn + volOut << endl;
- if (mag(volInside+volOutside-tetA.mag()) > SMALL)
+ if (mag(volIn + volOut - makeTetPointRef(tetA).mag()) > SMALL)
{
FatalErrorInFunction
<< "Tet volumes do not sum up to input tet."
<< exit(FatalError);
}
+
return 0;
}
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
index a36e0370cfb0d12b097b941af3598c6cac477298..1ab57746900c2ac4fb79a2148e32d22e7cacf88d 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -28,7 +28,7 @@ License
#include "pointIOField.H"
#include "scalarIOField.H"
#include "triadIOField.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "plane.H"
#include "transform.H"
#include "meshTools.H"
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/fileControl/fileControl.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/fileControl/fileControl.C
index a644134066fe54f1f7ab9c3aa9083ac4a37e715d..8f8aac17e988dc52095754762bbe7b6edcc99b31 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/fileControl/fileControl.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/fileControl/fileControl.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,7 +25,7 @@ License
#include "fileControl.H"
#include "addToRunTimeSelectionTable.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "scalarList.H"
#include "vectorTools.H"
#include "pointIOField.H"
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
index ef565a405f0b33a422d0b3b7acf384114e5d43f7..c0782e3a3c712e2d7e5c906a670a8ace0fef4fdc 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -28,7 +28,7 @@ License
#include "cellSizeFunction.H"
#include "triSurfaceMesh.H"
#include "searchableBox.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "vectorTools.H"
#include "quaternion.H"
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellChecks.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellChecks.C
index 49f6796781d235b8690a6a41819850e32a65cf06..570ad3e1845e0a15db87522c2adbfcc9e10b9e2c 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellChecks.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellChecks.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -24,7 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "plane.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "pointConversion.H"
#include "CGALTriangulation3DKernel.H"
diff --git a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H
index 18f1459f4e8565cb745b6fdd0c0bf51dc36f62fe..8669efbc3c3e8b285d54fa214e86cf28ec9193e4 100644
--- a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H
+++ b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -43,7 +43,7 @@ SourceFiles
#include "triFace.H"
#include "edge.H"
#include "pointField.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
index 0a3ed7a03b8624dbf3504abaf028d71559ff6c8d..18de73777c620b5bf63554d28e8539389ae48966 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -42,7 +42,7 @@ SourceFiles
#include "polyMesh.H"
#include "coupledPolyPatch.H"
#include "syncTools.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "tetIndices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
index 0c80210e636e62b6a25ee1335ea2efc7ae411d2f..8ad320403a8b22b00ed14a381a1c8ebc9d77b8a5 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -56,7 +56,7 @@ SourceFiles
#define tetIndices_H
#include "label.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "triPointRef.H"
#include "polyMesh.H"
#include "triFace.H"
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
new file mode 100644
index 0000000000000000000000000000000000000000..9e821c0c536050045c73445f6371e7416c92fa3e
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
@@ -0,0 +1,539 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+Description:
+ Functions which cut triangles and tetrahedra. Generic operations are
+ applied to each half of a cut.
+
+SourceFiles:
+ cutI.H
+ cutTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cut_H
+#define cut_H
+
+#include "FixedList.H"
+#include "nil.H"
+#include "plane.H"
+#include "tetPointRef.H"
+#include "triPointRef.H"
+#include "zero.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cut
+{
+
+/*---------------------------------------------------------------------------*\
+ Class uniformOp Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class uniformOp
+{
+private:
+
+ //- Data
+ Type data_;
+
+
+public:
+
+ // Constructors
+
+ //- Construct null
+ uniformOp()
+ {}
+
+ //- Construct from data
+ uniformOp(Type data)
+ :
+ data_(data)
+ {}
+
+
+ // Member functions
+
+ //- Access the data
+ Type data() const
+ {
+ return data_;
+ }
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class noOp Declaration
+\*---------------------------------------------------------------------------*/
+
+class noOp
+:
+ public uniformOp<nil>
+{
+public:
+
+ // Typedefs
+
+ //- Result type
+ typedef zero result;
+
+
+ // Constructors
+
+ //- Construct null
+ noOp()
+ {}
+
+ //- Construct from base
+ noOp(const uniformOp<nil>& op)
+ :
+ uniformOp<nil>(op)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return result();
+ }
+
+ //- Operate on a triangle or tetrahedron
+ template<unsigned Size>
+ result operator()(const FixedList<point, Size>& p) const
+ {
+ return result();
+ }
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class areaOp Declaration
+\*---------------------------------------------------------------------------*/
+
+class areaOp
+:
+ public uniformOp<nil>
+{
+public:
+
+ // Typedefs
+
+ //- Result type
+ typedef vector result;
+
+
+ // Constructors
+
+ //- Construct null
+ areaOp()
+ {}
+
+ //- Construct from base
+ areaOp(const uniformOp<nil>& op)
+ :
+ uniformOp<nil>(op)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return vector::zero;
+ }
+
+ //- Operate on a triangle
+ result operator()(const FixedList<point, 3>& p) const
+ {
+ return result(triPointRef(p[0], p[1], p[2]).normal());
+ }
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class volumeOp Declaration
+\*---------------------------------------------------------------------------*/
+
+class volumeOp
+:
+ public uniformOp<nil>
+{
+public:
+
+ // Typedefs
+
+ //- Result type
+ typedef scalar result;
+
+
+ // Constructors
+
+ //- Construct null
+ volumeOp()
+ {}
+
+ //- Construct from base
+ volumeOp(const uniformOp<nil>& op)
+ :
+ uniformOp<nil>(op)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return 0;
+ }
+
+ //- Operate on a tetrahedron
+ result operator()(const FixedList<point, 4>& p) const
+ {
+ return result(tetPointRef(p[0], p[1], p[2], p[3]).mag());
+ }
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class areaIntegrateOp Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class areaIntegrateOp
+:
+ public FixedList<Type, 3>
+{
+public:
+
+ // Typedefs
+
+ //- Result type
+ typedef typename outerProduct<Type, vector>::type result;
+
+
+ // Constructors
+
+ //- Construct from base
+ areaIntegrateOp(const FixedList<Type, 3>& x)
+ :
+ FixedList<Type, 3>(x)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return pTraits<result>::zero;
+ }
+
+ //- Operate on a triangle
+ result operator()(const FixedList<point, 3>& p) const
+ {
+ const FixedList<Type, 3>& x = *this;
+ return result(areaOp()(p)*(x[0] + x[1] + x[2])/3);
+ }
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class volumeIntegrateOp Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class volumeIntegrateOp
+:
+ public FixedList<Type, 4>
+{
+public:
+
+ // Typedefs
+
+ //- Result type
+ typedef Type result;
+
+
+ // Constructors
+
+ //- Construct from base
+ volumeIntegrateOp(const FixedList<Type, 4>& x)
+ :
+ FixedList<Type, 4>(x)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return pTraits<result>::zero;
+ }
+
+ //- Operate on a tetrahedron
+ result operator()(const FixedList<point, 4>& p) const
+ {
+ const FixedList<Type, 4>& x = *this;
+ return result(volumeOp()(p)*(x[0] + x[1] + x[2] + x[3])/4);
+ }
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class listOp Declaration
+\*---------------------------------------------------------------------------*/
+
+template<unsigned Size>
+class listOp
+:
+ public uniformOp<nil>
+{
+public:
+
+ // Classes
+
+ //- Result class
+ class result
+ :
+ public DynamicList<FixedList<point, Size>>
+ {
+ public:
+
+ // Constructors
+
+ //- Construct from a single element
+ result(const FixedList<point, Size>& x)
+ :
+ DynamicList<FixedList<point, Size>>(1, x)
+ {}
+
+
+ // Member operators
+
+ //- Add together two lists
+ result operator+(const result& x) const
+ {
+ result r(*this);
+ r.append(x);
+ return r;
+ }
+ };
+
+
+ // Constructors
+
+ //- Construct null
+ listOp()
+ {}
+
+ //- Construct from base
+ listOp(const uniformOp<nil>& op)
+ :
+ uniformOp<nil>(op)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return result();
+ }
+
+ //- Operate on a triangle or tetrahedron
+ result operator()(const FixedList<point, Size>& p) const
+ {
+ return result(p);
+ }
+};
+
+
+typedef listOp<3> listTriOp;
+
+
+typedef listOp<4> listTetOp;
+
+
+/*---------------------------------------------------------------------------*\
+ Class appendOp Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Container>
+class appendOp
+:
+ public uniformOp<Container&>
+{
+public:
+
+ // Typedefs
+
+ //- Result type
+ typedef zero result;
+
+
+ // Constructors
+
+ //- Construct from a container reference
+ appendOp(Container& x)
+ :
+ uniformOp<Container&>(x)
+ {}
+
+ //- Construct from base
+ appendOp(const uniformOp<Container&>& op)
+ :
+ uniformOp<Container&>(op)
+ {}
+
+
+ // Member operators
+
+ //- Operate on nothing
+ result operator()() const
+ {
+ return result();
+ }
+
+ //- Operate on a triangle or tetrahedron
+ template<unsigned Size>
+ result operator()(const FixedList<point, Size>& p) const
+ {
+ this->data().append(p);
+ return result();
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Trait to determine the result of the addition of two operations
+template<class AheadOp, class BehindOp>
+class opAddResult;
+
+template<class Op>
+class opAddResult<Op, Op>
+{
+public:
+
+ typedef typename Op::result type;
+};
+
+template<>
+class opAddResult<noOp, noOp>
+{
+public:
+
+ typedef typename noOp::result type;
+};
+
+template<class Op>
+class opAddResult<noOp, Op>
+{
+public:
+
+ typedef typename Op::result type;
+};
+
+template<class Op>
+class opAddResult<Op, noOp>
+{
+public:
+
+ typedef typename Op::result type;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace cut
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Cut a triangle along the zero plane defined by the given levels. Templated
+// on aboveOp and belowOp; the operations applied to the regions on either side
+// of the cut.
+template<class AboveOp, class BelowOp>
+typename cut::opAddResult<AboveOp, BelowOp>::type triCut
+(
+ const FixedList<point, 3>& tri,
+ const FixedList<scalar, 3>& level,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+);
+
+//- As above, but with a plane specifying the location of the cut
+template<class AboveOp, class BelowOp>
+typename cut::opAddResult<AboveOp, BelowOp>::type triCut
+(
+ const FixedList<point, 3>& tri,
+ const plane& s,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+);
+
+//- As triCut, but for a tetrahedron.
+template<class AboveOp, class BelowOp>
+typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
+(
+ const FixedList<point, 4>& tet,
+ const FixedList<scalar, 4>& level,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+);
+
+//- As above, but with a plane specifying the location of the cut
+template<class AboveOp, class BelowOp>
+typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
+(
+ const FixedList<point, 4>& tet,
+ const plane& s,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "cutI.H"
+
+#ifdef NoRepository
+ #include "cutTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
new file mode 100644
index 0000000000000000000000000000000000000000..aefc537f519929d3230f95d97ce9ad3b8aeecb51
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
@@ -0,0 +1,400 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cut.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Modify a uniform operation for reordering a tri (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& triReorder
+(
+ const cut::uniformOp<Type>& x,
+ const FixedList<label, 3>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for cutting a tri from a tri (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& triCutTri
+(
+ const cut::uniformOp<Type>& x,
+ const Pair<scalar>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for cutting a quad from a tri (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& triCutQuad
+(
+ const cut::uniformOp<Type>& x,
+ const Pair<scalar>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for reordering a tet (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& tetReorder
+(
+ const cut::uniformOp<Type>& x,
+ const FixedList<label, 4>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for cutting a tet from a tet (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& tetCutTet
+(
+ const cut::uniformOp<Type>& x,
+ const FixedList<scalar, 3>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for cutting prism0 from a tet (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& tetCutPrism0
+(
+ const cut::uniformOp<Type>& x,
+ const FixedList<scalar, 3>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for cutting prism01 from a tet (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& tetCutPrism01
+(
+ const cut::uniformOp<Type>& x,
+ const FixedList<scalar, 4>&
+)
+{
+ return x;
+}
+
+
+//- Modify a uniform operation for cutting prism23 from a tet (does nothing)
+template<class Type>
+inline const cut::uniformOp<Type>& tetCutPrism23
+(
+ const cut::uniformOp<Type>& x,
+ const FixedList<scalar, 4>&
+)
+{
+ return x;
+}
+
+
+//- Modify a fixed list for reordering a tri (does nothing)
+template<class Type, unsigned Size>
+inline FixedList<Type, 3> triReorder
+(
+ const FixedList<Type, Size>& x,
+ const FixedList<label, 3>& indices
+)
+{
+ FixedList<Type, 3> result;
+ for (unsigned i = 0; i < 3; ++ i)
+ {
+ result[i] = x[indices[i]];
+ }
+ return result;
+}
+
+
+//- Modify a list for cutting a tri from a tri
+template<class Type>
+inline FixedList<Type, 3> triCutTri
+(
+ const FixedList<Type, 3>& x,
+ const Pair<scalar>& f
+)
+{
+ FixedList<Type, 3> result;
+ result[0] = x[0];
+ for (label i = 0; i < 2; ++ i)
+ {
+ result[i+1] = x[0] + f[i]*(x[i+1] - x[0]);
+ }
+ return result;
+}
+
+
+//- Modify a list for cutting a quad from a tri
+template<class Type>
+inline FixedList<Type, 4> triCutQuad
+(
+ const FixedList<Type, 3>& x,
+ const Pair<scalar>& f
+)
+{
+ FixedList<Type, 4> result;
+ for (label i = 0; i < 2; ++ i)
+ {
+ result[i] = x[i+1];
+ result[3-i] = x[0] + f[i]*(x[i+1] - x[0]);
+ }
+ return result;
+}
+
+
+//- Modify a fixed list for reordering a tet (does nothing)
+template<class Type, unsigned Size>
+inline FixedList<Type, 4> tetReorder
+(
+ const FixedList<Type, Size>& x,
+ const FixedList<label, 4>& indices
+)
+{
+ FixedList<Type, 4> result;
+ for (unsigned i = 0; i < 4; ++ i)
+ {
+ result[i] = x[indices[i]];
+ }
+ return result;
+}
+
+
+//- Modify a list for cutting a tet from a tet
+template<class Type>
+inline FixedList<Type, 4> tetCutTet
+(
+ const FixedList<Type, 4>& x,
+ const FixedList<scalar, 3>& f
+)
+{
+ FixedList<Type, 4> result;
+ result[0] = x[0];
+ for (label i = 0; i < 3; ++ i)
+ {
+ result[i+1] = x[0] + f[i]*(x[i+1] - x[0]);
+ }
+ return result;
+}
+
+
+//- Modify a list for cutting prism0 from a tet
+template<class Type>
+inline FixedList<Type, 6> tetCutPrism0
+(
+ const FixedList<Type, 4>& x,
+ const FixedList<scalar, 3>& f
+)
+{
+ FixedList<Type, 6> result;
+ for (label i = 0; i < 3; ++ i)
+ {
+ result[i] = x[0] + f[i]*(x[i+1] - x[0]);
+ result[i+3] = x[i+1];
+ }
+ return result;
+}
+
+
+//- Modify a list for cutting prism01 from a tet
+template<class Type>
+inline FixedList<Type, 6> tetCutPrism01
+(
+ const FixedList<Type, 4>& x,
+ const FixedList<scalar, 4>& f
+)
+{
+ FixedList<Type, 6> result;
+ for (label i = 0; i < 2; ++ i)
+ {
+ result[3*i] = x[i];
+ for (label j = 0; j < 2; ++ j)
+ {
+ result[3*i+j+1] = x[i] + f[2*i+j]*(x[j+2] - x[i]);
+ }
+ }
+
+ return result;
+}
+
+
+//- Modify a list for cutting prism23 from a tet
+template<class Type>
+inline FixedList<Type, 6> tetCutPrism23
+(
+ const FixedList<Type, 4>& x,
+ const FixedList<scalar, 4>& f
+)
+{
+ FixedList<Type, 6> result = tetCutPrism01(x, f);
+ result[0] = x[2];
+ result[3] = x[3];
+ Swap(result[2], result[4]);
+ return result;
+}
+
+
+//- Cut a tri from a tri and apply an operation to the result. The cut is made
+// along the two edges connected to vertex 0, and the cut locations are given
+// as factors along these edges. The result is the side connected to vertex 0.
+template<class Op>
+inline typename Op::result triCutTri
+(
+ const Op& op,
+ const FixedList<point, 3>& p,
+ const Pair<scalar>& f
+)
+{
+ return Op(triCutTri(op, f))(triCutTri(p, f));
+}
+
+
+//- Apply an operation to a quad. Splits the quad into two tris.
+template<class Op, class OpData>
+inline typename Op::result quadOp
+(
+ const OpData& opData,
+ const FixedList<point, 4>& p
+)
+{
+ static const FixedList<FixedList<label, 3>, 2> i =
+ {{0, 1, 2}, {0, 2, 3}};
+ return
+ Op(triReorder(opData, i[0]))(triReorder(p, i[0]))
+ + Op(triReorder(opData, i[1]))(triReorder(p, i[1]));
+}
+
+
+//- Cut a quad from a tri and apply an operation to the result. The cuts are
+// the same as for triCutTri. The result is the side connected to vertices 1
+// and 2.
+template<class Op>
+inline typename Op::result triCutQuad
+(
+ const Op& op,
+ const FixedList<point, 3>& p,
+ const FixedList<scalar, 2>& f
+)
+{
+ return quadOp<Op>(triCutQuad(op, f), triCutQuad(p, f));
+}
+
+
+//- Cut a tet from a tet and apply an operation to the result. The cut is made
+// along the three edges connected to vertex 0, and the cut locations are given
+// as factors along these edges. The result is the side connected to vertex 0.
+template<class Op>
+inline typename Op::result tetCutTet
+(
+ const Op& op,
+ const FixedList<point, 4>& p,
+ const FixedList<scalar, 3>& f
+)
+{
+ return Op(tetCutTet(op, f))(tetCutTet(p, f));
+}
+
+
+//- Apply an operation to a prism. Splits the prism into three tets.
+template<class Op, class OpData>
+inline typename Op::result prismOp
+(
+ const OpData& opData,
+ const FixedList<point, 6>& p
+)
+{
+ static const FixedList<FixedList<label, 4>, 3> i =
+ {{0, 1, 2, 4}, {0, 2, 5, 4}, {0, 4, 5, 3}};
+ return
+ Op(tetReorder(opData, i[0]))(tetReorder(p, i[0]))
+ + Op(tetReorder(opData, i[1]))(tetReorder(p, i[1]))
+ + Op(tetReorder(opData, i[2]))(tetReorder(p, i[2]));
+}
+
+
+//- Cut a prism from a tet and apply an operation to the result. The cuts are
+// the same as for tetCutTet. The result is the side connected to vertices 1,
+// 2 and 3.
+template<class Op>
+inline typename Op::result tetCutPrism0
+(
+ const Op& op,
+ const FixedList<point, 4>& p,
+ const FixedList<scalar, 3>& f
+)
+{
+ return prismOp<Op>(tetCutPrism0(op, f), tetCutPrism0(p, f));
+}
+
+
+//- Cut a prism from a tet and apply an operation to the result. The cut is made
+// along four edges, not edges 01 or 23, and the cut locations are given as
+// factors along these edges. The result is the side connected to edge 01.
+template<class Op>
+inline typename Op::result tetCutPrism01
+(
+ const Op& op,
+ const FixedList<point, 4>& p,
+ const FixedList<scalar, 4>& f
+)
+{
+ return prismOp<Op>(tetCutPrism01(op, f), tetCutPrism01(p, f));
+}
+
+
+//- Cut a prism from a tet and apply an operation to the result. The cuts are
+// the same as for tetCutPrism01. The result is the side connected to edge 23.
+template<class Op>
+inline typename Op::result tetCutPrism23
+(
+ const Op& op,
+ const FixedList<point, 4>& p,
+ const FixedList<scalar, 4>& f
+)
+{
+ return prismOp<Op>(tetCutPrism23(op, f), tetCutPrism23(p, f));
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..b4dda7f4c04be47af17a1cfa824e058c940d0bf2
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cut.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class AboveOp, class BelowOp>
+typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
+(
+ const FixedList<point, 3>& tri,
+ const FixedList<scalar, 3>& level,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+)
+{
+ // If everything is positive or negative, then process the triangle as a
+ // whole, and do a quick return
+ if (level[0] >= 0 && level[1] >= 0 && level[2] >= 0)
+ {
+ return aboveOp(tri) + belowOp();
+ }
+ if (level[0] <= 0 && level[1] <= 0 && level[2] <= 0)
+ {
+ return aboveOp() + belowOp(tri);
+ }
+
+ // There will be just one edge without a sign change. Find it, and put it
+ // opposite the first vertex. This may change the sign of the tri.
+ FixedList<label, 3> indices({0, 1, 2});
+ label i;
+ for (i = 0; i < 3; ++ i)
+ {
+ if (level[(i + 1)%3]*level[(i + 2)%3] >= 0)
+ {
+ Swap(indices[0], indices[i]);
+ break;
+ }
+ }
+ if (i == 3)
+ {
+ FatalErrorInFunction
+ << "The number of tri vertices above the level set should always "
+ << "be 1" << exit(FatalError);
+ }
+
+ // Correct the sign
+ if (indices[0] != 0)
+ {
+ Swap(indices[1], indices[2]);
+ }
+
+ // Permute the data
+ const FixedList<point, 3> p = triReorder(tri, indices);
+ const FixedList<scalar, 3> l = triReorder(level, indices);
+ AboveOp a = triReorder(aboveOp, indices);
+ BelowOp b = triReorder(belowOp, indices);
+
+ // Slice off one corner to form a tri and a quad
+ Pair<scalar> f;
+ for (label i = 0; i < 2; ++ i)
+ {
+ f[i] = l[0]/(l[0] - l[i+1]);
+ }
+ if (l[0] > 0)
+ {
+ return triCutTri(a, p, f) + triCutQuad(b, p, f);
+ }
+ else
+ {
+ return triCutQuad(a, p, f) + triCutTri(b, p, f);
+ }
+}
+
+
+template<class AboveOp, class BelowOp>
+typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
+(
+ const FixedList<point, 3>& tri,
+ const plane& p,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+)
+{
+ // Set the level set to the signed distance from the plane
+ FixedList<scalar, 3> level;
+ for (label i = 0; i < 3; ++ i)
+ {
+ level[i] = (tri[i] - p.refPoint()) & p.normal();
+ }
+
+ // Run the level set function
+ return triCut(tri, level, aboveOp, belowOp);
+}
+
+
+template<class AboveOp, class BelowOp>
+typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
+(
+ const FixedList<point, 4>& tet,
+ const FixedList<scalar, 4>& level,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+)
+{
+ // Get the min and max over all four vertices and quick return if there is
+ // no change of sign
+ scalar levelMin = VGREAT, levelMax = - VGREAT;
+ for (label i = 0; i < 4; ++ i)
+ {
+ levelMin = min(levelMin, level[i]);
+ levelMax = max(levelMax, level[i]);
+ }
+ if (levelMin >= 0)
+ {
+ return aboveOp(tet) + belowOp();
+ }
+ if (levelMax <= 0)
+ {
+ return aboveOp() + belowOp(tet);
+ }
+
+ // Partition the level so that positive values are at the start. This is
+ // like a single iteration of quick-sort, except that the pivot is a hard-
+ // coded zero, rather than an element of the array. This can change the sign
+ // of the tet.
+ FixedList<label, 4> indices({0, 1, 2, 3});
+ bool signChange = false;
+ label i = 0, j = 3;
+ while (true)
+ {
+ while (i < j && level[indices[i]] > 0)
+ {
+ i ++;
+ }
+ while (j > i && level[indices[j]] <= 0)
+ {
+ j --;
+ }
+ if (i == j)
+ {
+ break;
+ }
+ Swap(indices[i], indices[j]);
+ signChange = !signChange;
+ }
+
+ // The number of vertices above the slice
+ label n = i;
+
+ // If there are more positives than negatives then reverse the order so that
+ // the negatives are at the start
+ if (n > 2)
+ {
+ n = 4 - n;
+ for (label i = 0; i < 2; ++ i)
+ {
+ Swap(indices[i], indices[3-i]);
+ }
+ }
+
+ // Correct the sign
+ if (signChange)
+ {
+ Swap(indices[2], indices[3]);
+ }
+
+ // Permute the data
+ const FixedList<point, 4> p = tetReorder(tet, indices);
+ const FixedList<scalar, 4> l = tetReorder(level, indices);
+ AboveOp a = tetReorder(aboveOp, indices);
+ BelowOp b = tetReorder(belowOp, indices);
+
+ // Calculate the integrals above and below the level set
+ if (n == 1)
+ {
+ // Slice off one corner to form a tet and a prism
+ FixedList<scalar, 3> f;
+ for (label i = 0; i < 3; ++ i)
+ {
+ f[i] = l[0]/(l[0] - l[i+1]);
+ }
+ if (l[0] > 0)
+ {
+ return tetCutTet(a, p, f) + tetCutPrism0(b, p, f);
+ }
+ else
+ {
+ return tetCutPrism0(a, p, f) + tetCutTet(b, p, f);
+ }
+ }
+ else if (n == 2)
+ {
+ // Slice off two corners to form two prisms
+ FixedList<scalar, 4> f;
+ for (label i = 0; i < 2; ++ i)
+ {
+ for (label j = 0; j < 2; ++ j)
+ {
+ f[2*i+j] = l[i]/(l[i] - l[j+2]);
+ }
+ }
+ if (l[0] > 0)
+ {
+ return tetCutPrism01(a, p, f) + tetCutPrism23(b, p, f);
+ }
+ else
+ {
+ return tetCutPrism23(a, p, f) + tetCutPrism01(b, p, f);
+ }
+ }
+
+ FatalErrorInFunction
+ << "The number of tet vertices above the level set should always be "
+ << "either 1 or 2" << exit(FatalError);
+
+ return aboveOp() + belowOp();
+}
+
+
+template<class AboveOp, class BelowOp>
+typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
+(
+ const FixedList<point, 4>& tet,
+ const plane& p,
+ const AboveOp& aboveOp,
+ const BelowOp& belowOp
+)
+{
+ // Set the level set to the signed distance from the plane
+ FixedList<scalar, 4> level;
+ for (label i = 0; i < 4; ++ i)
+ {
+ level[i] = (tet[i] - p.refPoint()) & p.normal();
+ }
+
+ // Run the level set function
+ return tetCut(tet, level, aboveOp, belowOp);
+}
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H
new file mode 100644
index 0000000000000000000000000000000000000000..8854eb46fb69bd3e2a3224d499c7f853dae048cb
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2017 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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
deleted file mode 100644
index 8ac906ca62072f5bf373f52d124dad126a75a3cf..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
+++ /dev/null
@@ -1,113 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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.C b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
index 163476a6b32e5166bca73e7e1402e17d9240d99a..45b10041f32c7f297d443ca83f4027111b7b75d5 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-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -29,96 +29,6 @@ License
// * * * * * * * * * * * * * * * 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);
- }
- }
-}
-
-
template<class Point, class PointRef>
Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
(
@@ -426,4 +336,16 @@ void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
}
+template<class Point, class PointRef>
+Foam::boundBox Foam::tetrahedron<Point, PointRef>::bounds() const
+{
+ return
+ boundBox
+ (
+ min(a(), min(b(), min(c(), d()))),
+ max(a(), max(b(), max(c(), d())))
+ );
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
index 2369e77ca81ced17a9d7d162bce139dd4432a45f..27e1ad63c306cf80add94460191037874bbaf5c7 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-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -46,6 +46,7 @@ SourceFiles
#include "FixedList.H"
#include "UList.H"
#include "triPointRef.H"
+#include "boundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -54,8 +55,6 @@ namespace Foam
class Istream;
class Ostream;
-class tetPoints;
-class plane;
// Forward declaration of friend functions and operators
@@ -75,8 +74,6 @@ inline Ostream& operator<<
const tetrahedron<Point, PointRef>&
);
-typedef tetrahedron<point, const point&> tetPointRef;
-
/*---------------------------------------------------------------------------*\
class tetrahedron Declaration
\*---------------------------------------------------------------------------*/
@@ -84,78 +81,12 @@ typedef tetrahedron<point, const point&> tetPointRef;
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:
@@ -254,26 +185,6 @@ 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
@@ -294,6 +205,9 @@ public:
void gradNiGradNj(tensorField& buffer) const;
+ //- Calculate the bounding box
+ boundBox bounds() const;
+
// IOstream operators
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
index b98fc15bb8088f36443445877435677f365ed3cf..8661a5487ee00307119e5bb0e7552fc335923673 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,7 +25,6 @@ License
#include "triangle.H"
#include "IOstreams.H"
-#include "tetPoints.H"
#include "plane.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -490,470 +489,6 @@ 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;
- 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;
- 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;
- 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
- {
- 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;
- 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>
diff --git a/src/OpenFOAM/primitives/zero/zeroI.H b/src/OpenFOAM/primitives/zero/zeroI.H
index 0b7b75f23766223f97e052a3bdc75f1a4aa77412..d69255f1599a4bd1ed018bdedbc91cc49c77dbe6 100644
--- a/src/OpenFOAM/primitives/zero/zeroI.H
+++ b/src/OpenFOAM/primitives/zero/zeroI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -32,6 +32,11 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+inline zero operator+(const zero&, const zero&)
+{
+ return Zero;
+}
+
template<class Type>
inline const Type& operator+(const Type& t, const zero&)
{
@@ -44,6 +49,11 @@ inline const Type& operator+(const zero&, const Type& t)
return t;
}
+inline zero operator-(const zero&, const zero&)
+{
+ return Zero;
+}
+
template<class Type>
inline const Type& operator-(const Type& t, const zero&)
{
@@ -56,6 +66,11 @@ inline Type operator-(const zero&, const Type& t)
return -t;
}
+inline zero operator*(const zero&, const zero&)
+{
+ return Zero;
+}
+
template<class Type>
inline zero operator*(const Type& t, const zero&)
{
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index b7cd414123fb3851d63a44e352c1e50d8217ea42..7c580ae692fc808c6708582a5a55853758c78bc1 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -26,7 +26,7 @@ License
#include "polyMeshGeometry.H"
#include "polyMeshTetDecomposition.H"
#include "pyramidPointFaceRef.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "syncTools.H"
#include "unitConversion.H"
#include "primitiveMeshTools.H"
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index ea678968ce23118d56c1c0e49bec1d0e956e33ca..24355e3a59353118ce0ae36b74397e777b1efb66 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -40,7 +40,7 @@ Description
#include "pointField.H"
#include "faceList.H"
#include "OFstream.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
#include "FixedList.H"
#include "polyMeshTetDecomposition.H"
#include "particleMacros.H"
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index 73a5b67fbd9046b4fc822dda311bd4e8eaef5de5..dcae8ae427c2a7d54e8323d1b2276ff2997d0707 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-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@@ -24,13 +24,12 @@ 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 * * * * * * * * * * * * * //
@@ -48,6 +47,67 @@ 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,
@@ -74,43 +134,122 @@ bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp
const scalar threshold
) const
{
- 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_;
+ 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;
}
@@ -125,18 +264,116 @@ Foam::tetOverlapVolume::cellCellOverlapMomentMinDecomp
const treeBoundBox& cellBbB
) const
{
- sumOverlapMomentOp overlapSumOp;
- cellCellOverlapMinDecomp<sumOverlapMomentOp>
- (
- meshA,
- cellAI,
- meshB,
- cellBI,
- cellBbB,
- overlapSumOp
- );
-
- return overlapSumOp.iop_.vol_;
+ 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;
}
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index a77e608c6e926a455275db78ae34225db50fbf59..06b315d47f96ff00547491abfbaf87f6adc71949 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-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -29,7 +29,6 @@ Description
SourceFiles
tetOverlapVolume.C
- tetOverlapVolumeTemplates.C
\*---------------------------------------------------------------------------*/
@@ -39,8 +38,7 @@ SourceFiles
#include "FixedList.H"
#include "labelList.H"
#include "treeBoundBox.H"
-#include "Tuple2.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
namespace Foam
{
@@ -54,126 +52,69 @@ class polyMesh;
class tetOverlapVolume
{
- // Private classes
+ // Private member functions
- //- tetPoints handling : sum resulting volumes
- class sumMomentOp
- {
- public:
- Tuple2<scalar, point> vol_;
+ //- Tet overlap volume
+ scalar tetTetOverlapVol
+ (
+ const tetPointRef& tetA,
+ const tetPointRef& tetB
+ ) const;
- inline sumMomentOp()
- :
- vol_(0.0, Zero)
- {}
+ //- Return a const treeBoundBox
+ treeBoundBox pyrBb
+ (
+ const pointField& points,
+ const face& f,
+ const point& fc
+ ) const;
- 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
+ // 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>
{
- public:
+ private:
- const scalar threshold_;
- tetPointRef::sumVolOp iop_;
- bool ok_;
+ //- The number of stored elements
+ label n_;
- inline hasOverlapOp(const scalar threshold)
- :
- threshold_(threshold),
- iop_(),
- ok_(false)
- {}
-
- //- Overlap two tets
- inline bool operator()(const tetPoints& A, const tetPoints& B)
- {
- tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
- ok_ = (iop_.vol_ > threshold_);
- return ok_;
- }
- };
- //- tetPoints combining : sum overlap volume
- class sumOverlapOp
- {
public:
- tetPointRef::sumVolOp iop_;
-
- inline sumOverlapOp()
+ //- Construct null
+ cutTetList()
:
- iop_()
+ n_(0)
{}
- //- Overlap two tets
- inline bool operator()(const tetPoints& A, const tetPoints& B)
+ //- Clear the array
+ void clear()
{
- tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
- return false;
+ n_ = 0;
}
- };
-
- //- tetPoints combining : sum overlap volume
- class sumOverlapMomentOp
- {
- public:
- sumMomentOp iop_;
-
- inline sumOverlapMomentOp()
- :
- iop_()
- {}
+ //- Get the current size
+ label size() const
+ {
+ return n_;
+ }
- //- Overlap two tets
- inline bool operator()(const tetPoints& A, const tetPoints& B)
+ //- Add a new tet to the end of the array
+ void append(const FixedList<point, 4>& t)
{
- tetTetOverlap<sumMomentOp>(A, B, iop_);
- return false;
+ this->operator[](n_) = t;
+ ++ n_;
}
};
- // 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
@@ -218,17 +159,6 @@ 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;
};
@@ -238,12 +168,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#ifdef NoRepository
-# include "tetOverlapVolumeTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
#endif
// ************************************************************************* //
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
deleted file mode 100644
index 34eded1ff313478509b37fc627b996297d358213..0000000000000000000000000000000000000000
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
+++ /dev/null
@@ -1,259 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / 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;
- }
- }
- }
- }
- }
-}
-
-
-// ************************************************************************* //