diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 387ae02d0bdb2257408b0e62db8c90b6d1e97e8e..887740cbd529ccf5b1955137faabcc648fc9da72 100644
--- a/Common/OpenFile.cpp
+++ b/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);
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 9eda7e6c0a91eb57b3a92207f8b6ff256c6453d6..2aac4d62041ee378e27f35eec6335a9b1a5efeff 100644
--- a/Geo/GFace.cpp
+++ b/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
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 5f4b827b5d8464e8bd8b0bc39476a31c290b463a..93dca7bc3396ce42dbe74d78e0fa48fe21c991e2 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/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);
}
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index b56dfcb6b7d20c9375dd2df7511ef1719e1ff36b..07258e90aa688808a256fda75987ec8bfd80a2cb 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/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);
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index caca1f867ef1383a34031c1e7b87e5862229e871..9943da7498c3456c5febc45e059d22d381083da5 100644
--- a/Geo/gmshFace.cpp
+++ b/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
{
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index a0d1666c3c6e134e75708540f47bd36b2fd1ad98..67d0a3ce357f7d79b1d9cf7dbfa76a0e49882c6f 100644
--- a/Geo/gmshFace.h
+++ b/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 ¶m) const;