Using BFGS to compute parametric coordinates

Common/OpenFile.cpp

@@ -362,7 +362,7 @@ int MergeFile(std::string fileName, bool warnIfMissing)
         GModel* tmp2 = GModel::current();
         GModel* tmp = new GModel();
         tmp->readMSH(fileName);
-        GeomMeshMatcher::instance()->match(tmp2, tmp);
+        status = GeomMeshMatcher::instance()->match(tmp2, tmp);
         delete tmp;
       } else
 	status = GModel::current()->readMSH(fileName);

Geo/GFace.cpp

@@ -22,8 +22,27 @@
 #include "meshGFaceOptimize.h"
 #include "meshGFaceLloyd.h"
 
+#if defined(HAVE_BFGS)
+#include "stdafx.h"
+#include "optimization.h"
+#endif
+
 #define SQU(a)      ((a)*(a))
 
+class data_wrapper{
+ private:
+  const GFace* gf;
+  SPoint3 point;
+ public:
+  data_wrapper() {gf = NULL; point = SPoint3();}
+  ~data_wrapper() {}
+  const GFace* get_face() {return gf;}
+  void set_face(const GFace* face) {gf = face;}
+  SPoint3 get_point() {return point;}
+  void set_point(SPoint3 _point) {point = SPoint3(_point);}
+};
+
+
 GFace::GFace(GModel *model, int tag)
   : GEntity(model, tag), r1(0), r2(0), compound(0), va_geom_triangles(0)
 {
@@ -906,11 +925,113 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p, bool onSurface) const
   return SPoint2(U, V);
 }
 
+#if defined(QUASINEWTON)
+// Callback function for BFGS
+void bfgs_callback(const alglib::real_1d_array& x, double& func, alglib::real_1d_array& grad, void* ptr) {
+  data_wrapper* w = static_cast<data_wrapper*>(ptr);
+  SPoint3 p = w->get_point();
+  const GFace* gf = w->get_face();
+
+  // Value of the objective
+  GPoint pnt = gf->point(x[0], x[1]);
+  SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
+  func = p.distance(spnt);
+  //printf("func : %f\n", func);
+
+  // Value of the gradient
+  std::vector<double> values(2);
+  Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(x[0], x[1]));
+  grad[0] = -2./(2.0*func) * (der.left().x() + der.left().y() + der.left().z());
+  grad[1] = -2./(2.0*func) * (der.right().x() + der.right().y() + der.right().z());
+  //printf("Gradients %f %f\n", grad[0], grad[1]);
+}
+#endif
+
 GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
 {
+#if defined(QUASINEWTON)
+
+  //
+  // Creating the optimisation problem
+  //
+  alglib::ae_int_t dim = 2;
+  alglib::ae_int_t corr = 2; // Num of corrections in the scheme in [3,7] 
+  alglib::minlbfgsstate state;
+  alglib::real_1d_array x; // = "[0.5,0.5]";
+  std::vector<double> initial_conditions(dim);
+
+  double min_dist = 1e18;
+  double min_u = 0.5;
+  double min_v = 0.5;
+
+  Range<double> uu = parBounds(0);
+  Range<double> vv = parBounds(1);
+
+  double initial_guesses = 1000.0;
+  for(double u = uu.low(); u<=uu.high(); u+=(uu.high()-uu.low())/initial_guesses) {
+    printf("%f\n", u);
+    for(double v = vv.low(); v<=vv.high(); v+=(vv.high()-vv.low())/initial_guesses) {
+      GPoint pnt = point(u, v);
+      printf("(point) : %f %f %f\n", pnt.x(), pnt.y(), pnt.z());
+      SPoint3 spnt(pnt.x(), pnt.y(), pnt.z());
+      double dist = queryPoint.distance(spnt);
+      if (dist < min_dist) {
+	printf("min_dist %f\n", dist);
+	min_dist = dist;
+	min_u = u;
+	min_v = v;
+	GPoint pnt = point(min_u, min_v);
+      }
+    }
+  }
+
+  initial_conditions[0] = min_u;
+  initial_conditions[1] = min_v;
+
+  printf("Initial conditions : %f %f\n", min_u, min_v);
+  GPoint pnt = point(min_u, min_v);
+  printf("Initial conditions (point) : %f %f %f\n", pnt.x(), pnt.y(), pnt.z());
+
+  x.setcontent(dim, &initial_conditions[0]);
+
+
+
+
+  minlbfgscreate(2, corr, x, state);
+
+  //
+  // Defining the stopping criterion
+  //
+  double epsg = 0.0;
+  double epsf = 0.0;
+  double epsx = 0.0;
+  alglib::ae_int_t maxits = 500;
+
+  minlbfgssetcond(state,epsg,epsf,epsx,maxits);
+
+  //
+  // Solving the problem
+  //
+
+  data_wrapper w;
+  w.set_point(queryPoint);
+  w.set_face(this);
+
+  minlbfgsoptimize(state,bfgs_callback,NULL,&w);
+
+  //
+  // Getting the results
+  //
+  alglib::minlbfgsreport rep;
+
+  minlbfgsresults(state,x,rep);
+
+  return point(x[0],x[1]);
+#else
   Msg::Error("Closest point not implemented for this type of surface");
   SPoint2 p = parFromPoint(queryPoint, false);
   return point(p);
+#endif
 }
 
 bool GFace::containsParam(const SPoint2 &pt) const

