zhyang-dev/OpenFOAM-5.x
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

primitiveShapes: Generalised tetrahedron and triangle cutting. Cuts are

Browse Source 
now possible with level-sets as well as planes. Removed tetPoints class
as this wasn't really used anywhere except for the old tet-cutting
routines. Restored tetPointRef.H to be consistent with other primitive
shapes. Re-wrote tet-overlap mapping in terms of the new cutting.
This commit is contained in:
Will Bainbridge
2017-05-22 11:40:37 +01:00
parent ac6f881502
commit df1f4be854
22 changed files with 1476 additions and 890 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
applications/test/momentOfInertia/Test-momentOfInertia.C
View File

@ -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"

126
applications/test/tetTetOverlap/Test-tetTetOverlap.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,9 +29,10 @@ Description
\*---------------------------------------------------------------------------*/
#include "tetrahedron.H"
#include "tetPointRef.H"
#include "OFstream.H"
#include "meshTools.H"
#include "cut.H"
using namespace Foam;
@ -41,7 +42,7 @@ void writeOBJ
(
Ostream& os,
label& vertI,
const tetPoints& tet
const FixedList<point, 4>& tet
)
{
forAll(tet, fp)
@ -58,105 +59,126 @@ void writeOBJ
}
tetPointRef makeTetPointRef(const FixedList<point, 4>& p)
{
return tetPointRef(p[0], p[1], p[2], p[3]);
}
int main(int argc, char *argv[])
{
tetPoints A
(
// Tets to test
FixedList<point, 4> tetA
({
point(0, 0, 0),
point(1, 0, 0),
point(1, 1, 0),
point(1, 1, 1)
);
const tetPointRef tetA = A.tet();
tetPoints B
(
});
FixedList<point, 4> tetB
({
point(0.1, 0.1, 0.1),
point(1.1, 0.1, 0.1),
point(1.1, 1.1, 0.1),
point(1.1, 1.1, 1.1)
);
const tetPointRef tetB = B.tet();
});
tetPointRef::tetIntersectionList insideTets;
label nInside = 0;
tetPointRef::tetIntersectionList outsideTets;
label nOutside = 0;
// Do intersection
typedef DynamicList<FixedList<point, 4>> tetList;
tetList tetsIn1, tetsIn2, tetsOut;
cut::appendOp<tetList> tetOpIn1(tetsIn1);
cut::appendOp<tetList> tetOpIn2(tetsIn2);
cut::appendOp<tetList> tetOpOut(tetsOut);
tetA.tetOverlap
(
tetB,
insideTets,
nInside,
outsideTets,
nOutside
);
const plane p0(tetB[1], tetB[3], tetB[2]);
tetsIn1.clear();
tetCut(tetA, p0, tetOpIn1, tetOpOut);
const plane p1(tetB[0], tetB[2], tetB[3]);
tetsIn2.clear();
forAll(tetsIn1, i)
{
tetCut(tetsIn1[i], p1, tetOpIn2, tetOpOut);
}
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)
{
tetCut(tetsIn1[i], p3, tetOpIn2, tetOpOut);
}
const tetList& tetsIn = tetsIn2;
// Dump to file
// ~~~~~~~~~~~~
{
OFstream str("tetA.obj");
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;
}

4
applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
View File

@ -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"

4
applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/fileControl/fileControl.C
View File

@ -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"

4
applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
View File

@ -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"

4
applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellChecks.C
View File

@ -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"

4
src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H
View File

@ -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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
View File

@ -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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
View File

@ -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"

539
src/OpenFOAM/meshes/primitiveShapes/cut/cut.H Normal file
View File

