zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues 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
committed by Andrew Heather
parent 04c11064b3
commit 62e3d37324
23 changed files with 1719 additions and 1268 deletions
Download Patch File Download Diff File Expand all files Collapse all files

329
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 | Copyright (C) 2016-2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@ -24,13 +24,12 @@ License
\*---------------------------------------------------------------------------*/
#include "tetOverlapVolume.H"
#include "tetrahedron.H"
#include "tetPoints.H"
#include "polyMesh.H"
#include "OFstream.H"
#include "treeBoundBox.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "cut.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -48,6 +47,67 @@ Foam::tetOverlapVolume::tetOverlapVolume()
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
(
const tetPointRef& tetA,
const tetPointRef& tetB
) const
{
// A maximum of three cuts are made (the tets that result from the final cut
// are not stored), and each cut can create at most three tets. The
// temporary storage must therefore extend to 3^3 = 27 tets.
typedef cutTetList<27> tetListType;
static tetListType cutTetList1, cutTetList2;
// face 0
const plane pl0(tetB.b(), tetB.d(), tetB.c());
const FixedList<point, 4> t({tetA.a(), tetA.b(), tetA.c(), tetA.d()});
cutTetList1.clear();
tetCut(t, pl0, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
if (cutTetList1.size() == 0)
{
return 0;
}
// face 1
const plane pl1(tetB.a(), tetB.c(), tetB.d());
cutTetList2.clear();
for (label i = 0; i < cutTetList1.size(); i++)
{
const FixedList<point, 4>& t = cutTetList1[i];
tetCut(t, pl1, cut::appendOp<tetListType>(cutTetList2), cut::noOp());
}
if (cutTetList2.size() == 0)
{
return 0;
}
// face 2
const plane pl2(tetB.a(), tetB.d(), tetB.b());
cutTetList1.clear();
for (label i = 0; i < cutTetList2.size(); i++)
{
const FixedList<point, 4>& t = cutTetList2[i];
tetCut(t, pl2, cut::appendOp<tetListType>(cutTetList1), cut::noOp());
}
if (cutTetList1.size() == 0)
{
return 0;
}
// face 3
const plane pl3(tetB.a(), tetB.b(), tetB.c());
scalar v = 0;
for (label i = 0; i < cutTetList1.size(); i++)
{
const FixedList<point, 4>& t = cutTetList1[i];
v += tetCut(t, pl3, cut::volumeOp(), cut::noOp());
}
return v;
}
Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
(
const pointField& points,
@ -74,43 +134,122 @@ bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp
const scalar threshold
) const
{
hasOverlapOp overlapCheckOp(threshold);
cellCellOverlapMinDecomp<hasOverlapOp>
(
meshA,
cellAI,
meshB,
cellBI,
cellBbB,
overlapCheckOp
);
const cell& cFacesA = meshA.cells()[cellAI];
const point& ccA = meshA.cellCentres()[cellAI];
return overlapCheckOp.ok_;
}
const cell& cFacesB = meshB.cells()[cellBI];
const point& ccB = meshB.cellCentres()[cellBI];
scalar vol = 0.0;
forAll(cFacesA, cFA)
{
label faceAI = cFacesA[cFA];
const face& fA = meshA.faces()[faceAI];
const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
if (!pyrA.overlaps(cellBbB))
{
continue;
}
bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
label tetBasePtAI = 0;
const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
{
label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
label otherFacePtAI = fA.fcIndex(facePtAI);
label pt0I = -1;
label pt1I = -1;
if (ownA)
{
pt0I = fA[facePtAI];
pt1I = fA[otherFacePtAI];
}
else
{
pt0I = fA[otherFacePtAI];
pt1I = fA[facePtAI];
}
const tetPointRef tetA
(
ccA,
tetBasePtA,
meshA.points()[pt0I],
meshA.points()[pt1I]
);
const treeBoundBox tetABb(tetA.bounds());
Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
(
const primitiveMesh& meshA,
const label cellAI,
// Loop over tets of cellB
forAll(cFacesB, cFB)
{
label faceBI = cFacesB[cFB];
const primitiveMesh& meshB,
const label cellBI,
const treeBoundBox& cellBbB
) const
{
sumOverlapOp overlapSumOp;
cellCellOverlapMinDecomp<sumOverlapOp>
(
meshA,
cellAI,
meshB,
cellBI,
cellBbB,
overlapSumOp
);
const face& fB = meshB.faces()[faceBI];
const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
if (!pyrB.overlaps(pyrA))
{
continue;
}
return overlapSumOp.iop_.vol_;
bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
label tetBasePtBI = 0;
const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
{
label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
label otherFacePtBI = fB.fcIndex(facePtBI);
label pt0I = -1;
label pt1I = -1;
if (ownB)
{
pt0I = fB[facePtBI];
pt1I = fB[otherFacePtBI];
}
else
{
pt0I = fB[otherFacePtBI];
pt1I = fB[facePtBI];
}
const tetPointRef tetB
(
ccB,
tetBasePtB,
meshB.points()[pt0I],
meshB.points()[pt1I]
);
if (!tetB.bounds().overlaps(tetABb))
{
continue;
}
vol += tetTetOverlapVol(tetA, tetB);
if (vol > threshold)
{
return true;
}
}
}
}
}
return false;
}
@ -125,18 +264,116 @@ Foam::tetOverlapVolume::cellCellOverlapMomentMinDecomp
const treeBoundBox& cellBbB
) const
{
sumOverlapMomentOp overlapSumOp;
cellCellOverlapMinDecomp<sumOverlapMomentOp>
(
meshA,
cellAI,
meshB,
cellBI,
cellBbB,
overlapSumOp
);
const cell& cFacesA = meshA.cells()[cellAI];
const point& ccA = meshA.cellCentres()[cellAI];
return overlapSumOp.iop_.vol_;
const cell& cFacesB = meshB.cells()[cellBI];
const point& ccB = meshB.cellCentres()[cellBI];
scalar vol = 0.0;
forAll(cFacesA, cFA)
{
label faceAI = cFacesA[cFA];
const face& fA = meshA.faces()[faceAI];
const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
if (!pyrA.overlaps(cellBbB))
{
continue;
}
bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
label tetBasePtAI = 0;
const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
{
label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
label otherFacePtAI = fA.fcIndex(facePtAI);
label pt0I = -1;
label pt1I = -1;
if (ownA)
{
pt0I = fA[facePtAI];
pt1I = fA[otherFacePtAI];
}
else
{
pt0I = fA[otherFacePtAI];
pt1I = fA[facePtAI];
}
const tetPointRef tetA
(
ccA,
tetBasePtA,
meshA.points()[pt0I],
meshA.points()[pt1I]
);
const treeBoundBox tetABb(tetA.bounds());
// Loop over tets of cellB
forAll(cFacesB, cFB)
{
label faceBI = cFacesB[cFB];
const face& fB = meshB.faces()[faceBI];
const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
if (!pyrB.overlaps(pyrA))
{
continue;
}
bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
label tetBasePtBI = 0;
const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
{
label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
label otherFacePtBI = fB.fcIndex(facePtBI);
label pt0I = -1;
label pt1I = -1;
if (ownB)
{
pt0I = fB[facePtBI];
pt1I = fB[otherFacePtBI];
}
else
{
pt0I = fB[otherFacePtBI];
pt1I = fB[facePtBI];
}
const tetPointRef tetB
(
ccB,
tetBasePtB,
meshB.points()[pt0I],
meshB.points()[pt1I]
);
if (!tetB.bounds().overlaps(tetABb))
{
continue;
}
vol += tetTetOverlapVol(tetA, tetB);
}
}
}
}
return vol;
}

