From b34b7fab821d44f6477c65a1752e795f1056c7ae Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 11 Aug 2016 12:14:43 +0100
Subject: [PATCH] ENH: tetOverlapVolume - improved robustness. Fixes #207
---
.../tetOverlapVolume/tetOverlapVolume.C | 6 --
.../tetOverlapVolume/tetOverlapVolume.H | 6 --
.../tetOverlapVolumeTemplates.C | 101 ++++++++++++------
3 files changed, 70 insertions(+), 43 deletions(-)
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index 209c2f7e375..cc54720dd33 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
@@ -39,12 +39,6 @@ namespace Foam
defineTypeNameAndDebug(tetOverlapVolume, 0);
}
-// When to consider a tet to be zero volume. We want to avoid doing clipping
-// against negative volume tets. Tet volume can be calculated incorrectly
-// due to truncation errors. The value below works for single and double
-// precision but could probably be tighter for double precision.
-Foam::scalar Foam::tetOverlapVolume::minTetVolume_ = SMALL*SMALL;
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index b027599218c..a77e608c6e9 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
@@ -54,12 +54,6 @@ class polyMesh;
class tetOverlapVolume
{
- // Private data
-
- //- Minimum tet volume to skip test
- static scalar minTetVolume_;
-
-
// Private classes
//- tetPoints handling : sum resulting volumes
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
index 3ff950311f1..34eded1ff31 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -45,51 +45,90 @@ void Foam::tetOverlapVolume::tetTetOverlap
tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
tetPointRef::dummyOp outside;
- if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_)
+ // 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)
{
- return;
+ tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
+ tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
+ if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
+ {
+ return;
+ }
}
- // face0
- plane pl0(tetB[1], tetB[3], tetB[2]);
- tetA.tet().sliceWithPlane(pl0, cutInside, outside);
- if (nCutInside == 0)
+ static FixedList<vector, 4> tetBFaceAreas;
+ static FixedList<scalar, 4> tetBMag2FaceAreas;
+ tetPointRef tetBTet = tetB.tet();
+ for (label facei = 0; facei < 4; ++facei)
{
- return;
+ tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
+ tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
+ if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
+ {
+ return;
+ }
}
- // face1
- plane pl1(tetB[0], tetB[2], tetB[3]);
- nInside = 0;
- for (label i = 0; i < nCutInside; i++)
- {
- const tetPointRef t = cutInsideTets[i].tet();
- t.sliceWithPlane(pl1, inside, outside);
- }
- if (nInside == 0)
+
+ // Face 0
{
- return;
+ 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;
+ }
}
- // face2
- plane pl2(tetB[0], tetB[3], tetB[1]);
- nCutInside = 0;
- for (label i = 0; i < nInside; i++)
+ // Face 1
{
- const tetPointRef t = insideTets[i].tet();
- t.sliceWithPlane(pl2, cutInside, outside);
+ 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;
+ }
}
- if (nCutInside == 0)
+
+ // Face 2
{
- return;
+ 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;
+ }
}
- // face3
- plane pl3(tetB[0], tetB[1], tetB[2]);
- for (label i = 0; i < nCutInside; i++)
+ // Face 3
{
- const tetPointRef t = cutInsideTets[i].tet();
- t.sliceWithPlane(pl3, insideOp, outside);
+ 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);
+ }
}
}
--
GitLab