From 4ee645493aa195c4bdf29457c02cc21a71b38d95 Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Sun, 3 Dec 2006 00:35:56 +0000 Subject: [PATCH] add extra arg to XYZtoUV so we can use it also for point that are not on the surface --- Geo/GFace.cpp | 47 ++++++++++++++++------------------- Geo/GFace.h | 8 +++--- Mesh/meshGFaceTransfinite.cpp | 7 +++--- 3 files changed, 29 insertions(+), 33 deletions(-) diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp index 7c1d108e5a..6ee192e4ba 100644 --- a/Geo/GFace.cpp +++ b/Geo/GFace.cpp @@ -1,4 +1,4 @@ -// $Id: GFace.cpp,v 1.24 2006-12-03 00:04:31 geuzaine Exp $ +// $Id: GFace.cpp,v 1.25 2006-12-03 00:35:56 geuzaine Exp $ // // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle // @@ -399,37 +399,34 @@ double GFace::curvature (const SPoint2 ¶m) const } void GFace::XYZtoUV(const double X, const double Y, const double Z, - double &U, double &V, - const double relax) const + double &U, double &V, const double relax, + bool onSurface) const { const double Precision = 1.e-8; const int MaxIter = 25; const int NumInitGuess = 11; - double Unew = 0., Vnew = 0., err, err_xyz; - int iter; + double Unew = 0., Vnew = 0.; double mat[3][3], jac[3][3]; - double umin, umax, vmin, vmax; - double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1,0,1}; - + double init[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0.9, 0.1, 0, 1}; Range<double> ru = parBounds(0); Range<double> rv = parBounds(1); - umin = ru.low(); - umax = ru.high(); - vmin = rv.low(); - vmax = rv.high(); + double umin = ru.low(); + double umax = ru.high(); + double vmin = rv.low(); + double vmax = rv.high(); for(int i = 0; i < NumInitGuess; i++){ for(int j = 0; j < NumInitGuess; j++){ U = init[i]; V = init[j]; - err = 1.0; - iter = 1; + int err = 1.0, err_xyz = 0.; + int iter = 1; while(err > Precision && iter < MaxIter) { - GPoint P = point(U,V); - Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V)); + GPoint P = point(U, V); + Pair<SVector3, SVector3> der = firstDer(SPoint2(U, V)); mat[0][0] = der.left().x(); mat[0][1] = der.left().y(); mat[0][2] = der.left().z(); @@ -449,7 +446,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z, jac[2][1] * (Z - P.z())); err = DSQR(Unew - U) + DSQR(Vnew - V); - err_xyz = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()); + if(onSurface) + err_xyz = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()); iter++; U = Unew; V = Vnew; @@ -458,7 +456,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z, if(iter < MaxIter && err <= Precision && Unew <= umax && Vnew <= vmax && Unew >= umin && Vnew >= vmin){ - if(err_xyz > 1.e-5) + if(onSurface && err_xyz > Precision * CTX.lc) Msg(WARNING,"converged for i=%d j=%d (err=%g iter=%d), but err_xyz = %g", i, j, err, iter, err_xyz); return; @@ -466,21 +464,20 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z, } } - if(relax < 1.e-6) + if(relax < 1.e-6){ Msg(GERROR, "Could not converge: surface mesh will be wrong"); + } else { Msg(INFO, "Relaxation factor = %g", 0.75 * relax); - XYZtoUV(X, Y, Z, U, V, 0.75 * relax); + XYZtoUV(X, Y, Z, U, V, 0.75 * relax, onSurface); } } - SPoint2 GFace::parFromPoint(const SPoint3 &p) const { - double U,V; + double U, V; - XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0); + XYZtoUV(p.x(), p.y(), p.z(), U, V, 1.0); - return SPoint2(U,V); + return SPoint2(U, V); } - diff --git a/Geo/GFace.h b/Geo/GFace.h index 94803e19ae..886e55e50d 100644 --- a/Geo/GFace.h +++ b/Geo/GFace.h @@ -76,10 +76,10 @@ class GFace : public GEntity virtual int dim() const {return 2;} virtual void setVisibility(char val, bool recursive=false); - // compute XYZ from parametric UV + // compute the parameters UV from a point XYZ void XYZtoUV(const double X, const double Y, const double Z, - double &U, double &V, - const double relax) const; + double &U, double &V, const double relax, + bool onSurface=true) const; // The bounding box virtual SBoundingBox3d bounds() const; @@ -94,7 +94,7 @@ class GFace : public GEntity // Return the parmater location on the face given a point in space // that is on the face. - virtual SPoint2 parFromPoint(const SPoint3 &) const; + virtual SPoint2 parFromPoint(const SPoint3 &) const; // True if the parameter value is interior to the face. virtual int containsParam(const SPoint2 &pt) const = 0; diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp index 335b9e2a3a..833d9fbb59 100644 --- a/Mesh/meshGFaceTransfinite.cpp +++ b/Mesh/meshGFaceTransfinite.cpp @@ -1,4 +1,4 @@ -// $Id: meshGFaceTransfinite.cpp,v 1.10 2006-12-02 22:48:25 geuzaine Exp $ +// $Id: meshGFaceTransfinite.cpp,v 1.11 2006-12-03 00:35:56 geuzaine Exp $ // // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle // @@ -221,9 +221,8 @@ int MeshTransfiniteSurface( GFace *gf) double zp = TRAN_TRI(m_vertices[iP1]->z(), m_vertices[iP2]->z(), m_vertices[iP3]->z(), m_vertices[N1]->z(), m_vertices[N2]->z(), m_vertices[N3]->z(), u, v); - SPoint2 param = gf->parFromPoint(SPoint3(xp, yp, zp)); - double Up = param.x(); - double Vp = param.y(); + double Up, Vp; + gf->XYZtoUV(xp, yp, zp, Up, Vp, 1.0, false); #endif GPoint gp = gf->point(SPoint2(Up, Vp)); MFaceVertex *newv = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, Up, Vp); -- GitLab