184
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-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,7 +29,6 @@ Description
SourceFiles
tetOverlapVolume.C
tetOverlapVolumeTemplates.C
\*---------------------------------------------------------------------------*/
@ -39,8 +38,7 @@ SourceFiles
#include "FixedList.H"
#include "labelList.H"
#include "treeBoundBox.H"
#include "Tuple2.H"
#include "tetrahedron.H"
#include "tetPointRef.H"
namespace Foam
{
@ -54,124 +52,67 @@ class polyMesh;
class tetOverlapVolume
{
// Private classes
//- tetPoints handling : sum resulting volumes
class sumMomentOp
{
public:
Tuple2<scalar, point> vol_;
inline sumMomentOp()
:
vol_(0.0, Zero)
{}
inline void operator()(const tetPoints& tet)
{
const tetPointRef t(tet.tet());
scalar tetVol = t.mag();
vol_.first() += tetVol;
vol_.second() += (tetVol*t.centre());
}
};
//- tetPoints combining : check for overlap
class hasOverlapOp
{
public:
const scalar threshold_;
tetPointRef::sumVolOp iop_;
bool ok_;
inline hasOverlapOp(const scalar threshold)
:
threshold_(threshold),
iop_(),
ok_(false)
{}
//- Overlap two tets
inline bool operator()(const tetPoints& A, const tetPoints& B)
{
tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
ok_ = (iop_.vol_ > threshold_);
return ok_;
}
};
//- tetPoints combining : sum overlap volume
class sumOverlapOp
{
public:
tetPointRef::sumVolOp iop_;
inline sumOverlapOp()
:
iop_()
{}
//- Overlap two tets
inline bool operator()(const tetPoints& A, const tetPoints& B)
{
tetTetOverlap<tetPointRef::sumVolOp>(A, B, iop_);
return false;
}
};
//- tetPoints combining : sum overlap volume
class sumOverlapMomentOp
{
public:
sumMomentOp iop_;
inline sumOverlapMomentOp()
:
iop_()
{}
//- Overlap two tets
inline bool operator()(const tetPoints& A, const tetPoints& B)
{
tetTetOverlap<sumMomentOp>(A, B, iop_);
return false;
}
};
// Private member functions
//- Tet overlap calculation
template<class tetPointsOp>
static void tetTetOverlap
//- Tet overlap volume
scalar tetTetOverlapVol
(
const tetPoints& tetA,
const tetPoints& tetB,
tetPointsOp& insideOp
);
//- Cell overlap calculation
template<class tetsOp>
static void cellCellOverlapMinDecomp
(
const primitiveMesh& meshA,
const label cellAI,
const primitiveMesh& meshB,
const label cellBI,
const treeBoundBox& cellBbB,
tetsOp& combineTetsOp
);
const tetPointRef& tetA,
const tetPointRef& tetB
) const;
//- Return a const treeBoundBox
static treeBoundBox pyrBb
treeBoundBox pyrBb
(
const pointField& points,
const face& f,
const point& fc
);
) const;
// Private classes
//- A fixed list of tets which simulates a dynamic list by incrementing
// a counter whenever its append method is called. This is used as an
// optimisation for the tetTetOverlapVol method.
template<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:
@ -218,17 +159,6 @@ public:
const label cellBI,
const treeBoundBox& cellBbB
) const;
//- Calculates the overlap volume and moment
Tuple2<scalar, point> cellCellOverlapMomentMinDecomp
(
const primitiveMesh& meshA,
const label cellAI,
const primitiveMesh& meshB,
const label cellBI,
const treeBoundBox& cellBbB
) const;
};
@ -238,12 +168,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "tetOverlapVolumeTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

