diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 058239095053ecfa780cbff1f1f8b2f451418855..f80e5d58b72aa923ea5ff86333a1166c32560b75 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.17 2006-11-14 17:11:33 remacle Exp $
+// $Id: GFace.cpp,v 1.18 2006-11-14 20:20:18 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -23,6 +23,7 @@
#include "GFace.h"
#include "GEdge.h"
#include "Message.h"
+#include "Utils.h"
#if defined(HAVE_GSL)
#include <gsl/gsl_vector.h>
@@ -393,3 +394,90 @@ void GFace :: computeDirs ()
{
throw;
}
+
+void GFace::XYZtoUV(const double X, const double Y, const double Z,
+ double &U, double &V,
+ const double relax) const
+{
+ const double Precision = 1.e-8;
+ const int MaxIter = 25;
+ const int NumInitGuess = 11;
+
+ 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);
+ 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];
+ err = 1.0;
+ iter = 1;
+
+ while(err > Precision && iter < MaxIter) {
+ 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();
+ mat[1][0] = der.right().x();
+ mat[1][1] = der.right().y();
+ mat[1][2] = der.right().z();
+ mat[2][0] = 0.;
+ mat[2][1] = 0.;
+ mat[2][2] = 0.;
+ invert_singular_matrix3x3(mat, jac);
+
+ Unew = U + relax *
+ (jac[0][0] * (X - P.x()) + jac[1][0] * (Y - P.y()) +
+ jac[2][0] * (Z - P.z()));
+ Vnew = V + relax *
+ (jac[0][1] * (X - P.x()) + jac[1][1] * (Y - P.y()) +
+ jac[2][1] * (Z - P.z()));
+
+ err = DSQR(Unew - U) + DSQR(Vnew - V);
+ // A BETTER TEST !! (JFR/AUG 2006)
+ err2 = DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z());
+ iter++;
+ U = Unew;
+ V = Vnew;
+ }
+
+
+ if(iter < MaxIter && err <= Precision && err2 <= 1.e-5 &&
+ Unew <= umax && Vnew <= vmax &&
+ Unew >= umin && Vnew >= vmin){
+ if (err2 > Precision)
+ Msg(WARNING,"converged for i=%d j=%d (err=%g iter=%d) BUT err2 = %g", i, j, err, iter,err2);
+ return;
+ }
+ }
+ }
+
+ 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);
+ }
+}
+
+
+SPoint2 GFace::parFromPoint(const SPoint3 &p) const
+{
+ double U,V;
+
+ XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
+
+ return SPoint2(U,V);
+}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 233b5913a5fb70fcddb025442fcdc6c187dce6c7..53037481eaea393a3cc2a4e052522671db3f7266 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -68,6 +68,11 @@ class GFace : public GEntity
virtual int dim() const {return 2;}
virtual void setVisibility(char val, bool recursive=false);
+ // compute XYZ from parametric UV
+ void XYZtoUV(const double X, const double Y, const double Z,
+ double &U, double &V,
+ const double relax) const;
+
// The bounding box
virtual SBoundingBox3d bounds() const;
@@ -81,7 +86,7 @@ class GFace : public GEntity
// Return the parmater location on the face given a point in space
// that is on the face.
- virtual SPoint2 parFromPoint(const SPoint3 &) const = 0;
+ virtual SPoint2 parFromPoint(const SPoint3 &) const;
// True if the parameter value is interior to the face.
virtual int containsParam(const SPoint2 &pt) const = 0;
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 52213ea47279531df7f7d79711c95e1a943cfdb9..16b5003d2a019b2f65c355c85c7074d670c9c5b6 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.1 2006-11-14 17:11:33 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.2 2006-11-14 20:20:18 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -75,7 +75,7 @@ GEntity::GeomType OCCEdge::geomType() const
int OCCEdge::minimumMeshSegments () const
{
- return GEdge::minimumMeshSegments () ;
+ return CTX.mesh.min_circ_points; ;
}
int OCCEdge::minimumDrawSegments () const
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 475d1cc5645bda855cd90e6ae4d4e56d994a3e1a..fee6ed52ffcaa9770740fc09621797463499d218 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.1 2006-11-14 17:11:33 remacle Exp $
+// $Id: OCCFace.cpp,v 1.2 2006-11-14 20:20:18 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -91,14 +91,23 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
Msg(INFO,"OCC Face %d with %d edges",num,l_edges.size());
ShapeAnalysis::GetFaceUVBounds (s, umin, umax, vmin, vmax);
+ // we do that for the projections to converge on the
+ // borders of the surface
+ umin -= fabs(umax-umin)/100.0;
+ vmin -= fabs(vmax-vmin)/100.0;
+ umax += fabs(umax-umin)/100.0;
+ vmax += fabs(vmax-vmin)/100.0;
occface = BRep_Tool::Surface(s);
}
Range<double> OCCFace::parBounds(int i) const
-{
+{
+ double umin2, umax2, vmin2, vmax2;
+
+ ShapeAnalysis::GetFaceUVBounds (s, umin2, umax2, vmin2, vmax2);
if (i==0)
- return Range<double>(umin, umax);
- return Range<double>(vmin, vmax);
+ return Range<double>(umin2, umax2);
+ return Range<double>(vmin2, vmax2);
}
SVector3 OCCFace::normal(const SPoint2 ¶m) const
@@ -164,9 +173,17 @@ int OCCFace::containsParam(const SPoint2 &pt) const
SPoint2 OCCFace::parFromPoint(const SPoint3 &qp) const
{
-
- GPoint gp = closestPoint(qp);
- return SPoint2(gp.u(), gp.v());
+ 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");
+ return GFace::parFromPoint(qp);
+ }
+ pnt = proj.NearestPoint();
+ double pp[2];
+ proj.LowerDistanceParameters (pp[0], pp[1]);
+ return SPoint2(pp[0],pp[1]);
}
GEntity::GeomType OCCFace::geomType() const
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index 8b0b1c9016a98597664127668432f0fcda3f8af5..3b74e43a5ca64bc41397dd4e4f2b1833b31a6b33 100644
--- a/Geo/OCCVertex.h
+++ b/Geo/OCCVertex.h
@@ -25,6 +25,7 @@
#include "GModel.h"
#include "OCCIncludes.h"
#include "GVertex.h"
+#include "Context.h"
class OCCVertex : public GVertex {
protected:
@@ -56,7 +57,11 @@ class OCCVertex : public GVertex {
return pnt.Z();
}
void * getNativePtr() const { return (void*) &v; }
- virtual double prescribedMeshSizeAtVertex() const { return 350; }
+ virtual double prescribedMeshSizeAtVertex() const {
+ SBoundingBox3d b = model()->bounds();
+ double lc = norm ( SVector3 ( b.max() , b.min() ) );
+ return lc*CTX.mesh.lc_factor;
+ }
};
#endif
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 64fcaba8db4a1bff8cbd219bcff2f730213647b5..4f41bf06cb7d4e92c9ac718d08103a3f25bd1d18 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.20 2006-11-14 17:11:33 remacle Exp $
+// $Id: gmshFace.cpp,v 1.21 2006-11-14 20:20:18 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -201,18 +201,19 @@ int gmshFace::containsParam(const SPoint2 &pt) const
SPoint2 gmshFace::parFromPoint(const SPoint3 &qp) const
{
- double u,v;
if(s->Typ == MSH_SURF_PLAN){
+ double u,v;
double x,y,z,VX[3],VY[3];
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);
+ return SPoint2(u, v);
}
else{
- XYZtoUV(s, qp.x(), qp.y(), qp.z(), &u, &v, 1.0);
+ // XYZtoUV(s, qp.x(), qp.y(), qp.z(), &u, &v, 1.0);
+ return GFace::parFromPoint(qp);
}
- return SPoint2(u, v);
}
GEntity::GeomType gmshFace::geomType() const
diff --git a/Mesh/Utils.h b/Mesh/Utils.h
index c3eeff146809268251f7609dfc29f976e4eef45f..31ea74d1cd9225dd5facbbfed553f995d23e2703 100644
--- a/Mesh/Utils.h
+++ b/Mesh/Utils.h
@@ -24,6 +24,7 @@
#include "Vertex.h"
#include "Mesh.h"
+void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
void direction (Vertex * v1, Vertex * v2, double d[3]);
void Projette (Vertex * v, double mat[3][3]);
void MeanPlane(List_T *point, Surface *s);
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index ac128c88a48d87b55e2e05da4001c67580c6314a..775b7271899e120e46c482373471cff6cad77fdc 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.125 2006-11-14 17:11:33 remacle Exp $
+// $Id: OpenFile.cpp,v 1.126 2006-11-14 20:20:18 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -293,7 +293,8 @@ int MergeProblem(char *name, int warn_if_missing)
if(!strcmp(ext, ".stl") || !strcmp(ext, ".STL")){
status = GMODEL->readSTL(name, CTX.mesh.stl_distance_tol);
}
- else if(!strcmp(ext, ".iges") || !strcmp(ext, ".IGES")){
+ else if(!strcmp(ext, ".iges") || !strcmp(ext, ".IGES") ||
+ !strcmp(ext, ".igs") || !strcmp(ext, ".IGS")){
GMODEL->readOCCIGES(std::string(name));
}
else if(!strcmp(ext, ".step") || !strcmp(ext, ".STEP")){