@ -0,0 +1,539 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description:
Functions which cut triangles and tetrahedra. Generic operations are
applied to each half of a cut.
SourceFiles:
cutI.H
cutTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef cut_H
#define cut_H
#include "FixedList.H"
#include "nil.H"
#include "plane.H"
#include "tetPointRef.H"
#include "triPointRef.H"
#include "zero.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace cut
{
/*---------------------------------------------------------------------------*\
Class uniformOp Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class uniformOp
{
private:
//- Data
Type data_;
public:
// Constructors
//- Construct null
uniformOp()
{}
//- Construct from data
uniformOp(Type data)
:
data_(data)
{}
// Member functions
//- Access the data
Type data() const
{
return data_;
}
};
/*---------------------------------------------------------------------------*\
Class noOp Declaration
\*---------------------------------------------------------------------------*/
class noOp
:
public uniformOp<nil>
{
public:
// Typedefs
//- Result type
typedef zero result;
// Constructors
//- Construct null
noOp()
{}
//- Construct from base
noOp(const uniformOp<nil>& op)
:
uniformOp<nil>(op)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return result();
}
//- Operate on a triangle or tetrahedron
template<unsigned Size>
result operator()(const FixedList<point, Size>& p) const
{
return result();
}
};
/*---------------------------------------------------------------------------*\
Class areaOp Declaration
\*---------------------------------------------------------------------------*/
class areaOp
:
public uniformOp<nil>
{
public:
// Typedefs
//- Result type
typedef vector result;
// Constructors
//- Construct null
areaOp()
{}
//- Construct from base
areaOp(const uniformOp<nil>& op)
:
uniformOp<nil>(op)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return vector::zero;
}
//- Operate on a triangle
result operator()(const FixedList<point, 3>& p) const
{
return result(triPointRef(p[0], p[1], p[2]).normal());
}
};
/*---------------------------------------------------------------------------*\
Class volumeOp Declaration
\*---------------------------------------------------------------------------*/
class volumeOp
:
public uniformOp<nil>
{
public:
// Typedefs
//- Result type
typedef scalar result;
// Constructors
//- Construct null
volumeOp()
{}
//- Construct from base
volumeOp(const uniformOp<nil>& op)
:
uniformOp<nil>(op)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return 0;
}
//- Operate on a tetrahedron
result operator()(const FixedList<point, 4>& p) const
{
return result(tetPointRef(p[0], p[1], p[2], p[3]).mag());
}
};
/*---------------------------------------------------------------------------*\
Class areaIntegrateOp Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class areaIntegrateOp
:
public FixedList<Type, 3>
{
public:
// Typedefs
//- Result type
typedef typename outerProduct<Type, vector>::type result;
// Constructors
//- Construct from base
areaIntegrateOp(const FixedList<Type, 3>& x)
:
FixedList<Type, 3>(x)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return pTraits<result>::zero;
}
//- Operate on a triangle
result operator()(const FixedList<point, 3>& p) const
{
const FixedList<Type, 3>& x = *this;
return result(areaOp()(p)*(x[0] + x[1] + x[2])/3);
}
};
/*---------------------------------------------------------------------------*\
Class volumeIntegrateOp Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class volumeIntegrateOp
:
public FixedList<Type, 4>
{
public:
// Typedefs
//- Result type
typedef Type result;
// Constructors
//- Construct from base
volumeIntegrateOp(const FixedList<Type, 4>& x)
:
FixedList<Type, 4>(x)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return pTraits<result>::zero;
}
//- Operate on a tetrahedron
result operator()(const FixedList<point, 4>& p) const
{
const FixedList<Type, 4>& x = *this;
return result(volumeOp()(p)*(x[0] + x[1] + x[2] + x[3])/4);
}
};
/*---------------------------------------------------------------------------*\
Class listOp Declaration
\*---------------------------------------------------------------------------*/
template<unsigned Size>
class listOp
:
public uniformOp<nil>
{
public:
// Classes
//- Result class
class result
:
public DynamicList<FixedList<point, Size>>
{
public:
// Constructors
//- Construct from a single element
result(const FixedList<point, Size>& x)
:
DynamicList<FixedList<point, Size>>(1, x)
{}
// Member operators
//- Add together two lists
result operator+(const result& x) const
{
result r(*this);
r.append(x);
return r;
}
};
// Constructors
//- Construct null
listOp()
{}
//- Construct from base
listOp(const uniformOp<nil>& op)
:
uniformOp<nil>(op)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return result();
}
//- Operate on a triangle or tetrahedron
result operator()(const FixedList<point, Size>& p) const
{
return result(p);
}
};
typedef listOp<3> listTriOp;
typedef listOp<4> listTetOp;
/*---------------------------------------------------------------------------*\
Class appendOp Declaration
\*---------------------------------------------------------------------------*/
template<class Container>
class appendOp
:
public uniformOp<Container&>
{
public:
// Typedefs
//- Result type
typedef zero result;
// Constructors
//- Construct from a container reference
appendOp(Container& x)
:
uniformOp<Container&>(x)
{}
//- Construct from base
appendOp(const uniformOp<Container&>& op)
:
uniformOp<Container&>(op)
{}
// Member operators
//- Operate on nothing
result operator()() const
{
return result();
}
//- Operate on a triangle or tetrahedron
template<unsigned Size>
result operator()(const FixedList<point, Size>& p) const
{
this->data().append(p);
return result();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Trait to determine the result of the addition of two operations
template<class AheadOp, class BehindOp>
class opAddResult;
template<class Op>
class opAddResult<Op, Op>
{
public:
typedef typename Op::result type;
};
template<>
class opAddResult<noOp, noOp>
{
public:
typedef typename noOp::result type;
};
template<class Op>
class opAddResult<noOp, Op>
{
public:
typedef typename Op::result type;
};
template<class Op>
class opAddResult<Op, noOp>
{
public:
typedef typename Op::result type;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace cut
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Cut a triangle along the zero plane defined by the given levels. Templated
// on aboveOp and belowOp; the operations applied to the regions on either side
// of the cut.
template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type triCut
(
const FixedList<point, 3>& tri,
const FixedList<scalar, 3>& level,
const AboveOp& aboveOp,
const BelowOp& belowOp
);
//- As above, but with a plane specifying the location of the cut
template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type triCut
(
const FixedList<point, 3>& tri,
const plane& s,
const AboveOp& aboveOp,
const BelowOp& belowOp
);
//- As triCut, but for a tetrahedron.
template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
(
const FixedList<point, 4>& tet,
const FixedList<scalar, 4>& level,
const AboveOp& aboveOp,
const BelowOp& belowOp
);
//- As above, but with a plane specifying the location of the cut
template<class AboveOp, class BelowOp>
typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
(
const FixedList<point, 4>& tet,
const plane& s,
const AboveOp& aboveOp,
const BelowOp& belowOp
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "cutI.H"
#ifdef NoRepository
#include "cutTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

400
src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H Normal file
View File

@ -0,0 +1,400 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cut.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Modify a uniform operation for reordering a tri (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& triReorder
(
const cut::uniformOp<Type>& x,
const FixedList<label, 3>&
)
{
return x;
}
//- Modify a uniform operation for cutting a tri from a tri (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& triCutTri
(
const cut::uniformOp<Type>& x,
const Pair<scalar>&
)
{
return x;
}
//- Modify a uniform operation for cutting a quad from a tri (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& triCutQuad
(
const cut::uniformOp<Type>& x,
const Pair<scalar>&
)
{
return x;
}
//- Modify a uniform operation for reordering a tet (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& tetReorder
(
const cut::uniformOp<Type>& x,
const FixedList<label, 4>&
)
{
return x;
}
//- Modify a uniform operation for cutting a tet from a tet (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& tetCutTet
(
const cut::uniformOp<Type>& x,
const FixedList<scalar, 3>&
)
{
return x;
}
//- Modify a uniform operation for cutting prism0 from a tet (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& tetCutPrism0
(
const cut::uniformOp<Type>& x,
const FixedList<scalar, 3>&
)
{
return x;
}
//- Modify a uniform operation for cutting prism01 from a tet (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& tetCutPrism01
(
const cut::uniformOp<Type>& x,
const FixedList<scalar, 4>&
)
{
return x;
}
//- Modify a uniform operation for cutting prism23 from a tet (does nothing)
template<class Type>
inline const cut::uniformOp<Type>& tetCutPrism23
(
const cut::uniformOp<Type>& x,
const FixedList<scalar, 4>&
)
{
return x;
}
//- Modify a fixed list for reordering a tri (does nothing)
template<class Type, unsigned Size>
inline FixedList<Type, 3> triReorder
(
const FixedList<Type, Size>& x,
const FixedList<label, 3>& indices
)
{
FixedList<Type, 3> result;
for (unsigned i = 0; i < 3; ++ i)
{
result[i] = x[indices[i]];
}
return result;
}
//- Modify a list for cutting a tri from a tri
template<class Type>
inline FixedList<Type, 3> triCutTri
(
const FixedList<Type, 3>& x,
const Pair<scalar>& f
)
{
FixedList<Type, 3> result;
result[0] = x[0];
for (label i = 0; i < 2; ++ i)
{
result[i+1] = x[0] + f[i]*(x[i+1] - x[0]);
}
return result;
}
//- Modify a list for cutting a quad from a tri
template<class Type>
inline FixedList<Type, 4> triCutQuad
(
const FixedList<Type, 3>& x,
const Pair<scalar>& f
)
{
FixedList<Type, 4> result;
for (label i = 0; i < 2; ++ i)
{
result[i] = x[i+1];
result[3-i] = x[0] + f[i]*(x[i+1] - x[0]);
}
return result;
}
//- Modify a fixed list for reordering a tet (does nothing)
template<class Type, unsigned Size>
inline FixedList<Type, 4> tetReorder
(
const FixedList<Type, Size>& x,
const FixedList<label, 4>& indices
)
{
FixedList<Type, 4> result;
for (unsigned i = 0; i < 4; ++ i)
{
result[i] = x[indices[i]];
}
return result;
}
//- Modify a list for cutting a tet from a tet
template<class Type>
inline FixedList<Type, 4> tetCutTet
(
const FixedList<Type, 4>& x,
const FixedList<scalar, 3>& f
)
{
FixedList<Type, 4> result;
result[0] = x[0];
for (label i = 0; i < 3; ++ i)
{
result[i+1] = x[0] + f[i]*(x[i+1] - x[0]);
}
return result;
}
//- Modify a list for cutting prism0 from a tet
template<class Type>
inline FixedList<Type, 6> tetCutPrism0
(
const FixedList<Type, 4>& x,
const FixedList<scalar, 3>& f
)
{
FixedList<Type, 6> result;
for (label i = 0; i < 3; ++ i)
{
result[i] = x[0] + f[i]*(x[i+1] - x[0]);
result[i+3] = x[i+1];
}
return result;
}
//- Modify a list for cutting prism01 from a tet
template<class Type>
inline FixedList<Type, 6> tetCutPrism01
(
const FixedList<Type, 4>& x,
const FixedList<scalar, 4>& f
)
{
FixedList<Type, 6> result;
for (label i = 0; i < 2; ++ i)
{
result[3*i] = x[i];
for (label j = 0; j < 2; ++ j)
{
result[3*i+j+1] = x[i] + f[2*i+j]*(x[j+2] - x[i]);
}
}
return result;
}
//- Modify a list for cutting prism23 from a tet
template<class Type>
inline FixedList<Type, 6> tetCutPrism23
(
const FixedList<Type, 4>& x,
const FixedList<scalar, 4>& f
)
{
FixedList<Type, 6> result = tetCutPrism01(x, f);
result[0] = x[2];
result[3] = x[3];
Swap(result[2], result[4]);
return result;
}
//- Cut a tri from a tri and apply an operation to the result. The cut is made
// along the two edges connected to vertex 0, and the cut locations are given
// as factors along these edges. The result is the side connected to vertex 0.
template<class Op>
inline typename Op::result triCutTri
(
const Op& op,
const FixedList<point, 3>& p,
const Pair<scalar>& f
)
{
return Op(triCutTri(op, f))(triCutTri(p, f));
}
//- Apply an operation to a quad. Splits the quad into two tris.
template<class Op, class OpData>
inline typename Op::result quadOp
(
const OpData& opData,
const FixedList<point, 4>& p
)
{
static const FixedList<FixedList<label, 3>, 2> i =
{{0, 1, 2}, {0, 2, 3}};
return
Op(triReorder(opData, i[0]))(triReorder(p, i[0]))
+ Op(triReorder(opData, i[1]))(triReorder(p, i[1]));
}
//- Cut a quad from a tri and apply an operation to the result. The cuts are
// the same as for triCutTri. The result is the side connected to vertices 1
// and 2.
template<class Op>
inline typename Op::result triCutQuad
(
const Op& op,
const FixedList<point, 3>& p,
const FixedList<scalar, 2>& f
)
{
return quadOp<Op>(triCutQuad(op, f), triCutQuad(p, f));
}
//- Cut a tet from a tet and apply an operation to the result. The cut is made
// along the three edges connected to vertex 0, and the cut locations are given
// as factors along these edges. The result is the side connected to vertex 0.
template<class Op>
inline typename Op::result tetCutTet
(
const Op& op,
const FixedList<point, 4>& p,
const FixedList<scalar, 3>& f
)
{
return Op(tetCutTet(op, f))(tetCutTet(p, f));
}
//- Apply an operation to a prism. Splits the prism into three tets.
template<class Op, class OpData>
inline typename Op::result prismOp
(
const OpData& opData,
const FixedList<point, 6>& p
)
{
static const FixedList<FixedList<label, 4>, 3> i =
{{0, 1, 2, 4}, {0, 2, 5, 4}, {0, 4, 5, 3}};
return
Op(tetReorder(opData, i[0]))(tetReorder(p, i[0]))
+ Op(tetReorder(opData, i[1]))(tetReorder(p, i[1]))
+ Op(tetReorder(opData, i[2]))(tetReorder(p, i[2]));
}
//- Cut a prism from a tet and apply an operation to the result. The cuts are
// the same as for tetCutTet. The result is the side connected to vertices 1,
// 2 and 3.
template<class Op>
inline typename Op::result tetCutPrism0
(
const Op& op,
const FixedList<point, 4>& p,
const FixedList<scalar, 3>& f
)
{
return prismOp<Op>(tetCutPrism0(op, f), tetCutPrism0(p, f));
}
//- Cut a prism from a tet and apply an operation to the result. The cut is made
// along four edges, not edges 01 or 23, and the cut locations are given as
// factors along these edges. The result is the side connected to edge 01.
template<class Op>
inline typename Op::result tetCutPrism01
(
const Op& op,
const FixedList<point, 4>& p,
const FixedList<scalar, 4>& f
)
{
return prismOp<Op>(tetCutPrism01(op, f), tetCutPrism01(p, f));
}
//- Cut a prism from a tet and apply an operation to the result. The cuts are
// the same as for tetCutPrism01. The result is the side connected to edge 23.
template<class Op>
inline typename Op::result tetCutPrism23
(
const Op& op,
const FixedList<point, 4>& p,
const FixedList<scalar, 4>& f
)
{
return prismOp<Op>(tetCutPrism23(op, f), tetCutPrism23(p, f));
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

263
src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C Normal file
View File

@ -0,0 +1,263 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cut.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class AboveOp, class BelowOp>
typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
(
const FixedList<point, 3>& tri,
const FixedList<scalar, 3>& level,
const AboveOp& aboveOp,
const BelowOp& belowOp
)
{
// If everything is positive or negative, then process the triangle as a
// whole, and do a quick return
if (level[0] >= 0 && level[1] >= 0 && level[2] >= 0)
{
return aboveOp(tri) + belowOp();
}
if (level[0] <= 0 && level[1] <= 0 && level[2] <= 0)
{
return aboveOp() + belowOp(tri);
}
// There will be just one edge without a sign change. Find it, and put it
// opposite the first vertex. This may change the sign of the tri.
FixedList<label, 3> indices({0, 1, 2});
label i;
for (i = 0; i < 3; ++ i)
{
if (level[(i + 1)%3]*level[(i + 2)%3] >= 0)
{
Swap(indices[0], indices[i]);
break;
}
}
if (i == 3)
{
FatalErrorInFunction
<< "The number of tri vertices above the level set should always "
<< "be 1" << exit(FatalError);
}
// Correct the sign
if (indices[0] != 0)
{
Swap(indices[1], indices[2]);
}
// Permute the data
const FixedList<point, 3> p = triReorder(tri, indices);
const FixedList<scalar, 3> l = triReorder(level, indices);
AboveOp a = triReorder(aboveOp, indices);
BelowOp b = triReorder(belowOp, indices);
// Slice off one corner to form a tri and a quad
Pair<scalar> f;
for (label i = 0; i < 2; ++ i)
{
f[i] = l[0]/(l[0] - l[i+1]);
}
if (l[0] > 0)
{
return triCutTri(a, p, f) + triCutQuad(b, p, f);
}
else
{
return triCutQuad(a, p, f) + triCutTri(b, p, f);
}
}
template<class AboveOp, class BelowOp>
typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
(
const FixedList<point, 3>& tri,
const plane& p,
const AboveOp& aboveOp,
const BelowOp& belowOp
)
{
// Set the level set to the signed distance from the plane
FixedList<scalar, 3> level;
for (label i = 0; i < 3; ++ i)
{
level[i] = (tri[i] - p.refPoint()) & p.normal();
}
// Run the level set function
return triCut(tri, level, aboveOp, belowOp);
}
template<class AboveOp, class BelowOp>
typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
(
const FixedList<point, 4>& tet,
const FixedList<scalar, 4>& level,
const AboveOp& aboveOp,
const BelowOp& belowOp
)
{
// Get the min and max over all four vertices and quick return if there is
// no change of sign
scalar levelMin = VGREAT, levelMax = - VGREAT;
for (label i = 0; i < 4; ++ i)
{
levelMin = min(levelMin, level[i]);
levelMax = max(levelMax, level[i]);
}
if (levelMin >= 0)
{
return aboveOp(tet) + belowOp();
}
if (levelMax <= 0)
{
return aboveOp() + belowOp(tet);
}
// Partition the level so that positive values are at the start. This is
// like a single iteration of quick-sort, except that the pivot is a hard-
// coded zero, rather than an element of the array. This can change the sign
// of the tet.
FixedList<label, 4> indices({0, 1, 2, 3});
bool signChange = false;
label i = 0, j = 3;
while (true)
{
while (i < j && level[indices[i]] > 0)
{
i ++;
}
while (j > i && level[indices[j]] <= 0)
{
j --;
}
if (i == j)
{
break;
}
Swap(indices[i], indices[j]);
signChange = !signChange;
}
// The number of vertices above the slice
label n = i;
// If there are more positives than negatives then reverse the order so that
// the negatives are at the start
if (n > 2)
{
n = 4 - n;
for (label i = 0; i < 2; ++ i)
{
Swap(indices[i], indices[3-i]);
}
}
// Correct the sign
if (signChange)
{
Swap(indices[2], indices[3]);
}
// Permute the data
const FixedList<point, 4> p = tetReorder(tet, indices);
const FixedList<scalar, 4> l = tetReorder(level, indices);
AboveOp a = tetReorder(aboveOp, indices);
BelowOp b = tetReorder(belowOp, indices);
// Calculate the integrals above and below the level set
if (n == 1)
{
// Slice off one corner to form a tet and a prism
FixedList<scalar, 3> f;
for (label i = 0; i < 3; ++ i)
{
f[i] = l[0]/(l[0] - l[i+1]);
}
if (l[0] > 0)
{
return tetCutTet(a, p, f) + tetCutPrism0(b, p, f);
}
else
{
return tetCutPrism0(a, p, f) + tetCutTet(b, p, f);
}
}
else if (n == 2)
{
// Slice off two corners to form two prisms
FixedList<scalar, 4> f;
for (label i = 0; i < 2; ++ i)
{
for (label j = 0; j < 2; ++ j)
{
f[2*i+j] = l[i]/(l[i] - l[j+2]);
}
}
if (l[0] > 0)
{
return tetCutPrism01(a, p, f) + tetCutPrism23(b, p, f);
}
else
{
return tetCutPrism23(a, p, f) + tetCutPrism01(b, p, f);
}
}
FatalErrorInFunction
<< "The number of tet vertices above the level set should always be "
<< "either 1 or 2" << exit(FatalError);
return aboveOp() + belowOp();
}
template<class AboveOp, class BelowOp>
typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
(
const FixedList<point, 4>& tet,
const plane& p,
const AboveOp& aboveOp,
const BelowOp& belowOp
)
{
// Set the level set to the signed distance from the plane
FixedList<scalar, 4> level;
for (label i = 0; i < 4; ++ i)
{
level[i] = (tet[i] - p.refPoint()) & p.normal();
}
// Run the level set function
return tetCut(tet, level, aboveOp, belowOp);
}
// ************************************************************************* //

52
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H Normal file
View File

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Typedef
Foam::tetPointRef
Description
\*---------------------------------------------------------------------------*/
#ifndef tetPointRef_H
#define tetPointRef_H
#include "point.H"
#include "tetrahedron.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
typedef tetrahedron<point, const point&> tetPointRef;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#endif
// ************************************************************************* //

114
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H
View File

@ -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 <http://www.gnu.org/licenses/>.
Class
Foam::tetPoints
Description
Tet storage. Null constructable (unfortunately tetrahedron<point, point>
is not)
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef tetPoints_H
#define tetPoints_H
#include "tetrahedron.H"
#include "FixedList.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class tetPoints Declaration
\*---------------------------------------------------------------------------*/
class tetPoints
:
public FixedList<point, 4>
{
public:
// Constructors
//- Construct null
inline tetPoints()
{}
//- Construct from four points
inline tetPoints
(
const point& a,
const point& b,
const point& c,
const point& d
)
{
operator[](0) = a;
operator[](1) = b;
operator[](2) = c;
operator[](3) = d;
}
// Member Functions
//- Return the tetrahedron
inline tetPointRef tet() const
{
return tetPointRef
(
operator[](0),
operator[](1),
operator[](2),
operator[](3)
);
}
//- Calculate the bounding box
inline treeBoundBox bounds() const
{
treeBoundBox bb(operator[](0), 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
// ************************************************************************* //

104
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,96 +29,6 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Point, class PointRef>
void Foam::tetrahedron<Point, PointRef>::tetOverlap
(
const tetrahedron<Point, PointRef>& tetB,
tetIntersectionList& insideTets,
label& nInside,
tetIntersectionList& outsideTets,
label& nOutside
) const
{
// Work storage
tetIntersectionList cutInsideTets;
label nCutInside = 0;
nInside = 0;
storeOp inside(insideTets, nInside);
storeOp cutInside(cutInsideTets, nCutInside);
nOutside = 0;
storeOp outside(outsideTets, nOutside);
// Cut tetA with all inwards pointing faces of tetB. Any tets remaining
// in aboveTets are inside tetB.
{
// face0
plane pl0(tetB.b_, tetB.d_, tetB.c_);
// Cut and insert subtets into cutInsideTets (either by getting
// an index from freeSlots or by appending to insideTets) or
// insert into outsideTets
sliceWithPlane(pl0, cutInside, outside);
}
if (nCutInside == 0)
{
nInside = nCutInside;
return;
}
{
// face1
plane pl1(tetB.a_, tetB.c_, tetB.d_);
nInside = 0;
for (label i = 0; i < nCutInside; i++)
{
cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
}
if (nInside == 0)
{
return;
}
}
{
// face2
plane pl2(tetB.a_, tetB.d_, tetB.b_);
nCutInside = 0;
for (label i = 0; i < nInside; i++)
{
insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
}
if (nCutInside == 0)
{
nInside = nCutInside;
return;
}
}
{
// face3
plane pl3(tetB.a_, tetB.b_, tetB.c_);
nInside = 0;
for (label i = 0; i < nCutInside; i++)
{
cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
}
}
}
template<class Point, class PointRef>
Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
(
@ -426,4 +336,16 @@ void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
}
template<class Point, class PointRef>
Foam::boundBox Foam::tetrahedron<Point, PointRef>::bounds() const
{
return
boundBox
(
min(a(), min(b(), min(c(), d()))),
max(a(), max(b(), max(c(), d())))
);
}
// ************************************************************************* //

96
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
View File

@ -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<Point, PointRef>&
);
typedef tetrahedron<point, const point&> tetPointRef;
/*---------------------------------------------------------------------------*\
class tetrahedron Declaration
\*---------------------------------------------------------------------------*/
@ -85,78 +82,12 @@ typedef tetrahedron<point, const point&> tetPointRef;
template<class Point, class PointRef>
class tetrahedron
{
public:
// Public typedefs
//- Storage type for tets originating from intersecting tets.
// (can possibly be smaller than 200)
typedef FixedList<tetPoints, 200> tetIntersectionList;
// Classes for use in sliceWithPlane. What to do with decomposition
// of tet.
//- Dummy
class dummyOp
{
public:
inline void operator()(const tetPoints&);
};
//- Sum resulting volumes
class sumVolOp
{
public:
scalar vol_;
inline sumVolOp();
inline void operator()(const tetPoints&);
};
//- Store resulting tets
class storeOp
{
tetIntersectionList& tets_;
label& nTets_;
public:
inline storeOp(tetIntersectionList&, label&);
inline void operator()(const tetPoints&);
};
private:
// Private data
PointRef a_, b_, c_, d_;
inline static point planeIntersection
(
const FixedList<scalar, 4>&,
const tetPoints&,
const label,
const label
);
template<class TetOp>
inline static void decomposePrism
(
const FixedList<point, 6>& points,
TetOp& op
);
template<class AboveTetOp, class BelowTetOp>
inline static void tetSliceWithPlane
(
const plane& pl,
const tetPoints& tet,
AboveTetOp& aboveOp,
BelowTetOp& belowOp
);
public:
@ -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<class AboveTetOp, class BelowTetOp>
inline void sliceWithPlane
(
const plane& pl,
AboveTetOp& aboveOp,
BelowTetOp& belowOp
) const;
//- Decompose tet into tets inside and outside other tet
inline void tetOverlap
(
const tetrahedron<Point, PointRef>& tetB,
tetIntersectionList& insideTets,
label& nInside,
tetIntersectionList& outsideTets,
label& nOutside
) const;
//- Return (min)containment sphere, i.e. the smallest sphere with
// all points inside. Returns pointHit with:
// - hit : if sphere is equal to circumsphere
@ -299,6 +210,9 @@ public:
void gradNiGradNj(tensorField& buffer) const;
//- Calculate the bounding box
boundBox bounds() const;
// IOstream operators

467
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
View File

@ -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<Point, PointRef>::inside(const point& pt) const
}
template<class Point, class PointRef>
inline void Foam::tetrahedron<Point, PointRef>::dummyOp::operator()
(
const tetPoints&
)
{}
template<class Point, class PointRef>
inline Foam::tetrahedron<Point, PointRef>::sumVolOp::sumVolOp()
:
vol_(0.0)
{}
template<class Point, class PointRef>
inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
(
const tetPoints& tet
)
{
vol_ += tet.tet().mag();
}
template<class Point, class PointRef>
inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
(
tetIntersectionList& tets,
label& nTets
)
:
tets_(tets),
nTets_(nTets)
{}
template<class Point, class PointRef>
inline void Foam::tetrahedron<Point, PointRef>::storeOp::operator()
(
const tetPoints& tet
)
{
tets_[nTets_++] = tet;
}
template<class Point, class PointRef>
inline Foam::point Foam::tetrahedron<Point, PointRef>::planeIntersection
(
const FixedList<scalar, 4>& d,
const tetPoints& t,
const label negI,
const label posI
)
{
return
(d[posI]*t[negI] - d[negI]*t[posI])
/ (-d[negI]+d[posI]);
}
template<class Point, class PointRef>
template<class TetOp>
inline void Foam::tetrahedron<Point, PointRef>::decomposePrism
(
const FixedList<point, 6>& points,
TetOp& op
)
{
op(tetPoints(points[1], points[3], points[2], points[0]));
op(tetPoints(points[1], points[2], points[3], points[4]));
op(tetPoints(points[4], points[2], points[3], points[5]));
}
template<class Point, class PointRef>
template<class AboveTetOp, class BelowTetOp>
inline void Foam::tetrahedron<Point, PointRef>::
tetSliceWithPlane
(
const plane& pl,
const tetPoints& tet,
AboveTetOp& aboveOp,
BelowTetOp& belowOp
)
{
// Distance to plane
FixedList<scalar, 4> d;
label nPos = 0;
forAll(tet, i)
{
d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
if (d[i] > 0)
{
nPos++;
}
}
if (nPos == 4)
{
aboveOp(tet);
}
else if (nPos == 3)
{
// Sliced into below tet and above prism. Prism gets split into
// two tets.
// Find the below tet
label i0 = -1;
forAll(d, i)
{
if (d[i] <= 0)
{
i0 = i;
break;
}
}
label i1 = d.fcIndex(i0);
label i2 = d.fcIndex(i1);
label i3 = d.fcIndex(i2);
point p01 = planeIntersection(d, tet, i0, i1);
point p02 = planeIntersection(d, tet, i0, i2);
point p03 = planeIntersection(d, tet, i0, i3);
// i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
// ,, 1 : ,, inwards pointing triad
// ,, 2 : ,, outwards pointing triad
// ,, 3 : ,, inwards pointing triad
//Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
if (i0 == 0 || i0 == 2)
{
tetPoints t(tet[i0], p01, p02, p03);
//Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
//checkTet(t, "nPos 3, belowTet i0==0 or 2");
belowOp(t);
// Prism
FixedList<point, 6> p;
p[0] = tet[i1];
p[1] = tet[i3];
p[2] = tet[i2];
p[3] = p01;
p[4] = p03;
p[5] = p02;
//Pout<< " aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
else
{
tetPoints t(p01, p02, p03, tet[i0]);
//Pout<< " belowtet:" << t << " around i0:" << i0 << endl;
//checkTet(t, "nPos 3, belowTet i0==1 or 3");
belowOp(t);
// Prism
FixedList<point, 6> p;
p[0] = tet[i3];
p[1] = tet[i1];
p[2] = tet[i2];
p[3] = p03;
p[4] = p01;
p[5] = p02;
//Pout<< " aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
}
else if (nPos == 2)
{
// Tet cut into two prisms. Determine the positive one.
label pos0 = -1;
label pos1 = -1;
forAll(d, i)
{
if (d[i] > 0)
{
if (pos0 == -1)
{
pos0 = i;
}
else
{
pos1 = i;
}
}
}
//Pout<< "Split 2pos tet " << tet << " d:" << d
// << " around pos0:" << pos0 << " pos1:" << pos1
// << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
const edge posEdge(pos0, pos1);
if (posEdge == edge(0, 1))
{
point p02 = planeIntersection(d, tet, 0, 2);
point p03 = planeIntersection(d, tet, 0, 3);
point p12 = planeIntersection(d, tet, 1, 2);
point p13 = planeIntersection(d, tet, 1, 3);
// Split the resulting prism
{
FixedList<point, 6> p;
p[0] = tet[0];
p[1] = p02;
p[2] = p03;
p[3] = tet[1];
p[4] = p12;
p[5] = p13;
//Pout<< " 01 aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
{
FixedList<point, 6> p;
p[0] = tet[2];
p[1] = p02;
p[2] = p12;
p[3] = tet[3];
p[4] = p03;
p[5] = p13;
//Pout<< " 01 belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else if (posEdge == edge(1, 2))
{
point p01 = planeIntersection(d, tet, 0, 1);
point p13 = planeIntersection(d, tet, 1, 3);
point p02 = planeIntersection(d, tet, 0, 2);
point p23 = planeIntersection(d, tet, 2, 3);
// Split the resulting prism
{
FixedList<point, 6> p;
p[0] = tet[1];
p[1] = p01;
p[2] = p13;
p[3] = tet[2];
p[4] = p02;
p[5] = p23;
//Pout<< " 12 aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
{
FixedList<point, 6> p;
p[0] = tet[3];
p[1] = p23;
p[2] = p13;
p[3] = tet[0];
p[4] = p02;
p[5] = p01;
//Pout<< " 12 belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else if (posEdge == edge(2, 0))
{
point p01 = planeIntersection(d, tet, 0, 1);
point p03 = planeIntersection(d, tet, 0, 3);
point p12 = planeIntersection(d, tet, 1, 2);
point p23 = planeIntersection(d, tet, 2, 3);
// Split the resulting prism
{
FixedList<point, 6> p;
p[0] = tet[2];
p[1] = p12;
p[2] = p23;
p[3] = tet[0];
p[4] = p01;
p[5] = p03;
//Pout<< " 20 aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
{
FixedList<point, 6> p;
p[0] = tet[1];
p[1] = p12;
p[2] = p01;
p[3] = tet[3];
p[4] = p23;
p[5] = p03;
//Pout<< " 20 belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else if (posEdge == edge(0, 3))
{
point p01 = planeIntersection(d, tet, 0, 1);
point p02 = planeIntersection(d, tet, 0, 2);
point p13 = planeIntersection(d, tet, 1, 3);
point p23 = planeIntersection(d, tet, 2, 3);
// Split the resulting prism
{
FixedList<point, 6> p;
p[0] = tet[3];
p[1] = p23;
p[2] = p13;
p[3] = tet[0];
p[4] = p02;
p[5] = p01;
//Pout<< " 03 aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
{
FixedList<point, 6> p;
p[0] = tet[2];
p[1] = p23;
p[2] = p02;
p[3] = tet[1];
p[4] = p13;
p[5] = p01;
//Pout<< " 03 belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else if (posEdge == edge(1, 3))
{
point p01 = planeIntersection(d, tet, 0, 1);
point p12 = planeIntersection(d, tet, 1, 2);
point p03 = planeIntersection(d, tet, 0, 3);
point p23 = planeIntersection(d, tet, 2, 3);
// Split the resulting prism
{
FixedList<point, 6> p;
p[0] = tet[1];
p[1] = p12;
p[2] = p01;
p[3] = tet[3];
p[4] = p23;
p[5] = p03;
//Pout<< " 13 aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
{
FixedList<point, 6> p;
p[0] = tet[2];
p[1] = p12;
p[2] = p23;
p[3] = tet[0];
p[4] = p01;
p[5] = p03;
//Pout<< " 13 belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else if (posEdge == edge(2, 3))
{
point p02 = planeIntersection(d, tet, 0, 2);
point p12 = planeIntersection(d, tet, 1, 2);
point p03 = planeIntersection(d, tet, 0, 3);
point p13 = planeIntersection(d, tet, 1, 3);
// Split the resulting prism
{
FixedList<point, 6> p;
p[0] = tet[2];
p[1] = p02;
p[2] = p12;
p[3] = tet[3];
p[4] = p03;
p[5] = p13;
//Pout<< " 23 aboveprism:" << p << endl;
decomposePrism(p, aboveOp);
}
{
FixedList<point, 6> p;
p[0] = tet[0];
p[1] = p02;
p[2] = p03;
p[3] = tet[1];
p[4] = p12;
p[5] = p13;
//Pout<< " 23 belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else
{
FatalErrorInFunction
<< "Missed edge:" << posEdge
<< abort(FatalError);
}
}
else if (nPos == 1)
{
// Find the positive tet
label i0 = -1;
forAll(d, i)
{
if (d[i] > 0)
{
i0 = i;
break;
}
}
label i1 = d.fcIndex(i0);
label i2 = d.fcIndex(i1);
label i3 = d.fcIndex(i2);
point p01 = planeIntersection(d, tet, i0, i1);
point p02 = planeIntersection(d, tet, i0, i2);
point p03 = planeIntersection(d, tet, i0, i3);
//Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
if (i0 == 0 || i0 == 2)
{
tetPoints t(tet[i0], p01, p02, p03);
//Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
//checkTet(t, "nPos 1, aboveTets i0==0 or 2");
aboveOp(t);
// Prism
FixedList<point, 6> p;
p[0] = tet[i1];
p[1] = tet[i3];
p[2] = tet[i2];
p[3] = p01;
p[4] = p03;
p[5] = p02;
//Pout<< " belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
else
{
tetPoints t(p01, p02, p03, tet[i0]);
//Pout<< " abovetet:" << t << " around i0:" << i0 << endl;
//checkTet(t, "nPos 1, aboveTets i0==1 or 3");
aboveOp(t);
// Prism
FixedList<point, 6> p;
p[0] = tet[i3];
p[1] = tet[i1];
p[2] = tet[i2];
p[3] = p03;
p[4] = p01;
p[5] = p02;
//Pout<< " belowprism:" << p << endl;
decomposePrism(p, belowOp);
}
}
else // nPos == 0
{
belowOp(tet);
}
}
template<class Point, class PointRef>
template<class AboveTetOp, class BelowTetOp>
inline void Foam::tetrahedron<Point, PointRef>::sliceWithPlane
(
const plane& pl,
AboveTetOp& aboveOp,
BelowTetOp& belowOp
) const
{
tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp);
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class Point, class PointRef>

17
src/OpenFOAM/primitives/zero/zeroI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -32,6 +32,11 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline zero operator+(const zero&, const zero&)
{
return Zero;
}
template<class Type>
inline const Type& operator+(const Type& t, const zero&)
{
@ -44,6 +49,11 @@ inline const Type& operator+(const zero&, const Type& t)
return t;
}
inline zero operator-(const zero&, const zero&)
{
return Zero;
}
template<class Type>
inline const Type& operator-(const Type& t, const zero&)
{
@ -56,6 +66,11 @@ inline Type operator-(const zero&, const Type& t)
return -t;
}
inline zero operator*(const zero&, const zero&)
{
return Zero;
}
template<class Type>
inline zero operator*(const Type& t, const zero&)
{

4
src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
View File

@ -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"

2
src/lagrangian/basic/particle/particle.H
View File

@ -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"

97
src/meshTools/tetOverlapVolume/tetOverlapVolume.C
View File

@ -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;
// 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;
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))
// face 0
const plane pl0(tetB.b(), tetB.d(), tetB.c());
const FixedList<point, 4> t({tetA.a(), tetA.b(), tetA.c(), tetA.d()});
cutTetList1.clear();
tetCut(t, pl0, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
if (cutTetList1.size() == 0)
{
return 0.0;
return 0;
}
// face0
plane pl0(tetB[1], tetB[3], tetB[2]);
tetA.tet().sliceWithPlane(pl0, cutInside, outside);
if (nCutInside == 0)
// face 1
const plane pl1(tetB.a(), tetB.c(), tetB.d());
cutTetList2.clear();
for (label i = 0; i < cutTetList1.size(); i++)
{
return 0.0;
const FixedList<point, 4>& t = cutTetList1[i];
tetCut(t, pl1, cut::appendOp<tetListType>(cutTetList2), cut::noOp());
}
if (cutTetList2.size() == 0)
{
return 0;
}
// face1
plane pl1(tetB[0], tetB[2], tetB[3]);
nInside = 0;
for (label i = 0; i < nCutInside; 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 = cutInsideTets[i].tet();
t.sliceWithPlane(pl1, inside, outside);
const FixedList<point, 4>& t = cutTetList2[i];
tetCut(t, pl2, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
}
if (nInside == 0)
if (cutTetList1.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 3
const plane pl3(tetB.a(), tetB.b(), tetB.c());
scalar v = 0;
for (label i = 0; i < cutTetList1.size(); i++)
{
const tetPointRef t = insideTets[i].tet();
t.sliceWithPlane(pl2, cutInside, outside);
}
if (nCutInside == 0)
{
return 0.0;
const FixedList<point, 4>& t = cutTetList1[i];
v += tetCut(t, pl3, cut::volumeOp(), cut::noOp());
}
// face3
plane pl3(tetB[0], tetB[1], tetB[2]);
for (label i = 0; i < nCutInside; i++)
{
const tetPointRef t = cutInsideTets[i].tet();
t.sliceWithPlane(pl3, volInside, outside);
}
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,

53
src/meshTools/tetOverlapVolume/tetOverlapVolume.H
View File

@ -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<unsigned Size>
class cutTetList
:
public FixedList<FixedList<point, 4>, 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<point, 4>& t)
{
this->operator[](n_) = t;
++ n_;
}
};
public:
//- Runtime type information