diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 3909261eab51f9cfb6e4350cb07d217f0c487079..8dfc35c287eb24a567dda1ddc2cce13915943f6e 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.30 2006-12-16 14:37:19 geuzaine Exp $
+// $Id: GFace.cpp,v 1.31 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -425,7 +425,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
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};
+ double initu[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 1., 0., 0., 1.};
+ double initv[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0., 1., 0., 1.};
Range<double> ru = parBounds(0);
Range<double> rv = parBounds(1);
@@ -437,13 +438,18 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
for(int i = 0; i < NumInitGuess; i++){
for(int j = 0; j < NumInitGuess; j++){
- U = init[i];
- V = init[j];
+ U = initu[i];
+ V = initv[j];
err = 1.0;
iter = 1;
-
+
+ GPoint P = point(U,V);
+ err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
+ if (err2 < 1.e-8 * CTX.lc)return;
+
+
while(err > Precision && iter < MaxIter) {
- GPoint P = point(U,V);
+ P = point(U,V);
Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V));
mat[0][0] = der.left().x();
mat[0][1] = der.left().y();
@@ -451,6 +457,12 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
mat[1][0] = der.right().x();
mat[1][1] = der.right().y();
mat[1][2] = der.right().z();
+
+// printf("X = %g Y = %g Z = %g U %g V %g deru = %g %g %g derv = %g %g %g\n",P.x(),P.y(),P.z(),U,V
+// ,der.left().x(),der.left().y(),der.left().z()
+// ,der.right().x(),der.right().y(),der.right().z());
+// getchar();
+
mat[2][0] = 0.;
mat[2][1] = 0.;
mat[2][2] = 0.;
@@ -484,7 +496,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
if(relax < 1.e-6)
Msg(GERROR, "Could not converge: surface mesh will be wrong");
else {
- Msg(INFO, "Relaxation factor = %g", 0.75 * relax);
+ Msg(INFO, "point %g %g %g : Relaxation factor = %g",X,Y,Z, 0.75 * relax);
XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
}
}
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index e72f63a5b5c0b4f2f1937751820d3a5cfb9f5c62..da739dc63a791b9842f414de76a5fbf9f88abe87 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.66 2006-12-02 22:05:15 geuzaine Exp $
+// $Id: Geo.cpp,v 1.67 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1562,14 +1562,17 @@ Curve *CreateReversedCurve(Curve * c)
newc->Typ = MSH_SEGM_ELLI_INV;
if(c->Typ == MSH_SEGM_ELLI_INV)
newc->Typ = MSH_SEGM_ELLI;
- newc->Method = c->Method;
- newc->degre = c->degre;
+
+
newc->beg = c->end;
newc->end = c->beg;
+ newc->Method = c->Method;
+ newc->degre = c->degre;
newc->ubeg = 1. - c->uend;
newc->uend = 1. - c->ubeg;
End_Curve(newc);
+
Curve **pc;
if((pc = (Curve **) Tree_PQuery(THEM->Curves, &newc))) {
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index b3cc8aecf32ff72843c60042231dad807d4ffc82..0cc31ca392801b798057390f9cc6596283c8cb7a 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.6 2006-12-02 19:29:36 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.7 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -262,6 +262,10 @@ Vertex InterpolateCurve(Curve * Curve, double u, int derivee)
#define TRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
(1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)
+#define DUTRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
+ (-1.)*c4 +c2 -( (-1.)*(1.-v)*s1+ (1.-v)*s2+ v*s3 -v*s4)
+#define DVTRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
+ (-1.)*c1+ c3-( (-1.)*(1.-u)*s1 -u*s2+ u*s3+ (1.-u)*s4)
Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
Vertex s1, Vertex s2, Vertex s3, Vertex s4,
@@ -280,11 +284,12 @@ Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, s4.Pos.Z, u, v);
return (V);
}
-
/* Interpolation transfinie sur un triangle : TRAN_QUA avec s1=s4=c4
f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3 */
-#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3);
+#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3);
+#define DUTRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) c2+ -( (1.-v)*s2+ v*s3);
+#define DVTRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) (-1.)*c1+ c3-(u*(-1) *s2+ u*s3);
Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
Vertex s1, Vertex s2, Vertex s3, double u, double v)
@@ -301,7 +306,6 @@ Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, u, v);
return (V);
}
-
void TransfiniteSph(Vertex S, Vertex center, Vertex * T)
{
double r, s, dirx, diry, dirz;
@@ -341,63 +345,16 @@ void Calcule_Z_Plan(void *data, Surface *THESURFACE)
v->Pos.Z = 0.0;
}
-Vertex InterpolateSurface(Surface * s, double u, double v,
- int derivee, int u_v)
+Vertex InterpolateRuledSurface(Surface * s, double u, double v,
+ int derivee, int u_v)
{
Vertex *c1, *c2, T, D[4], V[4], *S[4];
Curve *C[4];
int i, issphere;
double eps = 1.e-6;
- if(derivee) {
- if(u_v == 1) {
- if(u - eps < 0.0) {
- D[0] = InterpolateSurface(s, u, v, 0, 0);
- D[1] = InterpolateSurface(s, u + eps, v, 0, 0);
- }
- else {
- D[0] = InterpolateSurface(s, u - eps, v, 0, 0);
- D[1] = InterpolateSurface(s, u, v, 0, 0);
- }
- }
- else if(u_v == 2) {
- if(v - eps < 0.0) {
- D[0] = InterpolateSurface(s, u, v, 0, 0);
- D[1] = InterpolateSurface(s, u, v + eps, 0, 0);
- }
- else {
- D[0] = InterpolateSurface(s, u, v - eps, 0, 0);
- D[1] = InterpolateSurface(s, u, v, 0, 0);
- }
- }
- else {
- Msg(WARNING, "Arbitrary InterpolateSurface for derivative not done");
- }
- T.Pos.X = (D[1].Pos.X - D[0].Pos.X) / eps;
- T.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / eps;
- T.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / eps;
- return T;
- }
-
- Vertex x(u, v, .0);
- Vertex *xx = &x;
-
- if(s->Extrude && s->Extrude->geo.Mode == EXTRUDED_ENTITY &&
- s->Typ != MSH_SURF_PLAN) {
- Curve *c = FindCurve(s->Extrude->geo.Source);
- Vertex v1 = InterpolateCurve(c, u, 0);
- s->Extrude->Extrude(v, v1.Pos.X, v1.Pos.Y, v1.Pos.Z);
- return v1;
- }
-
switch (s->Typ) {
- case MSH_SURF_PLAN:
-
- Calcule_Z_Plan(&xx, s);
- //Projette_Inverse(&xx, &dum);
- return x;
-
case MSH_SURF_REGL:
issphere = 1;
for(i = 0; i < 4; i++) {
@@ -421,28 +378,67 @@ Vertex InterpolateSurface(Surface * s, double u, double v,
S[1] = C[1]->beg;
S[2] = C[2]->beg;
S[3] = C[3]->beg;
- V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
- V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
- V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
- V[3] = InterpolateCurve(C[3], C[3]->ubeg + (C[3]->uend - C[3]->ubeg) * (1. - v), 0);
+ if (C[0]->Num < 0)
+ {
+ Curve *C0 = FindCurve(-C[0]->Num);
+ V[0] = InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1.-u), 0);
+ }
+ else
+ V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
+ if (C[1]->Num < 0)
+ {
+ Curve *C1 = FindCurve(-C[1]->Num);
+ V[1] = InterpolateCurve(C1, C1->ubeg + (C1->uend - C1->ubeg) * (1.-v), 0);
+ }
+ else
+ V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
+ if (C[2]->Num < 0)
+ {
+ Curve *C2 = FindCurve(-C[2]->Num);
+ V[2] = InterpolateCurve(C2, C2->ubeg + (C2->uend - C2->ubeg) * u, 0);
+ }
+ else
+ V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
+ if (C[3]->Num < 0)
+ {
+ Curve *C3 = FindCurve(-C[3]->Num);
+ V[3] = InterpolateCurve(C3, C3->ubeg + (C3->uend - C3->ubeg) * v, 0);
+ }
+ else
+ V[3] = InterpolateCurve(C[3], C[3]->ubeg + (C[3]->uend - C[3]->ubeg) * (1. - v), 0);
+
T = TransfiniteQua(V[0], V[1], V[2], V[3], *S[0], *S[1], *S[2], *S[3], u, v);
- if(issphere)
- TransfiniteSph(*S[0], *c1, &T);
+// if (u == 1 && v == 0)
+// {
+// printf("beg and end of %d : %g %g\n", C[1]->Num, C[1]->ubeg ,C[1]->uend );
+// printf("V = %g %g %g , %g %g %g , %g %g %g , %g %g %g \n",
+// V[0].Pos.X,V[0].Pos.Y,V[0].Pos.Z,
+// V[1].Pos.X,V[1].Pos.Y,V[1].Pos.Z,
+// V[2].Pos.X,V[2].Pos.Y,V[2].Pos.Z,
+// V[3].Pos.X,V[3].Pos.Y,V[3].Pos.Z);
+// printf("S = %g %g %g , %g %g %g , %g %g %g , %g %g %g \n",
+// S[0]->Pos.X,S[0]->Pos.Y,S[0]->Pos.Z,
+// S[1]->Pos.X,S[1]->Pos.Y,S[1]->Pos.Z,
+// S[2]->Pos.X,S[2]->Pos.Y,S[2]->Pos.Z,
+// S[3]->Pos.X,S[3]->Pos.Y,S[3]->Pos.Z);
+// printf("%g %g %g for %g %g\n",T.Pos.X,T.Pos.Y,T.Pos.Z,u,v);
+// }
+
+ // if(issphere)
+ // TransfiniteSph(*S[0], *c1, &T);
+
return (T);
-
- case MSH_SURF_NURBS:
- return InterpolateNurbsSurface(s, u, v);
-
+
case MSH_SURF_TRIC:
issphere = 1;
for(i = 0; i < 3; i++) {
List_Read(s->Generatrices, i, &C[i]);
if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV) {
- issphere = 0;
+ issphere = 0;
}
else if(issphere) {
- if(!i) {
+ if(!i) {
List_Read(C[i]->Control_Points, 1, &c1);
}
else {
@@ -452,20 +448,83 @@ Vertex InterpolateSurface(Surface * s, double u, double v,
}
}
}
-
+
S[0] = C[0]->beg;
S[1] = C[1]->beg;
S[2] = C[2]->beg;
V[0] = InterpolateCurve(C[0], u, 0);
V[1] = InterpolateCurve(C[1], v, 0);
V[2] = InterpolateCurve(C[2], 1. - u, 0);
-
+
T = TransfiniteTri(V[0], V[1], V[2], *S[0], *S[1], *S[2], u, v);
if(issphere)
- TransfiniteSph(*S[0], *c1, &T);
-
+ TransfiniteSph(*S[0], *c1, &T);
return (T);
+ }
+}
+
+Vertex InterpolateSurface(Surface * s, double u, double v,
+ int derivee, int u_v)
+{
+ Vertex *c1, *c2, T, D[4], V[4], *S[4];
+ Curve *C[4];
+ int i, issphere;
+ double eps = 1.e-6;
+
+ if(s->Extrude && s->Extrude->geo.Mode == EXTRUDED_ENTITY &&
+ s->Typ != MSH_SURF_PLAN) {
+ Curve *c = FindCurve(s->Extrude->geo.Source);
+ Vertex v1 = InterpolateCurve(c, u, 0);
+ s->Extrude->Extrude(v, v1.Pos.X, v1.Pos.Y, v1.Pos.Z);
+ return v1;
+ }
+
+
+ if(derivee) {
+ if(u_v == 1) {
+ if(u - eps < 0.0) {
+ D[0] = InterpolateSurface(s, u, v, 0, 0);
+ D[1] = InterpolateSurface(s, u + eps, v, 0, 0);
+ }
+ else {
+ D[0] = InterpolateSurface(s, u - eps, v, 0, 0);
+ D[1] = InterpolateSurface(s, u, v, 0, 0);
+ }
+ }
+ else if(u_v == 2) {
+ if(v - eps < 0.0) {
+ D[0] = InterpolateSurface(s, u, v, 0, 0);
+ D[1] = InterpolateSurface(s, u, v + eps, 0, 0);
+ }
+ else {
+ D[0] = InterpolateSurface(s, u, v - eps, 0, 0);
+ D[1] = InterpolateSurface(s, u, v, 0, 0);
+ }
+ }
+ else {
+ Msg(WARNING, "Arbitrary InterpolateSurface for derivative not done");
+ }
+ T.Pos.X = (D[1].Pos.X - D[0].Pos.X) / eps;
+ T.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / eps;
+ T.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / eps;
+ return T;
+ }
+
+ if (s->Typ == MSH_SURF_REGL || s->Typ == MSH_SURF_TRIC) return InterpolateRuledSurface(s,u,v,derivee,u_v);
+
+ Vertex x(u, v, .0);
+ Vertex *xx = &x;
+
+ switch (s->Typ) {
+
+ case MSH_SURF_PLAN:
+
+ Calcule_Z_Plan(&xx, s);
+ //Projette_Inverse(&xx, &dum);
+ return x;
+ case MSH_SURF_NURBS:
+ return InterpolateNurbsSurface(s, u, v);
default:
Msg(GERROR, "Unknown surface type in interpolation");
T.Pos.X = T.Pos.Y = T.Pos.Z = 0.0;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index cb0c703c4956da41a28bcf390dfe935759f2f588..bbc667e492a9f193b5d8fa3d464c261c49b76a78 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.24 2006-11-27 22:22:13 geuzaine Exp $
+// $Id: MElement.cpp,v 1.25 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -23,6 +23,7 @@
#include "MElement.h"
#include "GEntity.h"
#include "Numeric.h"
+#include "Message.h"
int edges_tetra[6][2] = {
{0, 1},
@@ -402,3 +403,67 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
fprintf(fp, "\n");
}
}
+
+bool MTriangle::invertmappingXY(double *p, double *uv, double tol)
+{
+ double mat[2][2];
+ double b[2], dum;
+ getMat(mat);
+ b[0] = p[0] - getVertex(0)->x();
+ b[1] = p[1] - getVertex(0)->y();
+ sys2x2(mat, b, uv);
+ if(uv[0] >= -tol &&
+ uv[1] >= -tol &&
+ uv[0] <= 1. + tol &&
+ uv[1] <= 1. + tol &&
+ 1. - uv[0] - uv[1] > -tol) {
+ return true;
+ }
+ return false;
+}
+
+
+double MTriangle::getSurfaceXY() const
+{
+ const double x1 =_v[0]->x();
+ const double x2 =_v[1]->x();
+ const double x3 =_v[2]->x();
+ const double y1 =_v[0]->y();
+ const double y2 =_v[1]->y();
+ const double y3 =_v[2]->y();
+
+ const double v1 [2] = {x2-x1,y2-y1};
+ const double v2 [2] = {x3-x1,y3-y1};
+
+ double s = v1[0]*v2[1] - v1[1]*v2[0];
+ return s*0.5;
+
+}
+
+void MTriangle::circumcenterXY(double *res) const
+{
+ double d, a1, a2, a3;
+
+ const double x1 =_v[0]->x();
+ const double x2 =_v[1]->x();
+ const double x3 =_v[2]->x();
+ const double y1 =_v[0]->y();
+ const double y2 =_v[1]->y();
+ const double y3 =_v[2]->y();
+
+ d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
+ if(d == 0.0) {
+ Msg(WARNING, "Colinear points in circum circle computation");
+ res[0] = res[1] = -99999.;
+ return ;
+ }
+
+ a1 = x1 * x1 + y1 * y1;
+ a2 = x2 * x2 + y2 * y2;
+ a3 = x3 * x3 + y3 * y3;
+ res[0] = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
+ res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
+
+ return ;
+}
+
diff --git a/Geo/MElement.h b/Geo/MElement.h
index b14436311948d13e4deb1d9d5dfde5cbf9babf76..0d8bdc168eb54fa49eb29d6fa72e21fad80f1756 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -249,6 +249,16 @@ class MTriangle : public MElement {
}
~MTriangle(){}
virtual int getDim(){ return 2; }
+ virtual void getMat(double mat[2][2])
+ {
+ mat[0][0] = _v[1]->x() - _v[0]->x();
+ mat[0][1] = _v[2]->x() - _v[0]->x();
+ mat[1][0] = _v[1]->y() - _v[0]->y();
+ mat[1][1] = _v[2]->y() - _v[0]->y();
+ }
+ void circumcenterXY(double *res) const;
+ double getSurfaceXY() const;
+ bool invertmappingXY(double *p, double *uv, double tol = 1.e-8);
virtual int getNumVertices(){ return 3; }
virtual int getNumPrimaryVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
diff --git a/Geo/Makefile b/Geo/Makefile
index 5b49bdbd0c92ce9c04660f0acd7c6e4fdbb926fc..e304a138a0c72aee1a81e7b6bd96afeb924828c3 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.121 2006-12-03 03:19:55 geuzaine Exp $
+# $Id: Makefile,v 1.122 2006-12-20 15:50:57 remacle Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -27,9 +27,9 @@ INCLUDE = -I../Common -I../DataStr -I../Geo -I../Mesh -I../Numeric\
CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = GEntity.cpp\
- GVertex.cpp GEdge.cpp GEdgeLoop.cpp GFace.cpp GRegion.cpp\
- gmshEdge.cpp gmshFace.cpp gmshRegion.cpp\
- OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp\
+ GVertex.cpp GEdge.cpp GEdgeLoop.cpp GFace.cpp GRegion.cpp\
+ gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp\
+ OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp\
projectionFace.cpp\
GModel.cpp\
GModelIO_Geo.cpp\
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index bf61d26418c83c07eba136d09f1616bc981cc3de..cc7257fb15baddf6d45a61bdc537af30a9516ab4 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.16 2006-11-29 16:57:01 remacle Exp $
+// $Id: OCCFace.cpp,v 1.17 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -137,8 +137,10 @@ GPoint OCCFace::closestPoint(const SPoint3 & qp) const
gp_Pnt pnt(qp.x(),qp.y(),qp.z());
GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
if (!proj.NbPoints())
- Msg(GERROR,"OCC Project Point on Surface FAIL");
-
+ {
+ Msg(GERROR,"OCC Project Point on Surface FAIL");
+ return GPoint(0,0);
+ }
pnt = proj.NearestPoint();
double pp[2];
proj.LowerDistanceParameters (pp[0], pp[1]);
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 65e849a3ae55e4bc0c43e441108601f608ea53e1..4a05989c1a2248077eb4407c4d6b8106fc9fad6f 100644
--- a/Geo/gmshEdge.cpp
+++ b/Geo/gmshEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshEdge.cpp,v 1.22 2006-12-16 01:25:58 geuzaine Exp $
+// $Id: gmshEdge.cpp,v 1.23 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -19,6 +19,7 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
+#include "GFace.h"
#include "gmshEdge.h"
#include "Geo.h"
#include "GeoInterpolation.h"
@@ -136,3 +137,73 @@ int gmshEdge::minimumDrawSegments () const
else
return 10 * n;
}
+SPoint2 gmshEdge::reparamOnFace(GFace *face, double epar,int dir) const
+{
+ Surface *s = (Surface*) face->getNativePtr();
+ if (s->Typ == MSH_SURF_REGL)
+ {
+ double U,V;
+ Curve *C[4];
+ int iPos;
+ for(int i = 0; i < 4; i++) {
+ List_Read(s->Generatrices, i, &C[i]);
+ }
+
+ if (C[0]->Num== c->Num) {
+ U = (epar - C[0]->ubeg) / (C[0]->uend - C[0]->ubeg) ;
+ V = 0;
+ iPos = 0;
+ }
+ else if (C[0]->Num== -c->Num) {
+ U = (C[0]->uend - epar - C[0]->ubeg) / (C[0]->uend - C[0]->ubeg) ;
+ V = 0;
+ iPos = 0;
+ }
+ else if (C[1]->Num== c->Num) {
+ V = (epar - C[1]->ubeg) / (C[1]->uend - C[1]->ubeg) ;
+ U = 1;
+ iPos = 1;
+ }
+ else if (C[1]->Num== - c->Num ) {
+ V = (C[1]->uend - epar - C[1]->ubeg) / (C[1]->uend - C[1]->ubeg) ;
+ U = 1;
+ iPos = 1;
+ }
+ else if (C[2]->Num== c->Num) {
+ U = 1 - (epar - C[2]->ubeg) / (C[2]->uend - C[2]->ubeg) ;
+ V = 1;
+ iPos = 2;
+ }
+ else if (C[2]->Num== -c->Num) {
+ U = 1 - ( C[2]->uend -epar - C[2]->ubeg) / (C[2]->uend - C[2]->ubeg) ;
+ V = 1;
+ iPos = 2;
+ }
+ else if (C[3]->Num== c->Num) {
+ V = 1-(epar - C[3]->ubeg) / (C[3]->uend - C[3]->ubeg) ;
+ U = 0;
+ iPos = 3;
+ }
+ else if (C[3]->Num== -c->Num) {
+ V = 1-(C[3]->uend - epar - C[3]->ubeg) / (C[3]->uend - C[3]->ubeg) ;
+ U = 0;
+ iPos = 3;
+ }
+ else throw;
+
+// GPoint p1 = point(epar);
+// GPoint p2 = face->point(U,V);
+
+// printf("iPos %1d face %2d (%2d %2d %2d %2d) (%8.3f %8.3f) curve %2d (%8.3f) %8.3f %8.3f %8.3f vs %8.3f %8.3f %8.3f D = %8.3f\n",
+// iPos,face->tag(),C[0]->Num,C[1]->Num,C[2]->Num,C[3]->Num,U,V,
+// tag(),epar,
+// p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),
+// sqrt((p1.x()-p2.x())*(p1.x()-p2.x())+(p1.y()-p2.y())*(p1.y()-p2.y())+(p1.z()-p2.z())*(p1.z()-p2.z())));
+
+ return SPoint2(U,V);
+ }
+ else
+ return GEdge::reparamOnFace(face,epar,dir);
+
+
+}
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index fe214653bf792c2e4269a2108a13924a46ff8345..3151c4c7b0cf7a31f31ac9de4906c6ffd25e035b 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -50,6 +50,7 @@ class gmshEdge : public GEdge {
virtual int minimumMeshSegments () const;
virtual int minimumDrawSegments () const;
virtual void resetMeshAttributes ();
+ virtual SPoint2 reparamOnFace(GFace *face, double epar,int dir) const ;
};
#endif
diff --git a/Geo/gmshVertex.cpp b/Geo/gmshVertex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9fd1354286d9722f14e81b2f028a566afdc604c5
--- /dev/null
+++ b/Geo/gmshVertex.cpp
@@ -0,0 +1,74 @@
+#include "GFace.h"
+#include "gmshVertex.h"
+#include "Geo.h"
+#include "GeoInterpolation.h"
+#include "Context.h"
+
+SPoint2 gmshVertex::reparamOnFace ( GFace *face , int dir) const
+{
+ Surface *s = (Surface*) face->getNativePtr();
+
+
+ if (s->Typ == MSH_SURF_REGL)
+ {
+ const double corners[4][2] = {{0,0},{0,1},{1,1},{1,0}};
+
+ double U,V;
+ Curve *C[4];
+ int iPos;
+ for(int i = 0; i < 4; i++) {
+ List_Read(s->Generatrices, i, &C[i]);
+ }
+
+
+ if ((C[0]->beg == v && C[3]->beg == v) ||
+ (C[0]->end == v && C[3]->beg == v) ||
+ (C[0]->beg == v && C[3]->end == v) ||
+ (C[0]->end == v && C[3]->end == v))
+ {
+ U = V = 0;
+ }
+ else if ((C[0]->beg == v && C[1]->beg == v) ||
+ (C[0]->end == v && C[1]->beg == v) ||
+ (C[0]->beg == v && C[1]->end == v) ||
+ (C[0]->end == v && C[1]->end == v))
+ {
+ U = 1;
+ V = 0;
+ }
+ else if ((C[2]->beg == v && C[1]->beg == v) ||
+ (C[2]->end == v && C[1]->beg == v) ||
+ (C[2]->beg == v && C[1]->end == v) ||
+ (C[2]->end == v && C[1]->end == v))
+ {
+ U = 1;
+ V = 1;
+ }
+ else if ((C[2]->beg == v && C[3]->beg == v) ||
+ (C[2]->end == v && C[3]->beg == v) ||
+ (C[2]->beg == v && C[3]->end == v) ||
+ (C[2]->end == v && C[3]->end == v))
+ {
+ U = 0;
+ V = 1;
+ }
+ else
+ throw;
+// GPoint p1 = point();
+// GPoint p2 = face->point(U,V);
+
+// printf("face %2d (%8.3f %8.3f) point %2d %8.3f %8.3f %8.3f vs %8.3f %8.3f %8.3f D = %8.3f\n",
+// face->tag(),U,V,
+// tag(),
+// p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),
+// sqrt((p1.x()-p2.x())*(p1.x()-p2.x())+(p1.y()-p2.y())*(p1.y()-p2.y())+(p1.z()-p2.z())*(p1.z()-p2.z())));
+
+
+
+ return SPoint2(U,V);
+ }
+ else
+ {
+ return GVertex::reparamOnFace(face,dir);
+ }
+}
diff --git a/Geo/gmshVertex.h b/Geo/gmshVertex.h
index 256f438817415cf517c5a18c5dd0af0841f3b8aa..db24110a8ed5ed5a702236f9593105ce0eaaef2f 100644
--- a/Geo/gmshVertex.h
+++ b/Geo/gmshVertex.h
@@ -62,6 +62,7 @@ class gmshVertex : public GVertex {
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return v; }
virtual void setPrescribedMeshSizeAtVertex(double l) {meshSize = l;v->lc = meshSize;}
+ virtual SPoint2 reparamOnFace ( GFace *gf , int) const;
};
#endif
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 39ee4e2089e7cabae289668e864edf415ec5d233..521ae19c41683747028768429b3d568b18e33e18 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.156 2006-12-16 15:48:51 geuzaine Exp $
+# $Id: Makefile,v 1.157 2006-12-20 15:50:57 remacle Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -32,18 +32,19 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Generator.cpp \
meshGEdge.cpp \
- meshGEdgeExtruded.cpp \
+ meshGEdgeExtruded.cpp \
meshGFace.cpp \
- meshGFaceTransfinite.cpp \
- meshGFaceExtruded.cpp \
+ meshGFaceTransfinite.cpp \
+ meshGFaceExtruded.cpp \
meshGRegion.cpp \
- meshGRegionDelaunayInsertion.cpp \
- meshGRegionTransfinite.cpp \
- meshGRegionExtruded.cpp \
- DivideAndConquer.cpp \
- BackgroundMesh.cpp \
- BDS.cpp \
- SecondOrder.cpp
+ meshGRegionDelaunayInsertion.cpp \
+ meshGFaceDelaunayInsertion.cpp \
+ meshGRegionTransfinite.cpp \
+ meshGRegionExtruded.cpp \
+ DivideAndConquer.cpp \
+ BackgroundMesh.cpp \
+ BDS.cpp \
+ SecondOrder.cpp
OBJ = ${SRC:.cpp=.o}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 19fc7e6da7120bd085c618597b1f1f121fd89ebd..7864267e897dfc7e44269bc2ba9b0e5cdfbd1e2c 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.42 2006-12-15 03:15:32 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.43 2006-12-20 15:50:57 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -36,24 +36,6 @@
extern Context_T CTX;
-class fromCartesianToParametric
-{
- GFace *gf;
-public :
- fromCartesianToParametric ( GFace *_gf )
- : gf(_gf){}
- void operator () (MVertex * v)
- {
-
- GEntity *ge = v->onWhat();
-
- SPoint2 param = gf->parFromPoint (SPoint3(v->x(),v->y(),v->z()));
- v->x() = param.x();
- v->y() = param.y();
- v->z() = 0.0;
- }
-};
-
void computeEdgeLoops (const GFace *gf,
std::vector<MVertex*> & all_mvertices,
std::vector<int> &indices)
@@ -559,6 +541,8 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
bool gmsh2DMeshGenerator ( GFace *gf )
{
+ // if (gf->tag() != 575)return true;
+
typedef std::set<MVertex*> v_container ;
v_container all_vertices;
std::map<int, MVertex*>numbered_vertices;
@@ -599,7 +583,21 @@ bool gmsh2DMeshGenerator ( GFace *gf )
while(itv != all_vertices.end())
{
MVertex *here = *itv;
- SPoint2 param = gf->parFromPoint (SPoint3(here->x(),here->y(),here->z()));
+ SPoint2 param;
+ if (here->onWhat()->dim() == 0)
+ {
+ GVertex *gv = (GVertex*)here->onWhat();
+ param=gv->reparamOnFace (gf,1);
+ }
+ else if (here->onWhat()->dim() == 1)
+ {
+ GEdge *ge = (GEdge*)here->onWhat();
+ double UU;
+ here->getParameter(0,UU);
+ param=ge->reparamOnFace (gf,UU,1);
+ }
+ else
+ param = gf->parFromPoint (SPoint3(here->x(),here->y(),here->z()));
// fprintf(fdeb,"%d %g %g %g\n" ,here->getNum(),here->x(),here->y(),here->z());
U_[count] = param.x();
V_[count] = param.y();