Geo/GeoInterpolation.cpp

@@ -295,6 +295,9 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
   switch (c->Typ) {
 
   case MSH_SEGM_LINE:
+#if defined(QUASINEWTON)
+    printf("MSH_SEGM_LINE\n");
+#endif
     N = List_Nbr(c->Control_Points);
     i = (int)((double)(N - 1) * u);
     while(i >= N - 1)
@@ -415,6 +418,10 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
     break;
 
   case MSH_SEGM_DISCRETE:
+#if defined(QUASINEWTON)
+    printf("MSH_SEGM_DISCRETE\n");
+#endif
+
     Msg::Debug("Cannot interpolate discrete curve");
     break;
 
@@ -448,6 +455,9 @@ static Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
                      s1.Pos.Y, s2.Pos.Y, s3.Pos.Y, s4.Pos.Y, u, v);
   V.Pos.Z = TRAN_QUA(c1.Pos.Z, c2.Pos.Z, c3.Pos.Z, c4.Pos.Z,
                      s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, s4.Pos.Z, u, v);
+  //printf("TRANQUA %f %f %f %f %f %f %f %f\n", c1.Pos.Z, c2.Pos.Z, c3.Pos.Z, c4.Pos.Z,
+  //                   s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, s4.Pos.Z);
+  //printf("V.Pos.Z %f\n", V.Pos.Z);
   return (V);
 }
 
@@ -592,6 +602,17 @@ static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
     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 defined(QUASINEWTON)
+    printf("----- u %f v %f\n", u, v);
+    printf("S[0] %f %f %f\n", S[0]->Pos.X, S[0]->Pos.Y,S[0]->Pos.Z);
+    printf("S[1] %f %f %f\n", S[1]->Pos.X, S[1]->Pos.Y,S[1]->Pos.Z);
+    printf("S[2] %f %f %f\n", S[2]->Pos.X, S[2]->Pos.Y,S[2]->Pos.Z);
+    printf("S[3] %f %f %f\n", S[3]->Pos.X, S[3]->Pos.Y,S[3]->Pos.Z);
+    printf("V[0] %f %f %f\n", V[0].Pos.X, V[0].Pos.Y,V[0].Pos.Z);
+    printf("V[1] %f %f %f\n", V[1].Pos.X, V[1].Pos.Y,V[1].Pos.Z);
+    printf("V[2] %f %f %f\n", V[2].Pos.X, V[2].Pos.Y,V[2].Pos.Z);
+    printf("V[3] %f %f %f\n", V[3].Pos.X, V[3].Pos.Y,V[3].Pos.Z);
+#endif
     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], *O, &T);
   }

Geo/GeomMeshMatcher.cpp

