From 03c7ae885d63fc31a3379b02070ceba257073c06 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 14 Oct 2010 16:01:07 +0000
Subject: [PATCH] fix bug (using hack...) due to ambiguous distance function
(same abs value, different sign)
---
Numeric/cartesian.h | 3 ++-
utils/api_demos/mainCartesian.cpp | 20 +++++++++++++++++---
2 files changed, 19 insertions(+), 4 deletions(-)
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index 9fd1cf5cee..b90a934919 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -140,7 +140,8 @@ class cartesianBox {
_childBox = new cartesianBox<scalar>(X, Y, Z, dxi, deta, dzeta,
2 * Nxi, 2 * Neta, 2 * Nzeta,
level - 1);
- }
+ }
+ double getLC(){ return sqrt(_dxi * _dxi + _deta * _deta + _dzeta * _dzeta); }
int getNxi(){ return _Nxi; }
int getNeta(){ return _Neta; }
int getNzeta(){ return _Nzeta; }
diff --git a/utils/api_demos/mainCartesian.cpp b/utils/api_demos/mainCartesian.cpp
index 295fac8950..0187fb1b2b 100644
--- a/utils/api_demos/mainCartesian.cpp
+++ b/utils/api_demos/mainCartesian.cpp
@@ -29,6 +29,8 @@ void insertActiveCells(double x, double y, double z, double rmax,
void computeLevelset(GModel *gm, cartesianBox<double> &box)
{
+ // tolerance for desambiguation
+ const double tol = box.getLC() * 1.e-12;
std::vector<SPoint3> nodes;
std::vector<int> indices;
for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
@@ -60,9 +62,21 @@ void computeLevelset(GModel *gm, cartesianBox<double> &box)
signedDistancesPointsTriangle(localdist, dummy, nodes, P2, P1, P3);
if(dist.empty())
dist = localdist;
- else
- for (unsigned int j = 0; j < localdist.size(); j++)
- dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
+ else{
+ for (unsigned int j = 0; j < localdist.size(); j++){
+ // FIXME: if there is an ambiguity assume we are inside (to
+ // avoid holes in the structure). This is definitely just a
+ // hack, as it could create pockets of matter outside the
+ // structure...
+ if(dist[j] * localdist[j] < 0 &&
+ fabs(fabs(dist[j]) - fabs(localdist[j])) < tol){
+ dist[j] = std::max(dist[j], localdist[j]);
+ }
+ else{
+ dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
+ }
+ }
+ }
}
}
for (unsigned int j = 0; j < dist.size(); j++)
--
GitLab