diff --git a/applications/test/momentOfInertia/Test-momentOfInertia.C b/applications/test/momentOfInertia/Test-momentOfInertia.C index 5abba3d382..3961caaffc 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 e314931614..b85db9c681 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& tet ) { forAll(tet, fp) @@ -58,105 +59,126 @@ void writeOBJ } +tetPointRef makeTetPointRef(const FixedList& p) +{ + return tetPointRef(p[0], p[1], p[2], p[3]); +} + + int main(int argc, char *argv[]) { - tetPoints A - ( + // Tets to test + FixedList 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 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> tetList; + tetList tetsIn1, tetsIn2, tetsOut; + cut::appendOp tetOpIn1(tetsIn1); + cut::appendOp tetOpIn2(tetsIn2); + cut::appendOp 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 2c1c3695ad..46140b7ae5 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 | ------------------------------------------------------------------------------- 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 a644134066..8f8aac17e9 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 ef565a405f..c0782e3a3c 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 49f6796781..570ad3e184 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 d2f8d2a7e5..3eff6c526f 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 0a3ed7a03b..18de73777c 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 0c80210e63..8ad320403a 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 0000000000..9e821c0c53 --- /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 . + +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 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 +{ +public: + + // Typedefs + + //- Result type + typedef zero result; + + + // Constructors + + //- Construct null + noOp() + {} + + //- Construct from base + noOp(const uniformOp& op) + : + uniformOp(op) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return result(); + } + + //- Operate on a triangle or tetrahedron + template + result operator()(const FixedList& p) const + { + return result(); + } +}; + + +/*---------------------------------------------------------------------------*\ + Class areaOp Declaration +\*---------------------------------------------------------------------------*/ + +class areaOp +: + public uniformOp +{ +public: + + // Typedefs + + //- Result type + typedef vector result; + + + // Constructors + + //- Construct null + areaOp() + {} + + //- Construct from base + areaOp(const uniformOp& op) + : + uniformOp(op) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return vector::zero; + } + + //- Operate on a triangle + result operator()(const FixedList& p) const + { + return result(triPointRef(p[0], p[1], p[2]).normal()); + } +}; + + +/*---------------------------------------------------------------------------*\ + Class volumeOp Declaration +\*---------------------------------------------------------------------------*/ + +class volumeOp +: + public uniformOp +{ +public: + + // Typedefs + + //- Result type + typedef scalar result; + + + // Constructors + + //- Construct null + volumeOp() + {} + + //- Construct from base + volumeOp(const uniformOp& op) + : + uniformOp(op) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return 0; + } + + //- Operate on a tetrahedron + result operator()(const FixedList& p) const + { + return result(tetPointRef(p[0], p[1], p[2], p[3]).mag()); + } +}; + + +/*---------------------------------------------------------------------------*\ + Class areaIntegrateOp Declaration +\*---------------------------------------------------------------------------*/ + +template +class areaIntegrateOp +: + public FixedList +{ +public: + + // Typedefs + + //- Result type + typedef typename outerProduct::type result; + + + // Constructors + + //- Construct from base + areaIntegrateOp(const FixedList& x) + : + FixedList(x) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return pTraits::zero; + } + + //- Operate on a triangle + result operator()(const FixedList& p) const + { + const FixedList& x = *this; + return result(areaOp()(p)*(x[0] + x[1] + x[2])/3); + } +}; + + +/*---------------------------------------------------------------------------*\ + Class volumeIntegrateOp Declaration +\*---------------------------------------------------------------------------*/ + +template +class volumeIntegrateOp +: + public FixedList +{ +public: + + // Typedefs + + //- Result type + typedef Type result; + + + // Constructors + + //- Construct from base + volumeIntegrateOp(const FixedList& x) + : + FixedList(x) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return pTraits::zero; + } + + //- Operate on a tetrahedron + result operator()(const FixedList& p) const + { + const FixedList& x = *this; + return result(volumeOp()(p)*(x[0] + x[1] + x[2] + x[3])/4); + } +}; + + +/*---------------------------------------------------------------------------*\ + Class listOp Declaration +\*---------------------------------------------------------------------------*/ + +template +class listOp +: + public uniformOp +{ +public: + + // Classes + + //- Result class + class result + : + public DynamicList> + { + public: + + // Constructors + + //- Construct from a single element + result(const FixedList& x) + : + DynamicList>(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& op) + : + uniformOp(op) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return result(); + } + + //- Operate on a triangle or tetrahedron + result operator()(const FixedList& p) const + { + return result(p); + } +}; + + +typedef listOp<3> listTriOp; + + +typedef listOp<4> listTetOp; + + +/*---------------------------------------------------------------------------*\ + Class appendOp Declaration +\*---------------------------------------------------------------------------*/ + +template +class appendOp +: + public uniformOp +{ +public: + + // Typedefs + + //- Result type + typedef zero result; + + + // Constructors + + //- Construct from a container reference + appendOp(Container& x) + : + uniformOp(x) + {} + + //- Construct from base + appendOp(const uniformOp& op) + : + uniformOp(op) + {} + + + // Member operators + + //- Operate on nothing + result operator()() const + { + return result(); + } + + //- Operate on a triangle or tetrahedron + template + result operator()(const FixedList& p) const + { + this->data().append(p); + return result(); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Trait to determine the result of the addition of two operations +template +class opAddResult; + +template +class opAddResult +{ +public: + + typedef typename Op::result type; +}; + +template<> +class opAddResult +{ +public: + + typedef typename noOp::result type; +}; + +template +class opAddResult +{ +public: + + typedef typename Op::result type; +}; + +template +class opAddResult +{ +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 +typename cut::opAddResult::type triCut +( + const FixedList& tri, + const FixedList& level, + const AboveOp& aboveOp, + const BelowOp& belowOp +); + +//- As above, but with a plane specifying the location of the cut +template +typename cut::opAddResult::type triCut +( + const FixedList& tri, + const plane& s, + const AboveOp& aboveOp, + const BelowOp& belowOp +); + +//- As triCut, but for a tetrahedron. +template +typename cut::opAddResult::type tetCut +( + const FixedList& tet, + const FixedList& level, + const AboveOp& aboveOp, + const BelowOp& belowOp +); + +//- As above, but with a plane specifying the location of the cut +template +typename cut::opAddResult::type tetCut +( + const FixedList& 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 0000000000..aefc537f51 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#include "cut.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Modify a uniform operation for reordering a tri (does nothing) +template +inline const cut::uniformOp& triReorder +( + const cut::uniformOp& x, + const FixedList& +) +{ + return x; +} + + +//- Modify a uniform operation for cutting a tri from a tri (does nothing) +template +inline const cut::uniformOp& triCutTri +( + const cut::uniformOp& x, + const Pair& +) +{ + return x; +} + + +//- Modify a uniform operation for cutting a quad from a tri (does nothing) +template +inline const cut::uniformOp& triCutQuad +( + const cut::uniformOp& x, + const Pair& +) +{ + return x; +} + + +//- Modify a uniform operation for reordering a tet (does nothing) +template +inline const cut::uniformOp& tetReorder +( + const cut::uniformOp& x, + const FixedList& +) +{ + return x; +} + + +//- Modify a uniform operation for cutting a tet from a tet (does nothing) +template +inline const cut::uniformOp& tetCutTet +( + const cut::uniformOp& x, + const FixedList& +) +{ + return x; +} + + +//- Modify a uniform operation for cutting prism0 from a tet (does nothing) +template +inline const cut::uniformOp& tetCutPrism0 +( + const cut::uniformOp& x, + const FixedList& +) +{ + return x; +} + + +//- Modify a uniform operation for cutting prism01 from a tet (does nothing) +template +inline const cut::uniformOp& tetCutPrism01 +( + const cut::uniformOp& x, + const FixedList& +) +{ + return x; +} + + +//- Modify a uniform operation for cutting prism23 from a tet (does nothing) +template +inline const cut::uniformOp& tetCutPrism23 +( + const cut::uniformOp& x, + const FixedList& +) +{ + return x; +} + + +//- Modify a fixed list for reordering a tri (does nothing) +template +inline FixedList triReorder +( + const FixedList& x, + const FixedList& indices +) +{ + FixedList 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 +inline FixedList triCutTri +( + const FixedList& x, + const Pair& f +) +{ + FixedList 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 +inline FixedList triCutQuad +( + const FixedList& x, + const Pair& f +) +{ + FixedList 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 +inline FixedList tetReorder +( + const FixedList& x, + const FixedList& indices +) +{ + FixedList 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 +inline FixedList tetCutTet +( + const FixedList& x, + const FixedList& f +) +{ + FixedList 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 +inline FixedList tetCutPrism0 +( + const FixedList& x, + const FixedList& f +) +{ + FixedList 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 +inline FixedList tetCutPrism01 +( + const FixedList& x, + const FixedList& f +) +{ + FixedList 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 +inline FixedList tetCutPrism23 +( + const FixedList& x, + const FixedList& f +) +{ + FixedList 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 +inline typename Op::result triCutTri +( + const Op& op, + const FixedList& p, + const Pair& f +) +{ + return Op(triCutTri(op, f))(triCutTri(p, f)); +} + + +//- Apply an operation to a quad. Splits the quad into two tris. +template +inline typename Op::result quadOp +( + const OpData& opData, + const FixedList& p +) +{ + static const FixedList, 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 +inline typename Op::result triCutQuad +( + const Op& op, + const FixedList& p, + const FixedList& f +) +{ + return quadOp(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 +inline typename Op::result tetCutTet +( + const Op& op, + const FixedList& p, + const FixedList& f +) +{ + return Op(tetCutTet(op, f))(tetCutTet(p, f)); +} + + +//- Apply an operation to a prism. Splits the prism into three tets. +template +inline typename Op::result prismOp +( + const OpData& opData, + const FixedList& p +) +{ + static const FixedList, 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 +inline typename Op::result tetCutPrism0 +( + const Op& op, + const FixedList& p, + const FixedList& f +) +{ + return prismOp(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 +inline typename Op::result tetCutPrism01 +( + const Op& op, + const FixedList& p, + const FixedList& f +) +{ + return prismOp(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 +inline typename Op::result tetCutPrism23 +( + const Op& op, + const FixedList& p, + const FixedList& f +) +{ + return prismOp(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 0000000000..b4dda7f4c0 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#include "cut.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +typename Foam::cut::opAddResult::type Foam::triCut +( + const FixedList& tri, + const FixedList& 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 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 p = triReorder(tri, indices); + const FixedList 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 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 +typename Foam::cut::opAddResult::type Foam::triCut +( + const FixedList& tri, + const plane& p, + const AboveOp& aboveOp, + const BelowOp& belowOp +) +{ + // Set the level set to the signed distance from the plane + FixedList 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 +typename Foam::cut::opAddResult::type Foam::tetCut +( + const FixedList& tet, + const FixedList& 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 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 p = tetReorder(tet, indices); + const FixedList 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 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 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 +typename Foam::cut::opAddResult::type Foam::tetCut +( + const FixedList& tet, + const plane& p, + const AboveOp& aboveOp, + const BelowOp& belowOp +) +{ + // Set the level set to the signed distance from the plane + FixedList 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 0000000000..8854eb46fb --- /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 . + +Typedef + Foam::tetPointRef + +Description + +\*---------------------------------------------------------------------------*/ + +#ifndef tetPointRef_H +#define tetPointRef_H + +#include "point.H" +#include "tetrahedron.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +typedef tetrahedron 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 cb42d33d4b..0000000000 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H +++ /dev/null @@ -1,114 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 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 . - -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), operator[](0)); - for (label i = 1; i < size(); i++) - { - bb.min() = min(bb.min(), operator[](i)); - bb.max() = max(bb.max(), 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 163476a6b3..45b10041f3 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 -void Foam::tetrahedron::tetOverlap -( - const tetrahedron& tetB, - tetIntersectionList& insideTets, - label& nInside, - tetIntersectionList& outsideTets, - label& nOutside -) const -{ - // Work storage - tetIntersectionList cutInsideTets; - label nCutInside = 0; - - nInside = 0; - storeOp inside(insideTets, nInside); - storeOp cutInside(cutInsideTets, nCutInside); - - nOutside = 0; - storeOp outside(outsideTets, nOutside); - - - // Cut tetA with all inwards pointing faces of tetB. Any tets remaining - // in aboveTets are inside tetB. - - { - // face0 - plane pl0(tetB.b_, tetB.d_, tetB.c_); - - // Cut and insert subtets into cutInsideTets (either by getting - // an index from freeSlots or by appending to insideTets) or - // insert into outsideTets - sliceWithPlane(pl0, cutInside, outside); - } - - if (nCutInside == 0) - { - nInside = nCutInside; - return; - } - - { - // face1 - plane pl1(tetB.a_, tetB.c_, tetB.d_); - - nInside = 0; - - for (label i = 0; i < nCutInside; i++) - { - cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside); - } - - if (nInside == 0) - { - return; - } - } - - { - // face2 - plane pl2(tetB.a_, tetB.d_, tetB.b_); - - nCutInside = 0; - - for (label i = 0; i < nInside; i++) - { - insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside); - } - - if (nCutInside == 0) - { - nInside = nCutInside; - return; - } - } - - { - // face3 - plane pl3(tetB.a_, tetB.b_, tetB.c_); - - nInside = 0; - - for (label i = 0; i < nCutInside; i++) - { - cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside); - } - } -} - - template Foam::pointHit Foam::tetrahedron::containmentSphere ( @@ -426,4 +336,16 @@ void Foam::tetrahedron::gradNiGradNj } +template +Foam::boundBox Foam::tetrahedron::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 badea19542..80082ecce1 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 @@ -47,6 +47,7 @@ SourceFiles #include "FixedList.H" #include "UList.H" #include "triPointRef.H" +#include "boundBox.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -55,8 +56,6 @@ namespace Foam class Istream; class Ostream; -class tetPoints; -class plane; // Forward declaration of friend functions and operators @@ -76,8 +75,6 @@ inline Ostream& operator<< const tetrahedron& ); -typedef tetrahedron tetPointRef; - /*---------------------------------------------------------------------------*\ class tetrahedron Declaration \*---------------------------------------------------------------------------*/ @@ -85,78 +82,12 @@ typedef tetrahedron tetPointRef; 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: @@ -259,26 +190,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 - 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 @@ -299,6 +210,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 2c14d16121..08ed1816a9 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 * * * * * * * * * * * * * * // @@ -526,470 +525,6 @@ 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; - 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 -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/OpenFOAM/primitives/zero/zeroI.H b/src/OpenFOAM/primitives/zero/zeroI.H index 0b7b75f237..d69255f159 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 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 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 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 84220fa8f7..52db438141 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 | ------------------------------------------------------------------------------- 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 25081cfd6a..1486d6b0c2 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 db51ee66a6..f1406596a5 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 | ------------------------------------------------------------------------------- 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 * * * * * * * * * * * * * // @@ -50,68 +49,62 @@ Foam::tetOverlapVolume::tetOverlapVolume() Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol ( - const tetPoints& tetA, - const tetPoints& tetB + const tetPointRef& tetA, + const tetPointRef& tetB ) const { - 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::sumVolOp volInside; - tetPointRef::dummyOp outside; - - if ((tetA.tet().mag() < SMALL*SMALL) || (tetB.tet().mag() < SMALL*SMALL)) - { - return 0.0; - } - - // face0 - plane pl0(tetB[1], tetB[3], tetB[2]); - tetA.tet().sliceWithPlane(pl0, cutInside, outside); - if (nCutInside == 0) + // 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.0; + return 0; } - // face1 - plane pl1(tetB[0], tetB[2], tetB[3]); - nInside = 0; - for (label i = 0; i < nCutInside; i++) + // face 1 + const plane pl1(tetB.a(), tetB.c(), tetB.d()); + cutTetList2.clear(); + for (label i = 0; i < cutTetList1.size(); i++) { - const tetPointRef t = cutInsideTets[i].tet(); - t.sliceWithPlane(pl1, inside, outside); + const FixedList& t = cutTetList1[i]; + tetCut(t, pl1, cut::appendOp(cutTetList2), cut::noOp()); } - if (nInside == 0) + if (cutTetList2.size() == 0) { - return 0.0; + return 0; } - // face2 - plane pl2(tetB[0], tetB[3], tetB[1]); - nCutInside = 0; - for (label i = 0; i < nInside; i++) + // face 2 + const plane pl2(tetB.a(), tetB.d(), tetB.b()); + cutTetList1.clear(); + for (label i = 0; i < cutTetList2.size(); i++) { - const tetPointRef t = insideTets[i].tet(); - t.sliceWithPlane(pl2, cutInside, outside); + const FixedList& t = cutTetList2[i]; + tetCut(t, pl2, cut::appendOp(cutTetList1), cut::noOp()); } - if (nCutInside == 0) + if (cutTetList1.size() == 0) { - return 0.0; + return 0; } - // face3 - plane pl3(tetB[0], tetB[1], tetB[2]); - for (label i = 0; i < nCutInside; i++) + // face 3 + const plane pl3(tetB.a(), tetB.b(), tetB.c()); + scalar v = 0; + for (label i = 0; i < cutTetList1.size(); i++) { - const tetPointRef t = cutInsideTets[i].tet(); - t.sliceWithPlane(pl3, volInside, outside); + const FixedList& t = cutTetList1[i]; + v += tetCut(t, pl3, cut::volumeOp(), cut::noOp()); } - return volInside.vol_; + return v; } @@ -189,7 +182,7 @@ bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp pt1I = fA[facePtAI]; } - const tetPoints tetA + const tetPointRef tetA ( ccA, tetBasePtA, @@ -236,7 +229,7 @@ bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp pt1I = fB[facePtBI]; } - const tetPoints tetB + const tetPointRef tetB ( ccB, tetBasePtB, @@ -317,7 +310,7 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp pt1I = fA[facePtAI]; } - const tetPoints tetA + const tetPointRef tetA ( ccA, tetBasePtA, @@ -364,7 +357,7 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp pt1I = fB[facePtBI]; } - const tetPoints tetB + const tetPointRef tetB ( ccB, tetBasePtB, diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H index 87c65dcd89..06b315d47f 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-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -38,13 +38,13 @@ SourceFiles #include "FixedList.H" #include "labelList.H" #include "treeBoundBox.H" +#include "tetPointRef.H" namespace Foam { class primitiveMesh; class polyMesh; -class tetPoints; /*---------------------------------------------------------------------------*\ Class tetOverlapVolume Declaration @@ -57,8 +57,8 @@ class tetOverlapVolume //- Tet overlap volume scalar tetTetOverlapVol ( - const tetPoints& tetA, - const tetPoints& tetB + const tetPointRef& tetA, + const tetPointRef& tetB ) const; //- Return a const treeBoundBox @@ -70,6 +70,51 @@ class tetOverlapVolume ) 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: //- Runtime type information