@@ -155,7 +155,7 @@ GeomMeshMatcher::matchVertices(GModel* m1, GModel *m2, bool& ok)
     //m2->add(*vert);
 
   Msg::Info("Vertices matched : %i / %i (%i)\n", num_matched_vertices, num_total_vertices, num_mesh_vertices);
-  if(num_matched_vertices != num_total_vertices) ok = false;
+  //if(num_matched_vertices != num_total_vertices) ok = false;
   return (coresp_v);
 }
 
@@ -225,16 +225,17 @@ GeomMeshMatcher::matchEdges(GModel* m1, GModel* m2,
     choice->setTag(((GEdge*) *entity1)->tag());
 
     // This reverses the edge if it's not parametrized in the right direction
+/*
     if (choice->getBeginVertex() == findMatching<GVertex*>(*coresp_v,v2) &&
         choice->getEndVertex() == findMatching<GVertex*>(*coresp_v,v1)) {
       Msg::Info("Wrong parametrization direction, reversing.");
       choice->reverse();
     }
-
-    for (unsigned int v = 0; v < ((GEdge*) choice)->getNumMeshVertices(); v++) {
+*/
+    /*for (unsigned int v = 0; v < ((GEdge*) choice)->getNumMeshVertices(); v++) {
       if (((GEdge*) choice)->getMeshVertex(v)->onWhat()->dim() > 0)
 	((GEdge*) choice)->getMeshVertex(v)->setEntity((GEdge*) *entity1);
-    }
+    }*/
     num_matched_edges++;
   }
 
@@ -301,10 +302,10 @@ GeomMeshMatcher:: matchFaces(GModel* m1, GModel* m2,
     coresp_f->push_back(Pair<GFace*,GFace*>((GFace*) *entity1 ,
                                              choice));
     choice->setTag(((GFace*) *entity1)->tag());
-    for (unsigned int v = 0; v < ((GFace*) choice)->getNumMeshVertices(); v++) {
+    /*for (unsigned int v = 0; v < ((GFace*) choice)->getNumMeshVertices(); v++) {
       if(((GFace*) choice)->getMeshVertex(v)->onWhat()->dim() > 1)
 	((GFace*) choice)->getMeshVertex(v)->setEntity((GFace*) *entity1);
-    }
+    }*/
 
     num_matched_faces++;
   }
@@ -400,10 +401,10 @@ GeomMeshMatcher::matchRegions(GModel* m1, GModel* m2,
                                              choice));
        choice->setTag(((GRegion*) *entity1)->tag());
 
-    for (unsigned int v = 0; v < ((GRegion*) choice)->getNumMeshVertices(); v++) {
-      if ( ((GRegion*) choice)->getMeshVertex(v)->onWhat()->dim() > 2)
-	((GRegion*) choice)->getMeshVertex(v)->setEntity((GRegion*) *entity1);
-    }
+    //for (unsigned int v = 0; v < ((GRegion*) choice)->getNumMeshVertices(); v++) {
+    //  if ( ((GRegion*) choice)->getMeshVertex(v)->onWhat()->dim() > 2)
+    //    ((GRegion*) choice)->getMeshVertex(v)->setEntity((GRegion*) *entity1);
+    //}
        num_matched_regions++;
     }
   }
