diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp index 92165dc0c15f95c1e39b29132297c6dd4a098e3d..389389275c60869d4aa5f03cde87a4d447f76ce3 100644 --- a/Geo/GFace.cpp +++ b/Geo/GFace.cpp @@ -1,6 +1,16 @@ #include "GModel.h" #include "GFace.h" #include "GEdge.h" +#include "Message.h" + +#if defined(HAVE_GSL) +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> +#else +#define NRANSI +#include "nrutil.h" +void dsvdcmp(double **a, int m, int n, double w[], double **v); +#endif GFace::~GFace () { @@ -63,3 +73,225 @@ void GFace::recomputeMeshPartitions() if(part) model()->getMeshPartitions().insert(part); } } + +void GFace::computeMeanPlane() +{ + std::vector<SPoint3> pts; + std::list<GVertex*> verts = vertices(); + std::list<GVertex*>::const_iterator itv = verts.begin(); + for(; itv != verts.end(); itv++){ + const GVertex *v = *itv; + pts.push_back(SPoint3(v->x(), v->y(), v->z())); + } + computeMeanPlane(pts); +} + +void GFace::computeMeanPlane(const std::vector<MVertex*> &points) +{ + std::vector<SPoint3> pts; + for(unsigned int i = 0; i < points.size(); i++) + pts.push_back(SPoint3(points[i]->x(), points[i]->y(), points[i]->z())); + computeMeanPlane(pts); +} + +void GFace::computeMeanPlane(const std::vector<SPoint3> &points) +{ + // The concept of a mean plane computed in the sense of least + // squares is fine for plane surfaces(!), but not really the best + // one for non-plane surfaces. Indeed, imagine a quarter of a circle + // extruded along a very small height: the mean plane will be in the + // plane of the circle... Until we implement something better, there + // is a test in this routine that checks the coherence of the + // computation for non-plane surfaces and reverts to a basic mean + // plane approximatation if the check fails. + + double xm = 0., ym = 0., zm = 0.; + int ndata = points.size(); + int na = 3; + for(int i = 0; i < ndata; i++) { + xm += points[i].x(); + ym += points[i].y(); + zm += points[i].z(); + } + xm /= (double)ndata; + ym /= (double)ndata; + zm /= (double)ndata; + + int min; + double res[4], svd[3]; +#if defined(HAVE_GSL) + gsl_matrix *U = gsl_matrix_alloc(ndata, na); + gsl_matrix *V = gsl_matrix_alloc(na, na); + gsl_vector *W = gsl_vector_alloc(na); + gsl_vector *TMPVEC = gsl_vector_alloc(na); + for(int i = 0; i < ndata; i++) { + gsl_matrix_set(U, i, 0, points[i].x() - xm); + gsl_matrix_set(U, i, 1, points[i].y() - ym); + gsl_matrix_set(U, i, 2, points[i].z() - zm); + } + gsl_linalg_SV_decomp(U, V, W, TMPVEC); + svd[0] = gsl_vector_get(W, 0); + svd[1] = gsl_vector_get(W, 1); + svd[2] = gsl_vector_get(W, 2); + if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2])) + min = 0; + else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2])) + min = 1; + else + min = 2; + res[0] = gsl_matrix_get(V, 0, min); + res[1] = gsl_matrix_get(V, 1, min); + res[2] = gsl_matrix_get(V, 2, min); + norme(res); + gsl_matrix_free(U); + gsl_matrix_free(V); + gsl_vector_free(W); + gsl_vector_free(TMPVEC); +#else + double **U = dmatrix(1, ndata, 1, na); + double **V = dmatrix(1, na, 1, na); + double *W = dvector(1, na); + for(int i = 0; i < ndata; i++) { + U[i + 1][1] = points[i].x() - xm; + U[i + 1][2] = points[i].y() - ym; + U[i + 1][3] = points[i].z() - zm; + } + dsvdcmp(U, ndata, na, W, V); + if(fabs(W[1]) < fabs(W[2]) && fabs(W[1]) < fabs(W[3])) + min = 1; + else if(fabs(W[2]) < fabs(W[1]) && fabs(W[2]) < fabs(W[3])) + min = 2; + else + min = 3; + svd[0] = W[1]; + svd[1] = W[2]; + svd[2] = W[3]; + res[0] = V[1][min]; + res[1] = V[2][min]; + res[2] = V[3][min]; + norme(res); + free_dmatrix(U, 1, ndata, 1, na); + free_dmatrix(V, 1, na, 1, na); + free_dvector(W, 1, na); +#endif + + double ex[3], t1[3], t2[3]; + + // check coherence of results for non-plane surfaces + if(geomType() != GEntity::Plane) { + double res2[3], c[3], cosc, sinc, angplan; + double eps = 1.e-3; + + GPoint v1 = point(0.5, 0.5); + GPoint v2 = point(0.5 + eps, 0.5); + GPoint v3 = point(0.5, 0.5 + eps); + t1[0] = v2.x() - v1.x(); + t1[1] = v2.y() - v1.y(); + t1[2] = v2.z() - v1.z(); + t2[0] = v3.x() - v1.x(); + t2[1] = v3.y() - v1.y(); + t2[2] = v3.z() - v1.z(); + norme(t1); + norme(t2); + // prodve(t1, t2, res2); + // Warning: the rest of the code assumes res = t2 x t1, not t1 x t2 (WTF?) + prodve(t2, t1, res2); + norme(res2); + prodve(res, res2, c); + prosca(res, res2, &cosc); + sinc = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]); + angplan = myatan2(sinc, cosc); + angplan = angle_02pi(angplan) * 180. / Pi; + if((angplan > 70 && angplan < 110) || (angplan > 260 && angplan < 280)) { + Msg(INFO, "SVD failed (angle=%g): using rough algo...", angplan); + res[0] = res2[0]; + res[1] = res2[1]; + res[2] = res2[2]; + goto end; + } + } + + ex[0] = ex[1] = ex[2] = 0.0; + if(res[0] == 0.) + ex[0] = 1.0; + else if(res[1] == 0.) + ex[1] = 1.0; + else + ex[2] = 1.0; + + prodve(res, ex, t1); + norme(t1); + prodve(t1, res, t2); + norme(t2); + +end: + res[3] = (xm * res[0] + ym * res[1] + zm * res[2]); + + for(int i = 0; i < 3; i++) + meanPlane.plan[0][i] = t1[i]; + for(int i = 0; i < 3; i++) + meanPlane.plan[1][i] = t2[i]; + for(int i = 0; i < 3; i++) + meanPlane.plan[2][i] = res[i]; + + meanPlane.a = res[0]; + meanPlane.b = res[1]; + meanPlane.c = res[2]; + meanPlane.d = res[3]; + + meanPlane.x = meanPlane.y = meanPlane.z = 0.; + if(fabs(meanPlane.a) >= fabs(meanPlane.b) && + fabs(meanPlane.a) >= fabs(meanPlane.c) ){ + meanPlane.x = meanPlane.d / meanPlane.a; + } + else if(fabs(meanPlane.b) >= fabs(meanPlane.a) && + fabs(meanPlane.b) >= fabs(meanPlane.c)){ + meanPlane.y = meanPlane.d / meanPlane.b; + } + else{ + meanPlane.z = meanPlane.d / meanPlane.c; + } + + Msg(DEBUG1, "Surface: %d", tag()); + Msg(DEBUG2, "SVD : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min); + Msg(DEBUG2, "Plane : (%g x + %g y + %g z = %g)", + meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d); + Msg(DEBUG2, "Normal : (%g , %g , %g )", + meanPlane.a, meanPlane.b, meanPlane.c); + Msg(DEBUG3, "t1 : (%g , %g , %g )", t1[0], t1[1], t1[2]); + Msg(DEBUG3, "t2 : (%g , %g , %g )", t2[0], t2[1], t2[2]); + Msg(DEBUG3, "pt : (%g , %g , %g )", + meanPlane.x, meanPlane.y, meanPlane.z); + + //check coherence for plane surfaces + if(geomType() == GEntity::Plane) { + std::list<GVertex*> verts = vertices(); + std::list<GVertex*>::const_iterator itv = verts.begin(); + for(; itv != verts.end(); itv++){ + const GVertex *v = *itv; + double d = meanPlane.a * v->x() + meanPlane.b * v->y() + + meanPlane.c * v->z() - meanPlane.d; + if(fabs(d) > 1.e-3) { + Msg(GERROR1, "Plane surface %d (%gx+%gy+%gz+%g=0) is not plane!", + v->tag(), meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d); + Msg(GERROR3, "Control point %d = (%g,%g,%g), val=%g", + v->tag(), v->x(), v->y(), v->z(), d); + return; + } + } + } +} + +void GFace::getMeanPlaneData(double VX[3], double VY[3], + double &x, double &y, double &z) const +{ + VX[0] = meanPlane.plan[0][0]; + VX[1] = meanPlane.plan[0][1]; + VX[2] = meanPlane.plan[0][2]; + VY[0] = meanPlane.plan[1][0]; + VY[1] = meanPlane.plan[1][1]; + VY[2] = meanPlane.plan[1][2]; + x = meanPlane.x; + y = meanPlane.y; + z = meanPlane.z; +} diff --git a/Geo/GFace.h b/Geo/GFace.h index f93c1808a14f07cd92d4082c46f4d90ed0aa957e..9708d2c263f9bc02825e4629a18f76acba47653f 100644 --- a/Geo/GFace.h +++ b/Geo/GFace.h @@ -11,8 +11,8 @@ struct mean_plane { double plan[3][3]; - double a,b,c,d; - double x,y,z; + double a, b, c, d; + double x, y, z; }; class GRegion; @@ -25,6 +25,7 @@ class GFace : public GEntity std::list<GEdge *> l_edges; std::list<int> l_dirs; GRegion *r1, *r2; + mean_plane meanPlane; public: GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0) {} @@ -77,9 +78,23 @@ class GFace : public GEntity // recompute the mesh partitions defined on this face. void recomputeMeshPartitions(); + // recompute the mean plane of the surface from a list of points + void computeMeanPlane(const std::vector<MVertex*> &points); + void computeMeanPlane(const std::vector<SPoint3> &points); + + // recompute the mean plane of the surface from its bounding vertices + void computeMeanPlane(); + + // get the mean plane info + void getMeanPlaneData(double VX[3], double VY[3], + double &x, double &y, double &z) const; + + // a crude graphical representation using a "cross" defined by pairs + // of start/end points + std::vector<SPoint3> cross; + std::vector<MTriangle*> triangles; std::vector<MQuadrangle*> quadrangles; - mean_plane mp; }; #endif diff --git a/Geo/Makefile b/Geo/Makefile index b2256c5f110a41f97a976051c4684ffd2028b9c1..a941aba4abefc93a66698c40fff953fbdf23bd3f 100644 --- a/Geo/Makefile +++ b/Geo/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.88 2006-08-12 19:34:15 geuzaine Exp $ +# $Id: Makefile,v 1.89 2006-08-13 02:46:53 geuzaine Exp $ # # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle # @@ -136,7 +136,7 @@ GEdge.o: GEdge.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \ GFace.o: GFace.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \ SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h GPoint.h GEdge.h \ SVector3.h SPoint2.h MElement.h ../Numeric/Numeric.h GFace.h Pair.h \ - GRegion.h + GRegion.h ../Common/Message.h # 1 "/Users/geuzaine/.gmsh/Geo//" GRegion.o: GRegion.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \ SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h GPoint.h GEdge.h \ @@ -204,7 +204,7 @@ gmshFace.o: gmshFace.cpp gmshModel.h GModel.h GVertex.h GEntity.h Range.h \ ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h gmshFace.h \ ../Mesh/Interpolation.h ../Mesh/Vertex.h ../Mesh/Mesh.h CAD.h \ ExtrudeParams.h Geo.h ../Mesh/Create.h ../Mesh/Vertex.h ../Mesh/Mesh.h \ - ../Mesh/Utils.h ../Mesh/Vertex.h ../Mesh/Mesh.h + ../Mesh/Utils.h ../Mesh/Vertex.h ../Mesh/Mesh.h ../Common/Message.h # 1 "/Users/geuzaine/.gmsh/Geo//" gmshRegion.o: gmshRegion.cpp gmshModel.h GModel.h GVertex.h GEntity.h \ Range.h SPoint3.h SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h \ diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp index ef7f9b911d1603183cc081af4a35c2c66d2003f1..bb87d53ebd0a7051f5a6020a2f5e092de292aa9e 100644 --- a/Geo/gmshFace.cpp +++ b/Geo/gmshFace.cpp @@ -7,13 +7,14 @@ #include "Mesh.h" #include "Create.h" #include "Utils.h" +#include "Message.h" extern Mesh *THEM; gmshFace::gmshFace(GModel *m, Surface *face) : GFace(m, face->Num), s(face) { - for (int i=0 ; i < List_Nbr(s->Generatrices) ; i++){ + for(int i = 0 ; i < List_Nbr(s->Generatrices); i++){ Curve *c; List_Read(s->Generatrices, i, &c); GEdge *e = m->edgeByTag(abs(c->Num)); @@ -23,6 +24,9 @@ gmshFace::gmshFace(GModel *m, Surface *face) if(c->Num > 0) l_dirs.push_back(1); else l_dirs.push_back(-1); } + // always compute and store the mean plane for plane surfaces + // (simply using the bounding vertices) + if(s->Typ == MSH_SURF_PLAN) computeMeanPlane(); } gmshFace::gmshFace(GModel *m, int num) @@ -47,21 +51,59 @@ SBoundingBox3d gmshFace::bounds() const std::list<GEdge*>::const_iterator it = l_edges.begin(); SBoundingBox3d res = (*it)->bounds(); ++it; - while (it != l_edges.end()) - { - res += (*it)->bounds(); - ++it; - } + while(it != l_edges.end()){ + res += (*it)->bounds(); + ++it; + } return res; } SVector3 gmshFace::normal(const SPoint2 ¶m) const { - Vertex vu = InterpolateSurface(s, param[0], param[1], 1, 1); - Vertex vv = InterpolateSurface(s, param[0], param[1], 1, 2); - Vertex n = vu % vv; - n.norme(); - return SVector3(n.Pos.X, n.Pos.Y, n.Pos.Z); + if(geomType() != Plane){ + Vertex vu = InterpolateSurface(s, param[0], param[1], 1, 1); + Vertex vv = InterpolateSurface(s, param[0], param[1], 1, 2); + Vertex n = vu % vv; + n.norme(); + return SVector3(n.Pos.X, n.Pos.Y, n.Pos.Z); + } + else{ + // We cannot use InterpolateSurface() for plane surfaces since + // InterpolateSurface() relies on the mean plane, which does *not* + // respect the orientation + GPoint p = point(param); + double t1[3], t2[3], n[3] = {0., 0., 0.}; + Curve *c; + if(List_Nbr(s->Generatrices) > 1){ + List_Read(s->Generatrices, 0, &c); + if(c->beg){ + t1[0] = p.x() - c->beg->Pos.X; + t1[1] = p.y() - c->beg->Pos.Y; + t1[2] = p.z() - c->beg->Pos.Z; + // 1) try to get a point close to 'beg' on the same curve + // 2) if we are unlucky and these two points are aligned with + // (x,y,z), which we know is inside or on the boundary of + // the surface, then get a point from the next edge + // 3) repeat + for(int i = 0; i < List_Nbr(s->Generatrices); i++){ + List_Read(s->Generatrices, i, &c); + Vertex v = InterpolateCurve(c, 0.1, 0); + t2[0] = p.x() - v.Pos.X; + t2[1] = p.y() - v.Pos.Y; + t2[2] = p.z() - v.Pos.Z; + prodve(t1, t2, n); + if(norme(n)) + break; + } + } + } + if(norme(n)) + return SVector3(n[0], n[1], n[2]); + else{ + Msg(WARNING, "Using mean plane normal for surface %d", s->Num); + return SVector3(meanPlane.a, meanPlane.b, meanPlane.c); + } + } } Pair<SVector3,SVector3> gmshFace::firstDer(const SPoint2 ¶m) const @@ -72,40 +114,24 @@ Pair<SVector3,SVector3> gmshFace::firstDer(const SPoint2 ¶m) const SVector3(vv.Pos.X, vv.Pos.Y, vv.Pos.Z)); } - GPoint gmshFace::point(const SPoint2 &pt) const { return point(pt.x(),pt.y()); } - -void computePlaneDatas (const GFace *gf, double VX[3],double VY[3],double &x, double &y, double &z) -{ - VX[0] = gf->mp.plan[0][0]; - VX[1] = gf->mp.plan[0][1]; - VX[2] = gf->mp.plan[0][2]; - VY[0] = gf->mp.plan[1][0]; - VY[1] = gf->mp.plan[1][1]; - VY[2] = gf->mp.plan[1][2]; - x=gf->mp.x; - y=gf->mp.y; - z=gf->mp.z; -} - - -GPoint gmshFace::point(double par1,double par2) const +GPoint gmshFace::point(double par1, double par2) const { - double pp[2]={par1,par2}; + double pp[2] = {par1, par2}; if(s->Typ == MSH_SURF_PLAN){ - double x,y,z,VX[3],VY[3]; - computePlaneDatas (this,VX,VY,x,y,z); - return GPoint( x + VX[0] * par1 + VY[0] * par2, - y + VX[1] * par1 + VY[1] * par2, - z + VX[2] * par1 + VY[2] * par2, this,pp); + double x, y, z, VX[3], VY[3]; + getMeanPlaneData(VX, VY, x, y, z); + return GPoint(x + VX[0] * par1 + VY[0] * par2, + y + VX[1] * par1 + VY[1] * par2, + z + VX[2] * par1 + VY[2] * par2, this, pp); } else{ - Vertex v = InterpolateSurface( s, par1, par2,0,0); - return GPoint(v.Pos.X,v.Pos.Y,v.Pos.Z,this,pp); + Vertex v = InterpolateSurface(s, par1, par2, 0, 0); + return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, pp); } } @@ -125,19 +151,19 @@ int gmshFace::containsParam(const SPoint2 &pt) const { Range<double> uu = parBounds(0); Range<double> vv = parBounds(1); - if ((pt.x()>=uu.low() && pt.x()<=uu.high()) && (pt.y()>=vv.low() && pt.y()<=vv.high())) - return 1; + if((pt.x() >= uu.low() && pt.x() <= uu.high()) && + (pt.y() >= vv.low() && pt.y() <= vv.high())) + return 1; else - return 0; + return 0; } - SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const { double u,v; if (s->Typ == MSH_SURF_PLAN){ double x,y,z,VX[3],VY[3]; - computePlaneDatas (this,VX,VY,x,y,z); + getMeanPlaneData(VX, VY, x, y, z); double vec[3] = {qp.x()-x,qp.y()-y,qp.z()-z}; prosca(vec,VX,&u); prosca(vec,VY,&v); @@ -186,5 +212,32 @@ void * gmshFace::getNativePtr() const int gmshFace::containsPoint(const SPoint3 &pt) const { - throw; + if(geomType() == Plane){ + // OK to use the normal from the mean plane here: we compensate + // for the (possibly wrong) orientation at the end + double n[3] = {meanPlane.a, meanPlane.b, meanPlane.c}; + double angle = 0.; + Vertex v; + v.Pos.X = pt.x(); + v.Pos.Y = pt.y(); + v.Pos.Z = pt.z(); + for(int i = 0; i < List_Nbr(s->Generatrices); i++) { + Curve *c; + List_Read(s->Generatrices, i, &c); + int N = (c->Typ == MSH_SEGM_LINE) ? 1 : 10; + for(int j = 0; j < N; j++) { + double u1 = (double)j / (double)N; + double u2 = (double)(j + 1) / (double)N; + Vertex p1 = InterpolateCurve(c, u1, 0); + Vertex p2 = InterpolateCurve(c, u2, 0); + angle += angle_plan(&v, &p1, &p2, n); + } + } + // we're inside if angle equals 2 * pi + if(fabs(angle) > 2 * M_PI - 0.5 && fabs(angle) < 2 * M_PI + 0.5) + return true; + return false; + } + else + throw; } diff --git a/Geo/gmshModel.cpp b/Geo/gmshModel.cpp index 94f5081595d4c8192c6b6adb8e9c711d122c9a49..a8879d805fdd1e81a9374dd82ae076cfe9b1deb7 100644 --- a/Geo/gmshModel.cpp +++ b/Geo/gmshModel.cpp @@ -17,11 +17,12 @@ void gmshModel::import() for(int i = 0; i < List_Nbr(points); i++){ Vertex *p; List_Read(points, i, &p); - if(!vertexByTag(p->Num)){ - gmshVertex *v = new gmshVertex(this, p); - v->setVisibility(p->Visible); + GVertex *v = vertexByTag(p->Num); + if(!v){ + v = new gmshVertex(this, p); add(v); } + v->setVisibility(p->Visible); } List_Delete(points); } @@ -31,23 +32,14 @@ void gmshModel::import() Curve *c; List_Read(curves, i, &c); if(c->Num >= 0 && c->beg && c->end){ - if(!vertexByTag(c->beg->Num)){ - gmshVertex *v = new gmshVertex(this, c->beg); - v->setVisibility(c->beg->Visible); - add(v); - } - if(!vertexByTag(c->end->Num)){ - gmshVertex *v = new gmshVertex(this, c->end); - v->setVisibility(c->beg->Visible); - add(v); - } - if(!edgeByTag(c->Num)){ - gmshEdge *e = new gmshEdge(this, c, - vertexByTag(c->beg->Num), - vertexByTag(c->end->Num)); - e->setVisibility(c->Visible); + GEdge *e = edgeByTag(c->Num); + if(!e){ + e = new gmshEdge(this, c, + vertexByTag(c->beg->Num), + vertexByTag(c->end->Num)); add(e); } + e->setVisibility(c->Visible); } } List_Delete(curves); @@ -57,11 +49,12 @@ void gmshModel::import() for(int i = 0; i < List_Nbr(surfaces); i++){ Surface *s; List_Read(surfaces, i, &s); - if(!faceByTag(s->Num)){ - gmshFace *f = new gmshFace(this, s); - f->setVisibility(s->Visible); + GFace *f = faceByTag(s->Num); + if(!f){ + f = new gmshFace(this, s); add(f); } + f->setVisibility(s->Visible); } List_Delete(surfaces); } @@ -70,11 +63,12 @@ void gmshModel::import() for(int i = 0; i < List_Nbr(volumes); i++){ Volume *v; List_Read(volumes, i, &v); - if(!regionByTag(v->Num)){ - gmshRegion *r = new gmshRegion(this, v); - r->setVisibility(v->Visible); + GRegion *r = regionByTag(v->Num); + if(!r){ + r = new gmshRegion(this, v); add(r); } + r->setVisibility(v->Visible); } List_Delete(volumes); } diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp index e0c6de00a1f9ffc3639126afab47f2c92250b9f3..ab641b3a36ffb1b42a12de7d443c7e6108486b7d 100644 --- a/Graphics/Geom.cpp +++ b/Graphics/Geom.cpp @@ -1,4 +1,4 @@ -// $Id: Geom.cpp,v 1.107 2006-08-12 21:31:24 geuzaine Exp $ +// $Id: Geom.cpp,v 1.108 2006-08-13 02:46:53 geuzaine Exp $ // // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // @@ -25,6 +25,7 @@ #include "Context.h" #include "gl2ps.h" #include "GModel.h" +#include "SBoundingBox3d.h" extern Context_T CTX; extern GModel *GMODEL; @@ -214,28 +215,112 @@ private: } if(CTX.geom.normals) { - const double eps = 1.e-3; - GPoint p1 = f->point(0.5, 0.5); - GPoint p2 = f->point(0.5 + eps, 0.5); - GPoint p3 = f->point(0.5, 0.5 + eps); - double n[3]; - normal3points(p1.x(), p1.y(), p1.z(), - p2.x(), p2.y(), p2.z(), - p3.x(), p3.y(), p3.z(), n); - n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0]; - n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1]; - n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2]; + SPoint2 p2 = SPoint2(0.5, 0.5); + SVector3 nn = f->normal(p2); + GPoint p = f->point(p2); + double n[3] = {nn.x() * CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0], + nn.y() * CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1], + nn.z() * CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2]}; glColor4ubv((GLubyte *) & CTX.color.geom.normals); Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius, - p1.x(), p1.y(), p1.z(), n[0], n[1], n[2], CTX.geom.light); + p.x(), p.y(), p.z(), n[0], n[1], n[2], CTX.geom.light); } } - + void _drawPlaneGFace(GFace *f) { - } + // we must lock this part to avoid it being called recursively + // when redraw events are fired in rapid succession + if(!f->cross.size() && !CTX.threads_lock) { + CTX.threads_lock = 1; + std::list<GVertex*> pts = f->vertices(); + if(pts.size()){ + SBoundingBox3d bb; + for(std::list<GVertex*>::iterator it = pts.begin(); it != pts.end(); it++){ + SPoint3 p((*it)->x(), (*it)->y(), (*it)->z()); + SPoint2 uv = f->parFromPoint(p); + bb += SPoint3(uv.x(), uv.y(), 0.); + } + GPoint v0 = f->point(bb.min().x(), bb.min().y()); + GPoint v1 = f->point(bb.max().x(), bb.min().y()); + GPoint v2 = f->point(bb.max().x(), bb.max().y()); + GPoint v3 = f->point(bb.min().x(), bb.max().y()); + const int N = 100; + for(int dir = 0; dir < 2; dir++) { + int end_line = 0; + SPoint3 pt, pt_last_good; + for(int i = 0; i < N; i++) { + double t = (double)i / (double)(N - 1); + double x, y, z; + if(dir){ + x = t * 0.5 * (v0.x() + v1.x()) + (1. - t) * 0.5 * (v2.x() + v3.x()); + y = t * 0.5 * (v0.y() + v1.y()) + (1. - t) * 0.5 * (v2.y() + v3.y()); + z = t * 0.5 * (v0.z() + v1.z()) + (1. - t) * 0.5 * (v2.z() + v3.z()); + } + else{ + x = t * 0.5 * (v0.x() + v3.x()) + (1. - t) * 0.5 * (v2.x() + v1.x()); + y = t * 0.5 * (v0.y() + v3.y()) + (1. - t) * 0.5 * (v2.y() + v1.y()); + z = t * 0.5 * (v0.z() + v3.z()) + (1. - t) * 0.5 * (v2.z() + v1.z()); + } + pt.setPosition(x, y, z); + if(f->containsPoint(pt)){ + pt_last_good.setPosition(pt.x(), pt.y(), pt.z()); + if(!end_line) { f->cross.push_back(pt); end_line = 1; } + } + else { + if(end_line) { f->cross.push_back(pt_last_good); end_line = 0; } + } + } + if(end_line) f->cross.push_back(pt_last_good); + } + } + // if we couldn't determine a cross, add dummy point so that + // we won't try again + if(!f->cross.size()) f->cross.push_back(SPoint3(0., 0., 0.)); + CTX.threads_lock = 0; + } + + if(f->cross.size() < 2) return ; + + if(CTX.geom.surfaces) { + glEnable(GL_LINE_STIPPLE); + glLineStipple(1, 0x1F1F); + gl2psEnable(GL2PS_LINE_STIPPLE); + glBegin(GL_LINES); + for(unsigned int i = 0; i < f->cross.size(); i++) + glVertex3d(f->cross[i].x(), f->cross[i].y(), f->cross[i].z()); + glEnd(); + glDisable(GL_LINE_STIPPLE); + gl2psDisable(GL2PS_LINE_STIPPLE); + } + + if(CTX.geom.surfaces_num) { + char Num[100]; + sprintf(Num, "%d", f->tag()); + double offset = 0.3 * CTX.gl_fontsize * CTX.pixel_equiv_x; + glRasterPos3d(0.5 * (f->cross[0].x() + f->cross[1].x()) + offset / CTX.s[0], + 0.5 * (f->cross[0].y() + f->cross[1].y()) + offset / CTX.s[0], + 0.5 * (f->cross[0].z() + f->cross[1].z()) + offset / CTX.s[0]); + Draw_String(Num); + } + if(CTX.geom.normals) { + SPoint3 p(0.5 * (f->cross[0].x() + f->cross[1].x()), + 0.5 * (f->cross[0].y() + f->cross[1].y()), + 0.5 * (f->cross[0].z() + f->cross[1].z())); + SPoint2 uv = f->parFromPoint(p); + SVector3 nn = f->normal(uv); + double n[3] = {nn.x() * CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0], + nn.y() * CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1], + nn.z() * CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2]}; + glColor4ubv((GLubyte *) & CTX.color.geom.normals); + Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, + CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius, + p.x(), p.y(), p.z(), n[0], n[1], n[2], CTX.geom.light); + } + } + public : void operator () (GFace *f) { diff --git a/Graphics/Makefile b/Graphics/Makefile index 2ef72d7401b34cb9257e1f518a0c769f82bc22c9..0185be82a2e35024e0e6bbba8af66656acb2cd02 100644 --- a/Graphics/Makefile +++ b/Graphics/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.94 2006-08-12 19:47:57 geuzaine Exp $ +# $Id: Makefile,v 1.95 2006-08-13 02:46:53 geuzaine Exp $ # # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle # @@ -114,17 +114,15 @@ Geom.o: Geom.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \ ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \ ../Common/SmoothNormals.h ../Numeric/Numeric.h ../Common/GmshMatrix.h \ ../Common/AdaptiveViews.h ../Common/GmshMatrix.h ../Common/Context.h \ - ../Plugin/Plugin.h ../Common/Options.h ../Plugin/PluginManager.h \ - ../Plugin/Plugin.h gl2ps.h ../Geo/GModel.h ../Geo/GVertex.h \ - ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \ - ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/MVertex.h \ - ../Geo/SPoint3.h ../Common/GmshDefines.h ../Geo/MVertex.h \ - ../Geo/GPoint.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \ - ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \ - ../Geo/MElement.h ../Geo/MVertex.h ../Geo/GFace.h ../Geo/GPoint.h \ - ../Geo/GEntity.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \ - ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \ - ../Geo/SBoundingBox3d.h + gl2ps.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h \ + ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \ + ../Geo/SPoint3.h ../Geo/MVertex.h ../Geo/SPoint3.h \ + ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/GPoint.h ../Geo/GEdge.h \ + ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \ + ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/MVertex.h \ + ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/MElement.h \ + ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h \ + ../Geo/GEntity.h ../Geo/MElement.h ../Geo/SBoundingBox3d.h # 1 "/Users/geuzaine/.gmsh/Graphics//" Post.o: Post.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \ ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \ diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp index 9e5acaeb4a0c517cf5f5cf8b3aa9f3a1cb3e2915..13a42b1e4d3fdce999f6a368228cdf5f8fe988ec 100644 --- a/Mesh/meshGEdge.cpp +++ b/Mesh/meshGEdge.cpp @@ -106,12 +106,12 @@ void meshGEdge :: operator() (GEdge *ge) { if(ge->geomType() == GEntity::DiscreteCurve) return; + // Send a messsage to the GMSH environment + Msg(INFO, "Meshing curve %d", ge->tag()); + deMeshGEdge dem; dem(ge); - // Send a messsage to the GMSH environment - Msg(STATUS2, "Meshing curve %d", ge->tag()); - // Create a list of integration points List_T *Points = List_Create(10, 10, sizeof(IntPoint)); diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index f4c2a9fa89ea9ec0b2c35b1b48695143154f4cbe..65cf56073d1909270a1ef5896d4370c863f1283b 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -12,22 +12,6 @@ #include "Numeric.h" #include "BDS.h" -#if defined(HAVE_GSL) -#include <gsl/gsl_vector.h> -#include <gsl/gsl_linalg.h> -#else -#define NRANSI -#include "nrutil.h" -void dsvdcmp(double **a, int m, int n, double w[], double **v); -#endif - -/* - Le concept d'un plan moyen calcule au sens des moidres carres n'est - pas le bon pour les surfaces non-planes : imagine un quart de cercle - extrude d'une faible hauteur. Le plan moyen sera dans le plan du - cercle! En attendant mieux, il y a donc un test de coherence - supplementaire pour les surfaces non-planes. */ - int Orientation (std::vector<MVertex*> &cu) { int N, i, a, b, c; @@ -167,185 +151,6 @@ void computeEdgeLoops (const GFace *gf, } -void MeanPlane_bis(GFace *gf, const std::vector<MVertex*> &points) -{ - - int min, ndata, na; - double res[4], ex[3], t1[3], t2[3], svd[3]; - double xm = 0., ym = 0., zm = 0.; - - ndata = points.size(); - na = 3; - for(int i = 0; i < ndata; i++) { - xm += points[i]->x(); - ym += points[i]->y(); - zm += points[i]->z(); - } - xm /= (double)ndata; - ym /= (double)ndata; - zm /= (double)ndata; - -#if defined(HAVE_GSL) - gsl_matrix *U = gsl_matrix_alloc(ndata, na); - gsl_matrix *V = gsl_matrix_alloc(na, na); - gsl_vector *W = gsl_vector_alloc(na); - gsl_vector *TMPVEC = gsl_vector_alloc(na); - for(int i = 0; i < ndata; i++) { - gsl_matrix_set(U, i, 0, points[i]->x() - xm); - gsl_matrix_set(U, i, 1, points[i]->y() - ym); - gsl_matrix_set(U, i, 2, points[i]->z() - zm); - } - gsl_linalg_SV_decomp(U, V, W, TMPVEC); - svd[0] = gsl_vector_get(W, 0); - svd[1] = gsl_vector_get(W, 1); - svd[2] = gsl_vector_get(W, 2); - if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2])) - min = 0; - else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2])) - min = 1; - else - min = 2; - res[0] = gsl_matrix_get(V, 0, min); - res[1] = gsl_matrix_get(V, 1, min); - res[2] = gsl_matrix_get(V, 2, min); - norme(res); - gsl_matrix_free(U); - gsl_matrix_free(V); - gsl_vector_free(W); - gsl_vector_free(TMPVEC); -#else - double **U = dmatrix(1, ndata, 1, na); - double **V = dmatrix(1, na, 1, na); - double *W = dvector(1, na); - for(int i = 0; i < ndata; i++) { - U[i + 1][1] = points[i]->x() - xm; - U[i + 1][2] = points[i]->y() - ym; - U[i + 1][3] = points[i]->z() - zm; - } - dsvdcmp(U, ndata, na, W, V); - if(fabs(W[1]) < fabs(W[2]) && fabs(W[1]) < fabs(W[3])) - min = 1; - else if(fabs(W[2]) < fabs(W[1]) && fabs(W[2]) < fabs(W[3])) - min = 2; - else - min = 3; - svd[0] = W[1]; - svd[1] = W[2]; - svd[2] = W[3]; - res[0] = V[1][min]; - res[1] = V[2][min]; - res[2] = V[3][min]; - norme(res); - free_dmatrix(U, 1, ndata, 1, na); - free_dmatrix(V, 1, na, 1, na); - free_dvector(W, 1, na); -#endif - - // check coherence of results for non-plane surfaces - if(gf->geomType() != GEntity::Plane) { - double res2[3], c[3], cosc, sinc, angplan; - double eps = 1.e-3; - - GPoint v1 = gf->point( 0.5, 0.5); - GPoint v2 = gf->point( 0.5 + eps, 0.5); - GPoint v3 = gf->point( 0.5, 0.5 + eps); - t1[0] = v2.x() - v1.x(); - t1[1] = v2.y() - v1.y(); - t1[2] = v2.z() - v1.z(); - t2[0] = v3.x() - v1.x(); - t2[1] = v3.y() - v1.y(); - t2[2] = v3.z() - v1.z(); - norme(t1); - norme(t2); - // prodve(t1, t2, res2); - // Warning: the rest of the code assumes res = t2 x t1, not t1 x t2 (WTF?) - prodve(t2, t1, res2); - norme(res2); - prodve(res, res2, c); - prosca(res, res2, &cosc); - sinc = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]); - angplan = myatan2(sinc, cosc); - angplan = angle_02pi(angplan) * 180. / Pi; - if((angplan > 70 && angplan < 110) || (angplan > 260 && angplan < 280)) { - Msg(INFO, "SVD failed (angle=%g): using rough algo...", angplan); - res[0] = res2[0]; - res[1] = res2[1]; - res[2] = res2[2]; - goto end; - } - } - - ex[0] = ex[1] = ex[2] = 0.0; - if(res[0] == 0.) - ex[0] = 1.0; - else if(res[1] == 0.) - ex[1] = 1.0; - else - ex[2] = 1.0; - - prodve(res, ex, t1); - norme(t1); - prodve(t1, res, t2); - norme(t2); - -end: - res[3] = (xm * res[0] + ym * res[1] + zm * res[2]); - - for(int i = 0; i < 3; i++) - gf->mp.plan[0][i] = t1[i]; - for(int i = 0; i < 3; i++) - gf->mp.plan[1][i] = t2[i]; - for(int i = 0; i < 3; i++) - gf->mp.plan[2][i] = res[i]; - - gf->mp.a = res[0]; - gf->mp.b = res[1]; - gf->mp.c = res[2]; - gf->mp.d = res[3]; - - gf->mp.x=gf->mp.y=gf->mp.z=0; - if (fabs(gf->mp.a) >= fabs(gf->mp.b) && fabs(gf->mp.a) >= fabs(gf->mp.c) ) - { - gf->mp.x = gf->mp.d/gf->mp.a; - } - else if (fabs(gf->mp.b) >= fabs(gf->mp.a) && fabs(gf->mp.b) >= fabs(gf->mp.c) ) - { - gf->mp.y = gf->mp.d/gf->mp.b; - } - else - { - gf->mp.z = gf->mp.d/gf->mp.c; - } - - - Msg(DEBUG1, "Surface: %d", gf->tag()); - Msg(DEBUG2, "SVD : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min); - Msg(DEBUG2, "Plane : (%g x + %g y + %g z = %g)", gf->mp.a, gf->mp.b, gf->mp.c, gf->mp.d); - Msg(DEBUG2, "Normal : (%g , %g , %g )", gf->mp.a, gf->mp.b, gf->mp.c); - Msg(DEBUG3, "t1 : (%g , %g , %g )", t1[0], t1[1], t1[2]); - Msg(DEBUG3, "t2 : (%g , %g , %g )", t2[0], t2[1], t2[2]); - Msg(DEBUG3, "pt : (%g , %g , %g )", gf->mp.x, gf->mp.y, gf->mp.z); - - //check coherence for plane surfaces - if(gf->geomType() == GEntity::Plane) { - std::list<GVertex*> verts = gf->vertices(); - std::list<GVertex*>::const_iterator itv = verts.begin(); - - for (;itv!=verts.end();itv++) - { - const GVertex *v = *itv; - double d = - gf->mp.a * v->x() + gf->mp.b * v->y() + gf->mp.c * v->z() - gf->mp.d; - if(fabs(d) > CTX.lc * 1.e-3) { - Msg(GERROR1, "Plane surface %d (%gx+%gy+%gz+%g=0) is not plane!", - v->tag(), gf->mp.a, gf->mp.b, gf->mp.c, gf->mp.d); - Msg(GERROR3, "Control point %d = (%g,%g,%g), val=%g", - v->tag(), v->x(), v->y(), v->z(), d); - return; - } - } - } -} void computeEdgeParameters ( double x1, double y1, double x2, double y2, GFace *gf , const int numberOfTestPoints, double &coordMiddle, double &edgeLength ) @@ -759,6 +564,9 @@ void meshGFace :: operator() (GFace *gf) { if(gf->geomType() == GEntity::DiscreteSurface) return; + // Send a messsage to the GMSH environment + Msg(INFO, "Meshing surface %d", gf->tag()); + // destroy the mesh if it exists deMeshGFace dem; dem(gf); @@ -771,10 +579,10 @@ void meshGFace :: operator() (GFace *gf) // compute loops on the fly // indices indicate start and end points of a loop // loops are not yet oriented - computeEdgeLoops(gf,points,indices); + computeEdgeLoops(gf, points, indices); - //compute the mean plane, this is sometimes useful - MeanPlane_bis(gf,points); + // compute the mean plane, this is sometimes useful + gf->computeMeanPlane(points); Msg(DEBUG1, "Face %d type %d with %d edge loops and %d points", gf->tag(),gf->geomType(),indices.size()-1,points.size()); diff --git a/Parser/Makefile b/Parser/Makefile index 4c1cd3a0457254b9aae3aa69f0437e26c50e19c4..103ddf43fed258713d49520b5cac10c5c2ff7fa6 100644 --- a/Parser/Makefile +++ b/Parser/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.95 2006-08-12 19:34:16 geuzaine Exp $ +# $Id: Makefile,v 1.96 2006-08-13 02:46:53 geuzaine Exp $ # # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle # @@ -120,8 +120,7 @@ OpenFile.o: OpenFile.cpp ../Common/Gmsh.h ../Common/Message.h \ ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \ ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \ ../Geo/SBoundingBox3d.h ../Common/Context.h Parser.h OpenFile.h \ - ../Common/CommandLine.h ../Plugin/PluginManager.h ../Plugin/Plugin.h \ - ../Common/Options.h ../Common/Views.h ../Common/ColorTable.h \ + ../Common/CommandLine.h ../Common/Views.h ../Common/ColorTable.h \ ../Common/VertexArray.h ../Common/SmoothNormals.h \ ../Common/GmshMatrix.h ../Common/AdaptiveViews.h ../Common/GmshMatrix.h \ ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Vertex.h \ diff --git a/Plugin/Makefile b/Plugin/Makefile index eefdf276b6a8b17874ec73d5e53094ffd789e635..9e7944f34da3d9f4339bd69b0cd19def8ead7591 100644 --- a/Plugin/Makefile +++ b/Plugin/Makefile @@ -1,4 +1,4 @@ -# $Id: Makefile,v 1.105 2006-08-12 21:31:24 geuzaine Exp $ +# $Id: Makefile,v 1.106 2006-08-13 02:46:54 geuzaine Exp $ # # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle # @@ -79,12 +79,6 @@ Plugin.o: Plugin.cpp Plugin.h ../Common/Options.h ../Common/Message.h \ Extract.h ExtractElements.h ExtractEdges.h HarmonicToTime.h \ ModulusPhase.h Integrate.h Gradient.h Curl.h Divergence.h Annotate.h \ Remove.h DecomposeInSimplex.h Smooth.h Transform.h Triangulate.h Warp.h \ - ../Geo/Geo.h ../Mesh/Mesh.h ../Common/GmshDefines.h \ - ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h \ - ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h \ - ../Mesh/Element.h ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \ - ../Geo/ExtrudeParams.h ../Mesh/Metric.h ../Mesh/Vertex.h \ - ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h ../Common/GmshUI.h \ Eigenvectors.h Eigenvalues.h Lambda2.h Evaluate.h \ ../Common/OctreePost.h ../Common/Octree.h ../Common/OctreeInternals.h \ Probe.h ../Common/Context.h