diff --git a/Mesh/2D_Mesh.cpp b/Mesh/2D_Mesh.cpp
index 9ef24919d079f1290c411b3ea624d26751faeece..53e9aa51e4b7fc1ac84d20cc9206fe16f42f230a 100644
--- a/Mesh/2D_Mesh.cpp
+++ b/Mesh/2D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh.cpp,v 1.39 2002-02-01 14:34:05 remacle Exp $
+// $Id: 2D_Mesh.cpp,v 1.40 2002-02-12 20:11:34 geuzaine Exp $
/*
Maillage Delaunay d'une surface (Point insertion Technique)
@@ -145,7 +145,7 @@ void Calcule_Z (void *data, void *dum){
if (v->Frozen || v->Num < 0)
return;
- XYtoUV (THESUPPORT, &v->Pos.X, &v->Pos.Y, &U, &V, &Z);
+ XYtoUV (THESUPPORT, &v->Pos.X, &v->Pos.Y, &U, &V, &Z, 1.0);
v->Pos.Z = Z;
}
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index a00df881c4490c31332102b25ebe88eeccc7298e..deac822e13c222f1b802a58937542af5de279656 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: SecondOrder.cpp,v 1.7 2001-10-29 08:52:20 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.8 2002-02-12 20:11:34 geuzaine Exp $
#include "Gmsh.h"
#include "Geo.h"
@@ -45,8 +45,8 @@ Vertex *middleface (Vertex * v1, Vertex * v2){
if (THES->Typ == MSH_SURF_PLAN)
return NULL;
- XYZtoUV ( THES , v1->Pos.X , v1->Pos.Y , v1->Pos.Z, &U1 , &V1 );
- XYZtoUV ( THES , v2->Pos.X , v2->Pos.Y , v2->Pos.Z, &U2 , &V2 );
+ XYZtoUV ( THES , v1->Pos.X , v1->Pos.Y , v1->Pos.Z, &U1 , &V1 , 1.0);
+ XYZtoUV ( THES , v2->Pos.X , v2->Pos.Y , v2->Pos.Z, &U2 , &V2 , 1.0);
U = 0.5 *(U1+U2);
V = 0.5 *(V1+V2);
diff --git a/Mesh/Utils.cpp b/Mesh/Utils.cpp
index 9799e73d40c5ce5c8be56c65d73a35eb750f8f06..6d7d3b2b6ba4c13bbecb10be33b7df605c450873 100644
--- a/Mesh/Utils.cpp
+++ b/Mesh/Utils.cpp
@@ -1,4 +1,4 @@
-// $Id: Utils.cpp,v 1.8 2001-12-16 05:16:37 remacle Exp $
+// $Id: Utils.cpp,v 1.9 2002-02-12 20:11:34 geuzaine Exp $
#include "Gmsh.h"
#include "Numeric.h"
@@ -178,7 +178,7 @@ end:
#define Precision 1.e-10
-#define MaxIter 20
+#define MaxIter 50
void find_bestuv (Surface * s, double X, double Y,
double *U, double *V, double *Z, int N){
@@ -253,11 +253,24 @@ void invert_singular_matrix(double **M, int n, double **I){
free_dvector(W,1,n);
}
-void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
+void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
+ double relax) {
double Unew,Vnew,err;
int iter;
Vertex D_u,D_v,P;
double **mat, **jac ;
+ double umin, umax, vmin, vmax;
+
+ if (s->Typ == MSH_SURF_NURBS){
+ umin = s->ku[0];
+ umax = s->ku[s->OrderU + s->Nu];
+ vmin = s->kv[0];
+ vmax = s->kv[s->OrderV + s->Nv];
+ }
+ else{
+ umin = vmin = 0.0;
+ umax = vmax = 1.0;
+ }
mat = dmatrix(1,3,1,3);
jac = dmatrix(1,3,1,3);
@@ -282,8 +295,10 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
mat[3][3] = 0.;
invert_singular_matrix(mat,3,jac);
- Unew = *U + jac[1][1] * (X-P.Pos.X) + jac[2][1] * (Y-P.Pos.Y) + jac[3][1] * (Z-P.Pos.Z) ;
- Vnew = *V + jac[1][2] * (X-P.Pos.X) + jac[2][2] * (Y-P.Pos.Y) + jac[3][2] * (Z-P.Pos.Z) ;
+ Unew = *U + relax * (jac[1][1] * (X-P.Pos.X) + jac[2][1] * (Y-P.Pos.Y) +
+ jac[3][1] * (Z-P.Pos.Z)) ;
+ Vnew = *V + relax * (jac[1][2] * (X-P.Pos.X) + jac[2][2] * (Y-P.Pos.Y) +
+ jac[3][2] * (Z-P.Pos.Z)) ;
err = DSQR(Unew - *U) + DSQR(Vnew - *V) ;
@@ -292,18 +307,23 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V) {
*V = Vnew;
}
- if(iter > 10){
- if(iter == MaxIter) Msg(WARNING, "Could not converge in XYZtoUV");
- else Msg(WARNING, "Many (%d) iterations in XYZtoUV", iter);
- }
-
free_dmatrix(mat,1,3,1,3);
free_dmatrix(jac,1,3,1,3);
+ if(iter == MaxIter ||
+ (Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin)){
+ if(relax < 1.e-12)
+ Msg(GERROR, "Could not converge: surface mesh will be wrong");
+ else{
+ Msg(INFO, "Relaxation factor = %g", 0.75*relax);
+ XYZtoUV (s, X, Y, Z, U, V, 0.75*relax);
+ }
+ }
+
}
void XYtoUV (Surface * s, double *X, double *Y,
- double *U, double *V, double *Z){
+ double *U, double *V, double *Z, double relax){
double det, Unew, Vnew, err, mat[2][2], jac[2][2];
int iter;
@@ -345,8 +365,8 @@ void XYtoUV (Surface * s, double *X, double *Y,
jac[1][0] = -mat[1][0] / det;
jac[1][1] = mat[0][0] / det;
- Unew = *U + 1.0 * (jac[0][0] * (*X - P.Pos.X) + jac[1][0] * (*Y - P.Pos.Y));
- Vnew = *V + 1.0 * (jac[0][1] * (*X - P.Pos.X) + jac[1][1] * (*Y - P.Pos.Y));
+ Unew = *U + relax * (jac[0][0] * (*X - P.Pos.X) + jac[1][0] * (*Y - P.Pos.Y));
+ Vnew = *V + relax * (jac[0][1] * (*X - P.Pos.X) + jac[1][1] * (*Y - P.Pos.Y));
err = DSQR (Unew - *U) + DSQR (Vnew - *V);
@@ -357,33 +377,30 @@ void XYtoUV (Surface * s, double *X, double *Y,
*Z = P.Pos.Z;
- if(iter > 10){
- if(iter == MaxIter) Msg(WARNING, "Could not converge in XYtoUV");
- else Msg(WARNING, "Many (%d) iterations in XYtoUV...", iter);
- }
-
int thresh = Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin;
- if (thresh){
- //Msg(WARNING, "(U,V) thresholded in XYtoUV (surface mesh may be wrong)");
- if(Unew > umax) *U = umax;
- if(Vnew > vmax) *V = vmax;
- if(Unew < umin) *U = umin;
- if(Vnew < vmin) *V = vmin;
+ if(iter == MaxIter || thresh){
+ if(relax < 1.e-6){
+ Msg(GERROR, "Could not converge: surface mesh will probably be wrong");
+ if(Unew > umax) *U = umax;
+ if(Vnew > vmax) *V = vmax;
+ if(Unew < umin) *U = umin;
+ if(Vnew < vmin) *V = vmin;
+ find_bestuv (s, *X, *Y, U, V, Z, 30);
+ P = InterpolateSurface (s, *U, *V, 0, 0);
+ *X = P.Pos.X;
+ *Y = P.Pos.Y;
+ *Z = P.Pos.Z;
+ }
+ else{
+ Msg(INFO, "Relaxation factor = %g", 0.75*relax);
+ XYtoUV (s, X, Y, U, V, Z, 0.75*relax);
+ }
}
-#if 1
- if (iter == MaxIter || thresh){
- Msg(WARNING, "Entering rescue mode in XYtoUV: surface mesh may be wrong");
- find_bestuv (s, *X, *Y, U, V, Z, 30);
- P = InterpolateSurface (s, *U, *V, 0, 0);
- *X = P.Pos.X;
- *Y = P.Pos.Y;
- *Z = P.Pos.Z;
- }
-#endif
}
+
int Oriente (List_T * cu, double n[3]){
int N, i, a, b, c;
double cosa, sina, sum, v[3], w[3], u[3];
diff --git a/Mesh/Utils.h b/Mesh/Utils.h
index 12c5a589c0b7dd4a2245427f1b6d2c647e0907af..0b71f07145363a2d7ced3d9ca118612695d62719 100644
--- a/Mesh/Utils.h
+++ b/Mesh/Utils.h
@@ -7,8 +7,9 @@ void MeanPlane(List_T *point, Surface *s);
void find_bestuv (Surface * s, double X, double Y,
double *U, double *V, double *Z, int N);
void XYtoUV (Surface * s, double *X, double *Y,
- double *U, double *V, double *Z);
-void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V);
+ double *U, double *V, double *Z, double relax);
+void XYZtoUV (Surface *s, double X, double Y, double Z,
+ double *U, double *V, double relax);
int Oriente (List_T * cu, double n[3]);
double angle_3p (Vertex * V, Vertex * P1, Vertex * P2);
double angle_plan (Vertex * V, Vertex * P1, Vertex * P2, double n[3]);