diff --git a/applications/test/momentOfInertia/Test-momentOfInertia.C b/applications/test/momentOfInertia/Test-momentOfInertia.C index 0a825df73f..04987c6d1d 100644 --- a/applications/test/momentOfInertia/Test-momentOfInertia.C +++ b/applications/test/momentOfInertia/Test-momentOfInertia.C @@ -35,7 +35,7 @@ Description #include "polyMesh.H" #include "ListOps.H" #include "face.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "triFaceList.H" #include "OFstream.H" #include "meshTools.H" diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C index cd0c02ce64..b3d9d3f515 100644 --- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C +++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H index a237335ccc..d29dde8825 100644 --- a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H +++ b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H @@ -43,7 +43,7 @@ SourceFiles #include "triFace.H" #include "edge.H" #include "pointField.H" -#include "tetPointRef.H" +#include "tetrahedron.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H index 181a2d1122..f0b9be967b 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H @@ -42,7 +42,7 @@ SourceFiles #include "polyMesh.H" #include "coupledPolyPatch.H" #include "syncTools.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "tetIndices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H index 76f82fa78c..8da97e36e9 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H @@ -38,7 +38,7 @@ SourceFiles #define tetIndices_H #include "label.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "triPointRef.H" #include "polyMesh.H" #include "triFace.H" diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index e0998a6ca3..b5d12a042c 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C index f5f189dbc4..1200ed52e8 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C @@ -25,7 +25,7 @@ License #include "primitiveMesh.H" #include "pyramidPointFaceRef.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "ListOps.H" #include "unitConversion.H" #include "SortableList.H" diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H deleted file mode 100644 index 56382a55c7..0000000000 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H +++ /dev/null @@ -1,52 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011 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/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H similarity index 99% rename from src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H rename to src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H index b23f7540dc..2e1ecf288d 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H @@ -35,9 +35,9 @@ SourceFiles #ifndef tetPoints_H #define tetPoints_H +#include "tetrahedron.H" #include "FixedList.H" #include "treeBoundBox.H" -#include "tetPointRef.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index a2ccbc250f..306ba7ebc9 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H @@ -54,6 +54,8 @@ namespace Foam class Istream; class Ostream; +class tetPoints; +class plane; // Forward declaration of friend functions and operators @@ -73,6 +75,7 @@ inline Ostream& operator<< const tetrahedron& ); +typedef tetrahedron tetPointRef; /*---------------------------------------------------------------------------*\ class tetrahedron Declaration @@ -81,10 +84,71 @@ inline Ostream& operator<< template class tetrahedron { +public: + + // 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 + { + FixedList& tets_; + label& nTets_; + + public: + inline storeOp(FixedList&, 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: @@ -184,6 +248,16 @@ 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; + + //- Return (min)containment sphere, i.e. the smallest sphere with // all points inside. Returns pointHit with: // - hit : if sphere is equal to circumsphere diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 40b38f8067..c32418b663 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -25,7 +25,8 @@ License #include "triangle.H" #include "IOstreams.H" -#include "triPointRef.H" +#include "tetPoints.H" +#include "plane.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -492,6 +493,483 @@ 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 +( + FixedList& 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; + label neg0 = -1; + label neg1 = -1; + forAll(d, i) + { + if (d[i] > 0) + { + if (pos0 == -1) + { + pos0 = i; + } + else + { + pos1 = i; + } + } + else + { + if (neg0 == -1) + { + neg0 = i; + } + else + { + neg1 = 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 + { + FatalErrorIn("tetSliceWithPlane(..)") + << "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/functions/DataEntry/TableFile/TableFile.H b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H index f128f5ad51..60106598b1 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H +++ b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index afac722830..5c982cee59 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -26,7 +26,7 @@ License #include "polyMeshGeometry.H" #include "polyMeshTetDecomposition.H" #include "pyramidPointFaceRef.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "syncTools.H" #include "unitConversion.H" diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C index 13828ff7cf..cb449f5648 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C @@ -25,7 +25,6 @@ License #include "cellPointWeight.H" #include "polyMesh.H" -#include "tetPointRef.H" #include "polyMeshTetDecomposition.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H index 56b4759f1a..1832d84302 100644 --- a/src/lagrangian/basic/Cloud/Cloud.H +++ b/src/lagrangian/basic/Cloud/Cloud.H @@ -41,9 +41,6 @@ SourceFiles #include "IOField.H" #include "CompactIOField.H" #include "polyMesh.H" -#include "indexedOctree.H" -#include "treeDataCell.H" -#include "tetPointRef.H" #include "PackedBoolList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 08913b4145..b1938d40c7 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -35,11 +35,10 @@ Description #include "vector.H" #include "Cloud.H" #include "IDLList.H" -#include "labelList.H" #include "pointField.H" #include "faceList.H" #include "OFstream.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "FixedList.H" #include "polyMeshTetDecomposition.H" diff --git a/src/meshTools/momentOfInertia/momentOfInertia.H b/src/meshTools/momentOfInertia/momentOfInertia.H index a4d03ad808..e2fd291781 100644 --- a/src/meshTools/momentOfInertia/momentOfInertia.H +++ b/src/meshTools/momentOfInertia/momentOfInertia.H @@ -37,7 +37,6 @@ SourceFiles #ifndef momentOfInertia_H #define momentOfInertia_H -#include "tetPointRef.H" #include "triFaceList.H" #include "triSurface.H" #include "polyMesh.H" diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C index d36730a0ee..9beaa9952c 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C +++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C @@ -24,436 +24,27 @@ 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" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::tetOverlapVolume, 0); + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::tetOverlapVolume::tetOverlapVolume() {} -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::tetOverlapVolume::~tetOverlapVolume() -{} - - // * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // -Foam::point Foam::tetOverlapVolume::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 -inline void Foam::tetOverlapVolume::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 -inline void Foam::tetOverlapVolume::tetSliceWithPlane -( - const tetPoints& tet, - const plane& pl, - - 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; - label neg0 = -1; - label neg1 = -1; - forAll(d, i) - { - if (d[i] > 0) - { - if (pos0 == -1) - { - pos0 = i; - } - else - { - pos1 = i; - } - } - else - { - if (neg0 == -1) - { - neg0 = i; - } - else - { - neg1 = 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 - { - FatalErrorIn("tetSliceWithPlane(..)") << "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); - } -} - - void Foam::tetOverlapVolume::tetTetOverlap ( const tetPoints& tetA, @@ -462,16 +53,15 @@ void Foam::tetOverlapVolume::tetTetOverlap label& nInside, FixedList& outsideTets, label& nOutside -) +) const { // Work storage FixedList cutInsideTets; label nCutInside = 0; - storeTetOp inside(insideTets, nInside); - storeTetOp cutInside(cutInsideTets, nCutInside); - dummyTetOp outside; - + tetPointRef::storeOp inside(insideTets, nInside); + tetPointRef::storeOp cutInside(cutInsideTets, nCutInside); + tetPointRef::dummyOp outside; // Cut tetA with all inwards pointing faces of tetB. Any tets remaining @@ -484,7 +74,7 @@ void Foam::tetOverlapVolume::tetTetOverlap // Cut and insert subtets into cutInsideTets (either by getting // an index from freeSlots or by appending to insideTets) or // insert into outsideTets - tetSliceWithPlane(tetA, pl0, cutInside, outside); + tetA.tet().sliceWithPlane(pl0, cutInside, outside); } if (nCutInside == 0) @@ -501,7 +91,7 @@ void Foam::tetOverlapVolume::tetTetOverlap for (label i = 0; i < nCutInside; i++) { - tetSliceWithPlane(cutInsideTets[i], pl1, inside, outside); + cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside); } if (nInside == 0) @@ -518,7 +108,7 @@ void Foam::tetOverlapVolume::tetTetOverlap for (label i = 0; i < nInside; i++) { - tetSliceWithPlane(insideTets[i], pl2, cutInside, outside); + insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside); } if (nCutInside == 0) @@ -536,28 +126,27 @@ void Foam::tetOverlapVolume::tetTetOverlap for (label i = 0; i < nCutInside; i++) { - tetSliceWithPlane(cutInsideTets[i], pl3, inside, outside); + cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside); } } } -inline Foam::scalar -Foam::tetOverlapVolume::tetTetOverlapVol +Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol ( const tetPoints& tetA, const tetPoints& tetB -) +) const { FixedList insideTets; label nInside = 0; FixedList cutInsideTets; label nCutInside = 0; - storeTetOp inside(insideTets, nInside); - storeTetOp cutInside(cutInsideTets, nCutInside); - sumTetVolOp volInside; - dummyTetOp outside; + tetPointRef::storeOp inside(insideTets, nInside); + tetPointRef::storeOp cutInside(cutInsideTets, nCutInside); + tetPointRef::sumVolOp volInside; + tetPointRef::dummyOp outside; // face0 plane pl0(tetB[1], tetB[3], tetB[2]); @@ -605,12 +194,12 @@ Foam::tetOverlapVolume::tetTetOverlapVol } -inline const Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb +Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb ( const pointField& points, const face& f, const point& fc -) +) const { treeBoundBox bb(fc, fc); forAll(f, fp) @@ -750,8 +339,8 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp Foam::labelList Foam::tetOverlapVolume::overlappingCells ( - const fvMesh& fromMesh, - const fvMesh& toMesh, + const polyMesh& fromMesh, + const polyMesh& toMesh, const label iTo ) const { @@ -766,30 +355,4 @@ Foam::labelList Foam::tetOverlapVolume::overlappingCells } -/* - - forAll(cellsA, i) - { - label cellAI = cellsA[i]; - treeBoundBox bbA - ( - pointField(meshA.points(), meshA.cellPoints()[cellAI]) - ); - - scalar v = cellCellOverlapVolumeMinDecomp - ( - meshA, - cellAI, - bbA, - meshB, - cellBI, - bbB - ); - - overlapVol += v; - nOverlapTests++; - } - -*/ - // ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H index d0891c14d1..60fe428115 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H +++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H @@ -36,21 +36,17 @@ SourceFiles #ifndef tetOverlapVolume_H #define tetOverlapVolume_H -#include "tetrahedron.H" -#include "fvMesh.H" -#include "plane.H" -#include "tetPointRef.H" -#include "OFstream.H" -#include "meshTools.H" -#include "indexedOctree.H" -#include "treeDataCell.H" -#include "tetPoints.H" -#include "tetCell.H" -#include "EdgeMap.H" +#include "FixedList.H" +#include "labelList.H" +#include "treeBoundBox.H" namespace Foam { +class primitiveMesh; +class polyMesh; +class tetPoints; + /*---------------------------------------------------------------------------*\ Class tetOverlapVolume Declaration \*---------------------------------------------------------------------------*/ @@ -59,82 +55,6 @@ class tetOverlapVolume { // Private member functions - //- Plane intersection - inline point planeIntersection - ( - const FixedList& d, - const tetPoints& t, - const label negI, - const label posI - ); - - - //- Decompose prism - template inline void decomposePrism - ( - const FixedList& points, - TetOp& op - ); - - - //- Helping cľasses - class dummyTetOp - { - public: - - inline void operator()(const tetPoints&){} - }; - - - class sumTetVolOp - { - public: - scalar vol_; - - inline sumTetVolOp() - : - vol_(0.0) - {} - - inline void operator()(const tetPoints& tet) - { - vol_ += tet.tet().mag(); - } - }; - - - class storeTetOp - { - FixedList& tets_; - label& nTets_; - - public: - - inline storeTetOp(FixedList& tets, label& nTets) - : - tets_(tets), - nTets_(nTets) - {} - - inline void operator()(const tetPoints& tet) - { - tets_[nTets_++] = tet; - } - }; - - - //- Slice. Split tet into subtets above and below plane - template - inline void tetSliceWithPlane - ( - const tetPoints& tet, - const plane& pl, - - AboveTetOp& aboveOp, - BelowTetOp& belowOp - ); - - //- Tet overlap void tetTetOverlap ( @@ -144,31 +64,28 @@ class tetOverlapVolume label& nInside, FixedList& outsideTets, label& nOutside - ); - + ) const; //- Tet Overlap Vol - inline scalar tetTetOverlapVol + scalar tetTetOverlapVol ( const tetPoints& tetA, const tetPoints& tetB - ); - + ) const; //- Return a const treeBoundBox - inline const treeBoundBox pyrBb + treeBoundBox pyrBb ( const pointField& points, const face& f, const point& fc - ); - + ) const; public: //- Runtime type information - TypeName("tetOverlapVolume"); + ClassName("tetOverlapVolume"); // Constructors @@ -177,18 +94,14 @@ public: tetOverlapVolume(); - //- Destructor - virtual ~tetOverlapVolume(); - - // Public members //- Return a list of cells in meshA which overlaps with cellBI in // meshB labelList overlappingCells ( - const fvMesh& meshA, - const fvMesh& meshB, + const polyMesh& meshA, + const polyMesh& meshB, const label cellBI ) const; diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C deleted file mode 100644 index d82cf30a50..0000000000 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C +++ /dev/null @@ -1,341 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2012-2012 OpenCFD Ltd. - \\/ 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 - Calculation of shape function product for a tetrahedron - -\*---------------------------------------------------------------------------*/ - -#include "tetrahedron.H" -#include "triPointRef.H" -#include "scalarField.H" - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -// (Probably very inefficient) minimum containment sphere calculation. -// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf: -// Sphere ctr is smallest one of -// - tet circumcentre -// - triangle circumcentre -// - edge mids -template -Foam::pointHit Foam::tetrahedron::containmentSphere -( - const scalar tol -) const -{ - const scalar fac = 1 + tol; - - // Halve order of tolerance for comparisons of sqr. - const scalar facSqr = Foam::sqrt(fac); - - - // 1. Circumcentre itself. - - pointHit pHit(circumCentre()); - pHit.setHit(); - scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_); - - - // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and - // check if 4th point is inside. - - { - point ctr = triPointRef(a_, b_, c_).circumCentre(); - scalar radiusSqr = magSqr(ctr - a_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - { - point ctr = triPointRef(a_, b_, d_).circumCentre(); - scalar radiusSqr = magSqr(ctr - a_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - { - point ctr = triPointRef(a_, c_, d_).circumCentre(); - scalar radiusSqr = magSqr(ctr - a_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - { - point ctr = triPointRef(b_, c_, d_).circumCentre(); - scalar radiusSqr = magSqr(ctr - b_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - - // 3. Try midpoints of edges - - // mid of edge A-B - { - point ctr = 0.5*(a_ + b_); - scalar radiusSqr = magSqr(a_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(c_-ctr) <= testRadSrq - && magSqr(d_-ctr) <= testRadSrq) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge A-C - { - point ctr = 0.5*(a_ + c_); - scalar radiusSqr = magSqr(a_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(b_-ctr) <= testRadSrq - && magSqr(d_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge A-D - { - point ctr = 0.5*(a_ + d_); - scalar radiusSqr = magSqr(a_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(b_-ctr) <= testRadSrq - && magSqr(c_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge B-C - { - point ctr = 0.5*(b_ + c_); - scalar radiusSqr = magSqr(b_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(a_-ctr) <= testRadSrq - && magSqr(d_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge B-D - { - point ctr = 0.5*(b_ + d_); - scalar radiusSqr = magSqr(b_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(a_-ctr) <= testRadSrq - && magSqr(c_-ctr) <= testRadSrq) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge C-D - { - point ctr = 0.5*(c_ + d_); - scalar radiusSqr = magSqr(c_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(a_-ctr) <= testRadSrq - && magSqr(b_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - - pHit.setDistance(sqrt(minRadiusSqr)); - - return pHit; -} - - -template -void Foam::tetrahedron::gradNiSquared -( - scalarField& buffer -) const -{ - // Change of sign between face area vector and gradient - // does not matter because of square - - // Warning: Added a mag to produce positive coefficients even if - // the tetrahedron is twisted inside out. This is pretty - // dangerous, but essential for mesh motion. - scalar magVol = Foam::mag(mag()); - - buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol; - buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol; - buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol; - buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol; -} - - -template -void Foam::tetrahedron::gradNiDotGradNj -( - scalarField& buffer -) const -{ - // Warning. Ordering of edges needs to be the same for a tetrahedron - // class, a tetrahedron cell shape model and a tetCell - - // Warning: Added a mag to produce positive coefficients even if - // the tetrahedron is twisted inside out. This is pretty - // dangerous, but essential for mesh motion. - - // Double change of sign between face area vector and gradient - - scalar magVol = Foam::mag(mag()); - vector sa = Sa(); - vector sb = Sb(); - vector sc = Sc(); - vector sd = Sd(); - - buffer[0] = (1.0/9.0)*(sa & sb)/magVol; - buffer[1] = (1.0/9.0)*(sa & sc)/magVol; - buffer[2] = (1.0/9.0)*(sa & sd)/magVol; - buffer[3] = (1.0/9.0)*(sd & sb)/magVol; - buffer[4] = (1.0/9.0)*(sb & sc)/magVol; - buffer[5] = (1.0/9.0)*(sd & sc)/magVol; -} - - -template -void Foam::tetrahedron::gradNiGradNi -( - tensorField& buffer -) const -{ - // Change of sign between face area vector and gradient - // does not matter because of square - - scalar magVol = Foam::mag(mag()); - - buffer[0] = (1.0/9.0)*sqr(Sa())/magVol; - buffer[1] = (1.0/9.0)*sqr(Sb())/magVol; - buffer[2] = (1.0/9.0)*sqr(Sc())/magVol; - buffer[3] = (1.0/9.0)*sqr(Sd())/magVol; -} - - -template -void Foam::tetrahedron::gradNiGradNj -( - tensorField& buffer -) const -{ - // Warning. Ordering of edges needs to be the same for a tetrahedron - // class, a tetrahedron cell shape model and a tetCell - - // Double change of sign between face area vector and gradient - - scalar magVol = Foam::mag(mag()); - vector sa = Sa(); - vector sb = Sb(); - vector sc = Sc(); - vector sd = Sd(); - - buffer[0] = (1.0/9.0)*(sa * sb)/magVol; - buffer[1] = (1.0/9.0)*(sa * sc)/magVol; - buffer[2] = (1.0/9.0)*(sa * sd)/magVol; - buffer[3] = (1.0/9.0)*(sd * sb)/magVol; - buffer[4] = (1.0/9.0)*(sb * sc)/magVol; - buffer[5] = (1.0/9.0)*(sd * sc)/magVol; -} - - -// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H deleted file mode 100644 index 13e7acd4eb..0000000000 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H +++ /dev/null @@ -1,315 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2012-2012 OpenCFD Ltd. - \\/ 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::tetrahedron - -Description - A tetrahedron primitive. - - Ordering of edges needs to be the same for a tetrahedron - class, a tetrahedron cell shape model and a tetCell. - -SourceFiles - tetrahedronI.H - tetrahedron.C - -\*---------------------------------------------------------------------------*/ - -#ifndef tetrahedron_H -#define tetrahedron_H - -#include "point.H" -#include "primitiveFieldsFwd.H" -#include "pointHit.H" -#include "cachedRandom.H" -#include "Random.H" -#include "FixedList.H" -#include "UList.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -class Istream; -class Ostream; -class tetPoints; -class plane; - -// Forward declaration of friend functions and operators - -template class tetrahedron; - -template -inline Istream& operator>> -( - Istream&, - tetrahedron& -); - -template -inline Ostream& operator<< -( - Ostream&, - const tetrahedron& -); - - -/*---------------------------------------------------------------------------*\ - class tetrahedron Declaration -\*---------------------------------------------------------------------------*/ - -template -class tetrahedron -{ -public: - - // 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 - { - FixedList& tets_; - label& nTets_; - - public: - inline storeOp(FixedList&, 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: - - // Member constants - - enum - { - nVertices = 4, // Number of vertices in tetrahedron - nEdges = 6 // Number of edges in tetrahedron - }; - - - // Constructors - - //- Construct from points - inline tetrahedron - ( - const Point& a, - const Point& b, - const Point& c, - const Point& d - ); - - //- Construct from four points in the list of points - inline tetrahedron - ( - const UList&, - const FixedList& indices - ); - - //- Construct from Istream - inline tetrahedron(Istream&); - - - // Member Functions - - // Access - - //- Return vertices - inline const Point& a() const; - - inline const Point& b() const; - - inline const Point& c() const; - - inline const Point& d() const; - - - // Properties - - //- Return face normal - inline vector Sa() const; - - inline vector Sb() const; - - inline vector Sc() const; - - inline vector Sd() const; - - //- Return centre (centroid) - inline Point centre() const; - - //- Return volume - inline scalar mag() const; - - //- Return circum-centre - inline Point circumCentre() const; - - //- Return circum-radius - inline scalar circumRadius() const; - - //- Return quality: Ratio of tetrahedron and circum-sphere - // volume, scaled so that a regular tetrahedron has a - // quality of 1 - inline scalar quality() const; - - //- Return a random point in the tetrahedron from a - // uniform distribution - inline Point randomPoint(Random& rndGen) const; - - //- Return a random point in the tetrahedron from a - // uniform distribution - inline Point randomPoint(cachedRandom& rndGen) const; - - //- Calculate the barycentric coordinates of the given - // point, in the same order as a, b, c, d. Returns the - // determinant of the solution. - inline scalar barycentric - ( - const point& pt, - List& bary - ) const; - - //- Return nearest point to p on tetrahedron - inline pointHit nearestPoint(const point& p) const; - - //- 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; - - - //- Return (min)containment sphere, i.e. the smallest sphere with - // all points inside. Returns pointHit with: - // - hit : if sphere is equal to circumsphere - // (biggest sphere) - // - point : centre of sphere - // - distance : radius of sphere - // - eligiblemiss: false - // Tol (small compared to 1, e.g. 1E-9) is used to determine - // whether point is inside: mag(pt - ctr) < (1+tol)*radius. - pointHit containmentSphere(const scalar tol) const; - - //- Fill buffer with shape function products - void gradNiSquared(scalarField& buffer) const; - - void gradNiDotGradNj(scalarField& buffer) const; - - void gradNiGradNi(tensorField& buffer) const; - - void gradNiGradNj(tensorField& buffer) const; - - - // IOstream operators - - friend Istream& operator>> - ( - Istream&, - tetrahedron& - ); - - friend Ostream& operator<< - ( - Ostream&, - const tetrahedron& - ); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#include "tetrahedronI.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "tetrahedron.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H deleted file mode 100644 index 922ae0057e..0000000000 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H +++ /dev/null @@ -1,1012 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2012-2012 OpenCFD Ltd. - \\/ 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 "triangle.H" -#include "IOstreams.H" -#include "triPointRef.H" -#include "tetPoints.H" -#include "plane.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -inline Foam::tetrahedron::tetrahedron -( - const Point& a, - const Point& b, - const Point& c, - const Point& d -) -: - a_(a), - b_(b), - c_(c), - d_(d) -{} - - -template -inline Foam::tetrahedron::tetrahedron -( - const UList& points, - const FixedList& indices -) -: - a_(points[indices[0]]), - b_(points[indices[1]]), - c_(points[indices[2]]), - d_(points[indices[3]]) -{} - - -template -inline Foam::tetrahedron::tetrahedron(Istream& is) -{ - is >> *this; -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -inline const Point& Foam::tetrahedron::a() const -{ - return a_; -} - - -template -inline const Point& Foam::tetrahedron::b() const -{ - return b_; -} - - -template -inline const Point& Foam::tetrahedron::c() const -{ - return c_; -} - - -template -inline const Point& Foam::tetrahedron::d() const -{ - return d_; -} - - -template -inline Foam::vector Foam::tetrahedron::Sa() const -{ - return triangle(b_, c_, d_).normal(); -} - - -template -inline Foam::vector Foam::tetrahedron::Sb() const -{ - return triangle(a_, d_, c_).normal(); -} - - -template -inline Foam::vector Foam::tetrahedron::Sc() const -{ - return triangle(a_, b_, d_).normal(); -} - - -template -inline Foam::vector Foam::tetrahedron::Sd() const -{ - return triangle(a_, c_, b_).normal(); -} - - -template -inline Point Foam::tetrahedron::centre() const -{ - return 0.25*(a_ + b_ + c_ + d_); -} - - -template -inline Foam::scalar Foam::tetrahedron::mag() const -{ - return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_)); -} - - -template -inline Point Foam::tetrahedron::circumCentre() const -{ - vector a = b_ - a_; - vector b = c_ - a_; - vector c = d_ - a_; - - scalar lambda = magSqr(c) - (a & c); - scalar mu = magSqr(b) - (a & b); - - vector ba = b ^ a; - vector ca = c ^ a; - - vector num = lambda*ba - mu*ca; - scalar denom = (c & ba); - - if (Foam::mag(denom) < ROOTVSMALL) - { - // Degenerate tetrahedron, returning centre instead of circumCentre. - - return centre(); - } - - return a_ + 0.5*(a + num/denom); -} - - -template -inline Foam::scalar Foam::tetrahedron::circumRadius() const -{ - vector a = b_ - a_; - vector b = c_ - a_; - vector c = d_ - a_; - - scalar lambda = magSqr(c) - (a & c); - scalar mu = magSqr(b) - (a & b); - - vector ba = b ^ a; - vector ca = c ^ a; - - vector num = lambda*ba - mu*ca; - scalar denom = (c & ba); - - if (Foam::mag(denom) < ROOTVSMALL) - { - // Degenerate tetrahedron, returning GREAT for circumRadius. - return GREAT; - } - - return Foam::mag(0.5*(a + num/denom)); -} - - -template -inline Foam::scalar Foam::tetrahedron::quality() const -{ - return - mag() - /( - 8.0/(9.0*sqrt(3.0)) - *pow3(min(circumRadius(), GREAT)) - + ROOTVSMALL - ); -} - - -template -inline Point Foam::tetrahedron::randomPoint -( - Random& rndGen -) const -{ - // Adapted from - // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html - - scalar s = rndGen.scalar01(); - scalar t = rndGen.scalar01(); - scalar u = rndGen.scalar01(); - - if (s + t > 1.0) - { - s = 1.0 - s; - t = 1.0 - t; - } - - if (t + u > 1.0) - { - scalar tmp = u; - u = 1.0 - s - t; - t = 1.0 - tmp; - } - else if (s + t + u > 1.0) - { - scalar tmp = u; - u = s + t + u - 1.0; - s = 1.0 - t - tmp; - } - - return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_; -} - - -template -inline Point Foam::tetrahedron::randomPoint -( - cachedRandom& rndGen -) const -{ - // Adapted from - // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html - - scalar s = rndGen.sample01(); - scalar t = rndGen.sample01(); - scalar u = rndGen.sample01(); - - if (s + t > 1.0) - { - s = 1.0 - s; - t = 1.0 - t; - } - - if (t + u > 1.0) - { - scalar tmp = u; - u = 1.0 - s - t; - t = 1.0 - tmp; - } - else if (s + t + u > 1.0) - { - scalar tmp = u; - u = s + t + u - 1.0; - s = 1.0 - t - tmp; - } - - return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_; -} - - -template -Foam::scalar Foam::tetrahedron::barycentric -( - const point& pt, - List& bary -) const -{ - // From: - // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics) - - vector e0(a_ - d_); - vector e1(b_ - d_); - vector e2(c_ - d_); - - tensor t - ( - e0.x(), e1.x(), e2.x(), - e0.y(), e1.y(), e2.y(), - e0.z(), e1.z(), e2.z() - ); - - scalar detT = det(t); - - if (Foam::mag(detT) < SMALL) - { - // Degenerate tetrahedron, returning 1/4 barycentric coordinates. - - bary = List(4, 0.25); - - return detT; - } - - vector res = inv(t, detT) & (pt - d_); - - bary.setSize(4); - - bary[0] = res.x(); - bary[1] = res.y(); - bary[2] = res.z(); - bary[3] = (1.0 - res.x() - res.y() - res.z()); - - return detT; -} - - -template -inline Foam::pointHit Foam::tetrahedron::nearestPoint -( - const point& p -) const -{ - // Adapted from: - // Real-time collision detection, Christer Ericson, 2005, p142-144 - - // Assuming initially that the point is inside all of the - // halfspaces, so closest to itself. - - point closestPt = p; - - scalar minOutsideDistance = VGREAT; - - bool inside = true; - - if (((p - b_) & Sa()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(b_, c_, d_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - if (((p - a_) & Sb()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(a_, d_, c_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - if (((p - a_) & Sc()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(a_, b_, d_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - if (((p - a_) & Sd()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(a_, c_, b_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - // If the point is inside, then the distance to the closest point - // is zero - if (inside) - { - minOutsideDistance = 0; - } - - return pointHit - ( - inside, - closestPt, - minOutsideDistance, - !inside - ); -} - - -template -bool Foam::tetrahedron::inside(const point& pt) const -{ - // For robustness, assuming that the point is in the tet unless - // "definitively" shown otherwise by obtaining a positive dot - // product greater than a tolerance of SMALL. - - // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal - // vectors and base points for the half-space planes are: - // area[0] = Sa(); - // area[1] = Sb(); - // area[2] = Sc(); - // area[3] = Sd(); - // planeBase[0] = tetBasePt = b_ - // planeBase[1] = ptA = c_ - // planeBase[2] = tetBasePt = b_ - // planeBase[3] = tetBasePt = b_ - - vector n = vector::zero; - - { - // 0, a - const point& basePt = b_; - - n = Sa(); - n /= (Foam::mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - { - // 1, b - const point& basePt = c_; - - n = Sb(); - n /= (Foam::mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - { - // 2, c - const point& basePt = b_; - - n = Sc(); - n /= (Foam::mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - { - // 3, d - const point& basePt = b_; - - n = Sd(); - n /= (Foam::mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - return true; -} - - -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 -( - FixedList& 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; - label neg0 = -1; - label neg1 = -1; - forAll(d, i) - { - if (d[i] > 0) - { - if (pos0 == -1) - { - pos0 = i; - } - else - { - pos1 = i; - } - } - else - { - if (neg0 == -1) - { - neg0 = i; - } - else - { - neg1 = 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 - { - FatalErrorIn("tetSliceWithPlane(..)") - << "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 -inline Foam::Istream& Foam::operator>> -( - Istream& is, - tetrahedron& t -) -{ - is.readBegin("tetrahedron"); - is >> t.a_ >> t.b_ >> t.c_ >> t.d_; - is.readEnd("tetrahedron"); - - is.check("Istream& operator>>(Istream&, tetrahedron&)"); - - return is; -} - - -template -inline Foam::Ostream& Foam::operator<< -( - Ostream& os, - const tetrahedron& t -) -{ - os << nl - << token::BEGIN_LIST - << t.a_ << token::SPACE - << t.b_ << token::SPACE - << t.c_ << token::SPACE - << t.d_ - << token::END_LIST; - - return os; -} - - -// ************************************************************************* //