diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp index 933434a1a1fb271be246503057303e368a87f7f7..7b2f62dd65be7849a8053ee7a85d1462741e1c89 100644 --- a/Geo/Geo.cpp +++ b/Geo/Geo.cpp @@ -1,4 +1,4 @@ -// $Id: Geo.cpp,v 1.116 2008-06-20 12:15:44 remacle Exp $ +// $Id: Geo.cpp,v 1.117 2008-06-20 16:25:38 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -2963,43 +2963,23 @@ bool ProjectPointOnCurve(Curve *c, Vertex *v, Vertex *RES, Vertex *DER) bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2]) { double x[3] = { 0.5, u[0], u[1] }; - double res[3]; - Vertex vv; int check; SURFACE = s; VERTEX = &p; - double UMIN = 0.; - double UMAX = 1.; - double VMIN = 0.; - double VMAX = 1.; newt(x, 2, &check, projectPS); - vv = InterpolateSurface(s, x[1], x[2], 0, 0); - projectPS(2, x, res); - - double resid = sqrt (res[1]*res[1]+res[2]*res[2]); - /* - - double eps = 1.e-6; - int nn = 0; - while(++nn < 20) { - newt(x, 2, &check, projectPS); - vv = InterpolateSurface(s, x[1], x[2], 0, 0); - if(x[1] > UMIN-eps && x[1] < UMAX+eps && x[2] > VMIN-eps && x[2] < VMAX+eps) - break; - x[1] = UMIN + (UMAX - UMIN) * ((rand() % 10000) / 10000.); - x[2] = VMIN + (VMAX - VMIN) * ((rand() % 10000) / 10000.); - } - */ + Vertex vv = InterpolateSurface(s, x[1], x[2], 0, 0); + double res[3]; + projectPS(2, x, res); + double resid = sqrt(res[1] * res[1] + res[2] * res[2]); + p.Pos.X = vv.Pos.X; p.Pos.Y = vv.Pos.Y; p.Pos.Z = vv.Pos.Z; u[0] = x[1]; u[1] = x[2]; - if (resid > 1.e-6){ - // printf("error %12.5E\n",resid); + if(resid > 1.e-6) return false; - } return true; } diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp index 659348380ae69ecac217528b261325a1426834cc..0424544bb59eb1fe17c6f74c45239693786f1529 100644 --- a/Geo/MElement.cpp +++ b/Geo/MElement.cpp @@ -1,4 +1,4 @@ -// $Id: MElement.cpp,v 1.73 2008-06-20 12:15:44 remacle Exp $ +// $Id: MElement.cpp,v 1.74 2008-06-20 16:25:38 geuzaine Exp $ // // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle // @@ -788,44 +788,6 @@ void MTetrahedron::jac(int ord, MVertex *vs[], double uu, double vv, double ww, #endif } - - -int MTriangle6::getNumEdgesRep(){ return 3 * 6; } - -void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) -{ - n[0] = n[1] = getFace(0).normal(); - int N = getNumEdgesRep() / 3; - if (num < N){ - SPoint3 pnt1, pnt2; - pnt((double)num / N, 0., 0,pnt1); - pnt((double)(num + 1) / N, 0., 0,pnt2); - x[0] = pnt1.x(); x[1] = pnt2.x(); - y[0] = pnt1.y(); y[1] = pnt2.y(); - z[0] = pnt1.z(); z[1] = pnt2.z(); - return; - } - if (num < 2 * N){ - SPoint3 pnt1, pnt2; - num -= N; - pnt(1. - (double)num / N, (double)num / N, 0,pnt1); - pnt(1. - (double)(num + 1) / N, (double)(num + 1) / N, 0,pnt2); - x[0] = pnt1.x(); x[1] = pnt2.x(); - y[0] = pnt1.y(); y[1] = pnt2.y(); - z[0] = pnt1.z(); z[1] = pnt2.z(); - return ; - } - { - SPoint3 pnt1, pnt2; - num -= 2 * N; - pnt(0, (double)num / N, 0,pnt1); - pnt(0, (double)(num + 1) / N, 0,pnt2); - x[0] = pnt1.x(); x[1] = pnt2.x(); - y[0] = pnt1.y(); y[1] = pnt2.y(); - z[0] = pnt1.z(); z[1] = pnt2.z(); - } -} - const int numSubEdges = 12; int MTriangleN::getNumFacesRep(){ return numSubEdges * numSubEdges; } diff --git a/Geo/MElement.h b/Geo/MElement.h index 707c9ab7c6eb98de7cca22f31297e7648986039f..449edf5d5d615f3ca85cdce2127c4f7a9a03a459 100644 --- a/Geo/MElement.h +++ b/Geo/MElement.h @@ -513,16 +513,16 @@ class MTriangle6 : public MTriangle { return getVertex(map[num]); } virtual int getNumEdgeVertices(){ return 3; } - virtual int getNumEdgesRep(); - virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n); - //{ - // static const int e[6][2] = { - // {0, 3}, {3, 1}, - // {1, 4}, {4, 2}, - // {2, 5}, {5, 0} - // }; - // _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0); - //} + virtual int getNumEdgesRep(){ return 6; } + virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) + { + static const int e[6][2] = { + {0, 3}, {3, 1}, + {1, 4}, {4, 2}, + {2, 5}, {5, 0} + }; + _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0); + } virtual int getNumFacesRep(){ return 4; } virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n) {