@@ -597,6 +598,7 @@ int GeomMeshMatcher::forceTomatch(GModel *geom, GModel *mesh, const double TOL)
 
 static void copy_vertices (GVertex *to, GVertex *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
   if (from) {
+    to->deleteMesh();
     for (int i=0;i<from->mesh_vertices.size();i++){
       MVertex *v_from = from->mesh_vertices[i];
       MVertex *v_to = new MVertex (v_from->x(),v_from->y(),v_from->z(), to);
@@ -614,19 +616,25 @@ static void copy_vertices (GRegion *to, GRegion *from, std::map<MVertex*,MVertex
   }
 }
 static void copy_vertices (GEdge *to, GEdge *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+  to->deleteMesh();
   for (int i=0;i<from->mesh_vertices.size();i++){
     MVertex *v_from = from->mesh_vertices[i];
-    double t = to->parFromPoint ( SPoint3(v_from->x(),v_from->y(),v_from->z() ) );
+    double t;
+    GPoint gp = to->closestPoint(SPoint3(v_from->x(),v_from->y(),v_from->z()), t );
     MEdgeVertex *v_to = new MEdgeVertex (v_from->x(),v_from->y(),v_from->z(), to, t );
     to->mesh_vertices.push_back(v_to);
     _mesh_to_geom[v_from] = v_to;
   }
 }
 static void copy_vertices (GFace *to, GFace *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+  printf("Starting Face %d, with %d vertices\n", to->tag(),  from->mesh_vertices.size());
   for (int i=0;i<from->mesh_vertices.size();i++){
     MVertex *v_from = from->mesh_vertices[i];
-    SPoint2 uv = to->parFromPoint ( SPoint3(v_from->x(),v_from->y(),v_from->z() ) );
-    MFaceVertex *v_to = new MFaceVertex (v_from->x(),v_from->y(),v_from->z(), to, uv.x(),uv.y() );
+    double uv[2];
+    GPoint gp = to->closestPoint ( SPoint3(v_from->x(),v_from->y(),v_from->z()), uv );
+    printf("Original point %f %f %f\n", v_from->x(), v_from->y(), v_from->z());
+    printf("New point %f %f %f\n", gp.x(), gp.y(), gp.z());
+    MFaceVertex *v_to = new MFaceVertex (v_from->x(),v_from->y(),v_from->z(), to, uv[0],uv[1] );
     to->mesh_vertices.push_back(v_to);
     _mesh_to_geom[v_from] = v_to;
   }
@@ -635,12 +643,14 @@ static void copy_vertices (GFace *to, GFace *from, std::map<MVertex*,MVertex*> &
 template <class ELEMENT>
 static void copy_elements (std::vector<ELEMENT*> &to, std::vector<ELEMENT*> &from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
   MElementFactory toto;
-  for (int i=0;i < to.size();i++){
+  to.clear();
+  for (int i=0;i < from.size();i++){
     ELEMENT *e = from[i];
     std::vector<MVertex*> nodes;
-    for(int j=0;j<e->getNumVertices();j++)
+    for(int j=0;j<e->getNumVertices();j++) {
       nodes.push_back(_mesh_to_geom[e->getVertex(j)]);
-    to.push_back( (ELEMENT*)(toto.create(e->getType(), nodes) ));
+    }
+    to.push_back( (ELEMENT*)(toto.create(e->getTypeForMSH(), nodes) ));
   }
 }
 
@@ -693,7 +703,6 @@ int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
         //std::vector<Pair<GRegion*, GRegion*> >* coresp_r =
         matchRegions(geom, mesh, coresp_f,ok);
         if (ok) {
-          //mesh->writeMSH("out.msh",2.0,false,true);
 	  std::map<MVertex*,MVertex*> _mesh_to_geom;
 	  copy_vertices(geom, mesh, _mesh_to_geom);
 	  copy_elements(geom, mesh, _mesh_to_geom);

Geo/gmshFace.cpp

@@ -217,6 +217,7 @@ GPoint gmshFace::point(double par1, double par2) const
   }
 }
 
+#if not defined(QUASINEWTON)
 GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
 {
   if (s->Typ == MSH_SURF_PLAN && !s->geometry){
@@ -257,6 +258,7 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
     return GPoint(-1.e22, -1.e22, -1.e22, 0, u);
   return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
 }
+#endif
 
 SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
 {

Geo/gmshFace.h

@@ -21,9 +21,11 @@ class gmshFace : public GFace {
   virtual ~gmshFace(){}
   Range<double> parBounds(int i) const; 
   void setModelEdges(std::list<GEdge*> &);
-  virtual GPoint point(double par1, double par2) const; 
+  virtual GPoint point(double par1, double par2) const;
+#if not defined(QUASINEWTON)
   virtual GPoint closestPoint(const SPoint3 &queryPoint, 
                               const double initialGuess[2]) const; 
+#endif
   virtual bool containsPoint(const SPoint3 &pt) const;  
   virtual double getMetricEigenvalue(const SPoint2 &);  
   virtual SVector3 normal(const SPoint2 &param) const;