259
src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
View File

@ -1,259 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "tetOverlapVolume.H"
#include "primitiveMesh.H"
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
template<class tetPointsOp>
void Foam::tetOverlapVolume::tetTetOverlap
(
const tetPoints& tetA,
const tetPoints& tetB,
tetPointsOp& insideOp
)
{
static tetPointRef::tetIntersectionList insideTets;
label nInside = 0;
static tetPointRef::tetIntersectionList cutInsideTets;
label nCutInside = 0;
tetPointRef::storeOp inside(insideTets, nInside);
tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
tetPointRef::dummyOp outside;
// Precompute the tet face areas and exit early if any face area is
// too small
static FixedList<vector, 4> tetAFaceAreas;
static FixedList<scalar, 4> tetAMag2FaceAreas;
tetPointRef tetATet = tetA.tet();
for (label facei = 0; facei < 4; ++facei)
{
tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
{
return;
}
}
static FixedList<vector, 4> tetBFaceAreas;
static FixedList<scalar, 4> tetBMag2FaceAreas;
tetPointRef tetBTet = tetB.tet();
for (label facei = 0; facei < 4; ++facei)
{
tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
{
return;
}
}
// Face 0
{
vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
plane pl0(tetBTet.tri(0).a(), n, false);
tetA.tet().sliceWithPlane(pl0, cutInside, outside);
if (nCutInside == 0)
{
return;
}
}
// Face 1
{
vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
plane pl1(tetBTet.tri(1).a(), n, false);
nInside = 0;
for (label i = 0; i < nCutInside; i++)
{
const tetPointRef t = cutInsideTets[i].tet();
t.sliceWithPlane(pl1, inside, outside);
}
if (nInside == 0)
{
return;
}
}
// Face 2
{
vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
plane pl2(tetBTet.tri(2).a(), n, false);
nCutInside = 0;
for (label i = 0; i < nInside; i++)
{
const tetPointRef t = insideTets[i].tet();
t.sliceWithPlane(pl2, cutInside, outside);
}
if (nCutInside == 0)
{
return;
}
}
// Face 3
{
vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
plane pl3(tetBTet.tri(3).a(), n, false);
for (label i = 0; i < nCutInside; i++)
{
const tetPointRef t = cutInsideTets[i].tet();
t.sliceWithPlane(pl3, insideOp, outside);
}
}
}
template<class tetsOp>
void Foam::tetOverlapVolume::cellCellOverlapMinDecomp
(
const primitiveMesh& meshA,
const label cellAI,
const primitiveMesh& meshB,
const label cellBI,
const treeBoundBox& cellBbB,
tetsOp& combineTetsOp
)
{
const cell& cFacesA = meshA.cells()[cellAI];
const point& ccA = meshA.cellCentres()[cellAI];
const cell& cFacesB = meshB.cells()[cellBI];
const point& ccB = meshB.cellCentres()[cellBI];
forAll(cFacesA, cFA)
{
label faceAI = cFacesA[cFA];
const face& fA = meshA.faces()[faceAI];
const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
if (!pyrA.overlaps(cellBbB))
{
continue;
}
bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
label tetBasePtAI = 0;
const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
{
label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
label otherFacePtAI = fA.fcIndex(facePtAI);
label pt0I = -1;
label pt1I = -1;
if (ownA)
{
pt0I = fA[facePtAI];
pt1I = fA[otherFacePtAI];
}
else
{
pt0I = fA[otherFacePtAI];
pt1I = fA[facePtAI];
}
const tetPoints tetA
(
ccA,
tetBasePtA,
meshA.points()[pt0I],
meshA.points()[pt1I]
);
const treeBoundBox tetABb(tetA.bounds());
// Loop over tets of cellB
forAll(cFacesB, cFB)
{
label faceBI = cFacesB[cFB];
const face& fB = meshB.faces()[faceBI];
const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
if (!pyrB.overlaps(pyrA))
{
continue;
}
bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
label tetBasePtBI = 0;
const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
{
label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
label otherFacePtBI = fB.fcIndex(facePtBI);
label pt0I = -1;
label pt1I = -1;
if (ownB)
{
pt0I = fB[facePtBI];
pt1I = fB[otherFacePtBI];
}
else
{
pt0I = fB[otherFacePtBI];
pt1I = fB[facePtBI];
}
const tetPoints tetB
(
ccB,
tetBasePtB,
meshB.points()[pt0I],
meshB.points()[pt1I]
);
if (!tetB.bounds().overlaps(tetABb))
{
continue;
}
if (combineTetsOp(tetA, tetB))
{
return;
}
}
}
}
}
}
// ************************************************************************* //