From ee690f3cae75930f37a5a0de3199cea3e3fb4f5b Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 24 Nov 2015 11:17:49 +0000
Subject: [PATCH] Change default computation of mean plane (and thus of the
parametrization) of plane surfaces, to make it more robust with respect to
perturbations of the boundary.
This is crucial for sensitivity analyses and GFace::relocateMeshVertices(): after perturbation of the
boundary, we want a parametrization of the surface that is "close" to the original one. This cannot be
guaranteed by the automatic SVD-based algorithm.
Unfortunately the new algorithm is also not guaranteed to work in all cases -- but unless we explicitely
prescribe the parametrization of plane surfaces (at least the normal) this is better.
I will revert this as soon as we have something better.
---
Geo/GFace.cpp | 76 +++++++++++++++++++++++++++++++--------------
Numeric/Numeric.cpp | 52 +++++++++++++++++--------------
Numeric/Numeric.h | 7 +++--
3 files changed, 85 insertions(+), 50 deletions(-)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 2b362254c9..aaf2ec1072 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -410,6 +410,57 @@ void GFace::writeGEO(FILE *fp)
void GFace::computeMeanPlane()
{
std::vector<SPoint3> pts;
+
+ if(geomType() == Plane) {
+ // Special case for planar surfaces: we first try to compute the
+ // parametrization in a way that is more robust with respect to
+ // perturbations of the boundary. This is crucial for sensitivity analyses
+ // and GFace::relocateMeshVertices(): after perturbation of the boundary, we
+ // want a parametrization of the surface that is "close" to the original
+ // one. If this fails, we fallback to the classical (SVD-based) algorithm.
+ std::list<GEdge*> edg = edges();
+ for(std::list<GEdge*>::const_iterator ite = edg.begin(); ite != edg.end(); ite++){
+ const GEdge *e = *ite;
+ Range<double> b = e->parBounds(0);
+ GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low()));
+ pts.push_back(SPoint3(p1.x(), p1.y(), p1.z()));
+ GPoint p2 = e->point(b.low() + 0.666 * (b.high() - b.low()));
+ pts.push_back(SPoint3(p2.x(), p2.y(), p2.z()));
+ }
+ bool ok = false;
+ double res[4] = {0., 0., 0., 0.}, xm = 0., ym = 0., zm = 0.;
+ if(pts.size() >= 3){
+ SVector3 d01(pts[0], pts[1]);
+ for(unsigned int i = 2; i < pts.size(); i++){
+ SVector3 d0i(pts[0], pts[i]);
+ SVector3 n = crossprod(d01, d0i);
+ if(norm(n) > 1e-12){
+ res[0] = n.x(); res[1] = n.y(); res[2] = n.z();
+ xm = pts[0].x(); ym = pts[0].y(); zm = pts[0].z();
+ ok = true;
+ break;
+ }
+ }
+ }
+ if(ok){
+ double ex[3], t1[3], t2[3];
+ ex[0] = ex[1] = ex[2] = 0.0;
+ if(res[0] == 0.)
+ ex[0] = 1.0;
+ else if(res[1] == 0.)
+ ex[1] = 1.0;
+ else
+ ex[2] = 1.0;
+ prodve(res, ex, t1);
+ norme(t1);
+ prodve(t1, res, t2);
+ norme(t2);
+ res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
+ fillMeanPlane(res, t1, t2, meanPlane);
+ return;
+ }
+ }
+
std::list<GVertex*> verts = vertices();
std::list<GVertex*>::const_iterator itv = verts.begin();
for(; itv != verts.end(); itv++){
@@ -552,30 +603,7 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
end:
res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
- for(int i = 0; i < 3; i++)
- meanPlane.plan[0][i] = t1[i];
- for(int i = 0; i < 3; i++)
- meanPlane.plan[1][i] = t2[i];
- for(int i = 0; i < 3; i++)
- meanPlane.plan[2][i] = res[i];
-
- meanPlane.a = res[0];
- meanPlane.b = res[1];
- meanPlane.c = res[2];
- meanPlane.d = res[3];
-
- meanPlane.x = meanPlane.y = meanPlane.z = 0.;
- if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
- fabs(meanPlane.a) >= fabs(meanPlane.c) ){
- meanPlane.x = meanPlane.d / meanPlane.a;
- }
- else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
- fabs(meanPlane.b) >= fabs(meanPlane.c)){
- meanPlane.y = meanPlane.d / meanPlane.b;
- }
- else{
- meanPlane.z = meanPlane.d / meanPlane.c;
- }
+ fillMeanPlane(res, t1, t2, meanPlane);
Msg::Debug("Surface: %d", tag());
Msg::Debug("SVD : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 784c023f33..1810e3fd00 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1378,6 +1378,33 @@ int intersection_segments(const SPoint3 &p1, const SPoint3 &p2,
return false;
}
+void fillMeanPlane(double res[4], double t1[3], double t2[3], mean_plane &meanPlane)
+{
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[0][i] = t1[i];
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[1][i] = t2[i];
+ for(int i = 0; i < 3; i++)
+ meanPlane.plan[2][i] = res[i];
+
+ meanPlane.a = res[0];
+ meanPlane.b = res[1];
+ meanPlane.c = res[2];
+ meanPlane.d = res[3];
+
+ meanPlane.x = meanPlane.y = meanPlane.z = 0.;
+ if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
+ fabs(meanPlane.a) >= fabs(meanPlane.c) ){
+ meanPlane.x = meanPlane.d / meanPlane.a;
+ }
+ else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
+ fabs(meanPlane.b) >= fabs(meanPlane.c)){
+ meanPlane.y = meanPlane.d / meanPlane.b;
+ }
+ else{
+ meanPlane.z = meanPlane.d / meanPlane.c;
+ }
+}
void computeMeanPlaneSimple(const std::vector<SPoint3> &points, mean_plane &meanPlane)
{
@@ -1434,30 +1461,7 @@ void computeMeanPlaneSimple(const std::vector<SPoint3> &points, mean_plane &mean
res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
- for(int i = 0; i < 3; i++)
- meanPlane.plan[0][i] = t1[i];
- for(int i = 0; i < 3; i++)
- meanPlane.plan[1][i] = t2[i];
- for(int i = 0; i < 3; i++)
- meanPlane.plan[2][i] = res[i];
-
- meanPlane.a = res[0];
- meanPlane.b = res[1];
- meanPlane.c = res[2];
- meanPlane.d = res[3];
-
- meanPlane.x = meanPlane.y = meanPlane.z = 0.;
- if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
- fabs(meanPlane.a) >= fabs(meanPlane.c) ){
- meanPlane.x = meanPlane.d / meanPlane.a;
- }
- else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
- fabs(meanPlane.b) >= fabs(meanPlane.c)){
- meanPlane.y = meanPlane.d / meanPlane.b;
- }
- else{
- meanPlane.z = meanPlane.d / meanPlane.c;
- }
+ fillMeanPlane(res, t1, t2, meanPlane);
}
void projectPointToPlane(const SPoint3 &pt, SPoint3 &ptProj, const mean_plane &meanPlane)
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 17aad6e177..795ddcb9b8 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -160,10 +160,13 @@ int intersection_segments (const SPoint2 &p1, const SPoint2 &p2,
double x[2]);
//tools for projection onto plane
+void fillMeanPlane(double res[4], double t1[3], double t2[3], mean_plane &meanPlane);
void computeMeanPlaneSimple(const std::vector<SPoint3> &points, mean_plane &meanPlane);
void projectPointToPlane(const SPoint3 &pt, SPoint3 &ptProj, const mean_plane &meanPlane);
-void projectPointsToPlane(const std::vector<SPoint3> &pts, std::vector<SPoint3> &ptsProj, const mean_plane &meanPlane);
-void transformPointsIntoOrthoBasis(const std::vector<SPoint3> &ptsProj, std::vector<SPoint3> &pointsUV,
+void projectPointsToPlane(const std::vector<SPoint3> &pts, std::vector<SPoint3> &ptsProj,
+ const mean_plane &meanPlane);
+void transformPointsIntoOrthoBasis(const std::vector<SPoint3> &ptsProj,
+ std::vector<SPoint3> &pointsUV,
const SPoint3 &ptCG, const mean_plane &meanPlane);
#endif
--
GitLab