diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 6ee192e4bad1d6ca1219ece0742013076568a234..a6476f5467f2f14bbeec71adc50a5779433ebb97 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.25 2006-12-03 00:35:56 geuzaine Exp $
+// $Id: GFace.cpp,v 1.26 2006-12-03 01:09:34 geuzaine Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -399,34 +399,37 @@ 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,
- bool onSurface) const
+ double &U, double &V, const double relax,
+ const bool onSurface) const
{
const double Precision = 1.e-8;
const int MaxIter = 25;
const int NumInitGuess = 11;
- double Unew = 0., Vnew = 0.;
+ double Unew = 0., Vnew = 0., err, err2;
+ int iter;
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};
+
Range<double> ru = parBounds(0);
Range<double> rv = parBounds(1);
- double umin = ru.low();
- double umax = ru.high();
- double vmin = rv.low();
- double vmax = rv.high();
+ umin = ru.low();
+ umax = ru.high();
+ vmin = rv.low();
+ vmax = rv.high();
for(int i = 0; i < NumInitGuess; i++){
for(int j = 0; j < NumInitGuess; j++){
U = init[i];
V = init[j];
- int err = 1.0, err_xyz = 0.;
- int iter = 1;
+ err = 1.0;
+ 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();
@@ -446,38 +449,39 @@ 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);
- if(onSurface)
- err_xyz = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z());
+ err2 = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z());
iter++;
U = Unew;
V = Vnew;
}
- if(iter < MaxIter && err <= Precision &&
+
+ if(iter < MaxIter && err <= Precision &&
Unew <= umax && Vnew <= vmax &&
Unew >= umin && Vnew >= vmin){
- 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);
+ if (onSurface && err2 > Precision * CTX.lc)
+ Msg(WARNING,"Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g",
+ i, j, err, iter, err2);
return;
}
}
}
- 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, onSurface);
+ XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
}
}
+
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 886e55e50d5a3b206b5ace865173bd9189062799..ed4f9b4d08916fb9f79bb35c1cd92224098d9285 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -78,8 +78,8 @@ class GFace : public GEntity
// 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,
- bool onSurface=true) const;
+ double &U, double &V, const double relax,
+ const bool onSurface=true) const;
// The bounding box
virtual SBoundingBox3d bounds() const;
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 001879993d589696a916f68fa1554f51716500be..8e9ec7abe2be284604c3f16c18a4e30518839771 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.153 2006-12-02 19:29:36 geuzaine Exp $
+# $Id: Makefile,v 1.154 2006-12-03 01:09:34 geuzaine Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -195,17 +195,18 @@ meshGRegionDelaunayInsertion.o: meshGRegionDelaunayInsertion.cpp \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/MElement.h \
../Geo/ExtrudeParams.h ../Common/Message.h
meshGRegionTransfinite.o: meshGRegionTransfinite.cpp meshGFace.h \
- ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h \
- ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h \
+ ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/Range.h \
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+ ../Common/GmshDefines.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
+ ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/MVertex.h \
+ ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h \
../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
- ../Geo/ExtrudeParams.h ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h \
- ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h \
- ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h \
+ ../Geo/ExtrudeParams.h ../Geo/MElement.h ../Geo/SPoint2.h \
+ ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h \
+ ../Geo/GEntity.h ../Geo/MElement.h ../Geo/ExtrudeParams.h \
../Common/Message.h
meshGRegionExtruded.o: meshGRegionExtruded.cpp ../Geo/ExtrudeParams.h \
../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index 833d9fbb5979e52b6392997cd7f5bbeb11947ed4..782689039ca0e5d97e773e710c65b6e167c1a4b4 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceTransfinite.cpp,v 1.11 2006-12-03 00:35:56 geuzaine Exp $
+// $Id: meshGFaceTransfinite.cpp,v 1.12 2006-12-03 01:09:34 geuzaine Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -208,7 +208,7 @@ int MeshTransfiniteSurface( GFace *gf)
int iP1 = N1 + i;
int iP2 = N2 + j;
int iP3 = ((N3 + N2) - i) % m_vertices.size();
-#if 0 // FIXME: this is buggy, so let's just do it in real space instead
+#if 0 // FIXME: this is buggy, so let's just do it in real space instead for now
double Up = TRAN_TRI(U[iP1], U[iP2], U[iP3], UC1, UC2, UC3, u, v);
double Vp = TRAN_TRI(V[iP1], V[iP2], V[iP3], VC1, VC2, VC3, u, v);
#else
@@ -222,7 +222,14 @@ int MeshTransfiniteSurface( GFace *gf)
m_vertices[iP3]->z(), m_vertices[N1]->z(),
m_vertices[N2]->z(), m_vertices[N3]->z(), u, v);
double Up, Vp;
- gf->XYZtoUV(xp, yp, zp, Up, Vp, 1.0, false);
+ if(gf->geomType() == GEntity::Plane){
+ SPoint2 param = gf->parFromPoint(SPoint3(xp, yp, zp));
+ Up = param.x();
+ Vp = param.y();
+ }
+ else{ // xp, yp, zp is usually not on the surface
+ 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);