diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 7fb5b4b2809e0e01789505ff2e6e51e3cfabe73c..0e623b4bd5bd6a169c27d9abe1d36e0b25db88d2 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -833,64 +833,166 @@ int GModel::writeDistanceMSH(const std::string &name, double version, bool binar
int numnodes = 0;
for(unsigned int i = 0; i < entities.size()-1; i++)
numnodes += entities[i]->mesh_vertices.size();
-
int totNumNodes = numnodes + entities[entities.size()-1]->mesh_vertices.size();
- std::map<MVertex*,double > distance;
-
- // EMI: TODO COMPUTE this with EDP distanceTERM instead of ANN write
- // solution for elements
// CG: Emi -- pourquoi ne pas en faire un plugin, plutot ?
-#if defined(HAVE_ANN)
- Msg::Info("Computing Distance function to the boundaries with ANN");
-
- // add all points from boundaries into kdtree
- ANNpointArray kdnodes = annAllocPts(numnodes, 4);
- int k = 0;
- for(unsigned int i = 0; i < entities.size()-1; i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
- MVertex *v = entities[i]->mesh_vertices[j];
- distance.insert(std::make_pair(v, 0.0));
- kdnodes[k][0] = v->x();
- kdnodes[k][1] = v->y();
- kdnodes[k][2] = v->z();
- k++;
- }
- ANNkd_tree *kdtree = new ANNkd_tree(kdnodes, numnodes, 3);
-
- // compute the distances with the kdtree
- ANNidxArray index = new ANNidx[1];
- ANNdistArray dist = new ANNdist[1];
- GEntity* g = entities[entities.size()-1];
- for(unsigned int j = 0; j < g->mesh_vertices.size(); j++){
- MVertex *v = g->mesh_vertices[j];
- double xyz[3] = {v->x(), v->y(), v->z()};
- kdtree->annkSearch(xyz, 1, index, dist);
- double d = sqrt(dist[0]);
- distance.insert(std::make_pair(v, d));
- }
- delete [] index;
- delete [] dist;
- delete kdtree;
- annDeallocPts(kdnodes);
-#endif
+ //geometrical distance
+
+ std::map<MVertex*,double* > distance_map;
+ std::vector<SPoint3> pts;
+ std::vector<double> distances;
+ pts.clear();
+ distances.clear();
+ distances.reserve(totNumNodes);
+ pts.reserve(totNumNodes);
+ for (int i = 0; i< totNumNodes; i++) distances.push_back(1.e22);
- //write distance in msh file
- Msg::Info("Writing distance.pos");
- FILE * f = fopen("distance.pos","w");
- fprintf(f,"View \"distance\"{\n");
- for(std::map<MVertex*, double>::iterator it = distance.begin(); it != distance.end(); it++) {
- fprintf(f,"SP(%g,%g,%g){%g};\n",
- it->first->x(), it->first->y(), it->first->z(),
- it->second);
+ int k=0;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity* ge = entities[i];
+ for(unsigned int j = 0; j < ge->mesh_vertices.size(); j++){
+ MVertex *v = ge->mesh_vertices[j];
+ pts.push_back(SPoint3(v->x(), v->y(),v->z()));
+ distance_map.insert(std::make_pair(v, &(distances[k])));
+ k++;
+ }
}
- fprintf(f,"};\n");
- fclose(f);
+
+// GEntity* g2 = entities[entities.size()-2];//faces
+// for(unsigned int i = 0; i < g2->getNumMeshElements(); i++){
+// std::vector<double> iDistances;
+// std::vector<SPoint3> iClosePts;
+// MElement *e = g2->getMeshElement(i);
+// MVertex *v1 = e->getVertex(0);
+// MVertex *v2 = e->getVertex(1);
+// MVertex *v3 = e->getVertex(2);
+// SPoint3 p1 (v1->x(),v1->y(),v1->z());
+// SPoint3 p2 (v2->x(),v2->y(),v2->z());
+// SPoint3 p3 (v3->x(),v3->y(),v3->z());
+// //signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1,p2,p3);
+// signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+// signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p3);
+// signedDistancesPointsLine(iDistances, iClosePts, pts, p2,p3);
+// for (int k = 0; k< pts.size(); k++) {
+// if (std::abs(iDistances[k]) < distances[k] )
+// distances[k] = std::abs(iDistances[k]);
+// }
+// }
+
+ std::map<int, std::vector<GEntity*> > groups[4];
+ getPhysicalGroups(groups);
+ std::map<int, std::vector<GEntity*> >::const_iterator it = groups[1].find(100);//Physical Line(100)={...}
+ if(it == groups[1].end()) return 0;
+ const std::vector<GEntity *> &physEntities = it->second;
+ for(unsigned int i = 0; i < physEntities.size(); i++){
+ GEntity* g2 = physEntities[i];
+ for(unsigned int i = 0; i < g2->getNumMeshElements(); i++){
+ std::vector<double> iDistances;
+ std::vector<SPoint3> iClosePts;
+ MElement *e = g2->getMeshElement(i);
+ MVertex *v1 = e->getVertex(0);
+ MVertex *v2 = e->getVertex(1);
+ SPoint3 p1 (v1->x(),v1->y(),v1->z());
+ SPoint3 p2 (v2->x(),v2->y(),v2->z());
+ signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+ for (int k = 0; k< pts.size(); k++) {
+ if (std::abs(iDistances[k]) < distances[k] )
+ distances[k] = std::abs(iDistances[k]);
+ }
+ }
+ }
+
+ Msg::Info("Writing distance-GEOM.pos");
+ FILE * f2 = fopen("distance-GEOM.pos","w");
+ fprintf(f2,"View \"distance GEOM\"{\n");
+ int j = 0;
+ for(std::vector<double>::iterator it = distances.begin(); it != distances.end(); it++) {
+ fprintf(f2,"SP(%g,%g,%g){%g};\n",
+ pts[j].x(), pts[j].y(), pts[j].z(),
+ *it);
+ j++;
+ }
+ fprintf(f2,"};\n");
+ fclose(f2);
+
+ Msg::Info("Writing distance-bgm.pos");
+ FILE * f3 = fopen("distance-bgm.pos","w");
+ fprintf(f3,"View \"distance\"{\n");
+ GEntity* ge = entities[entities.size()-1];
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement(i);
+ fprintf(f3,"SS(");
+ std::vector<double> dist;
+ for(int j = 0; j < e->getNumVertices(); j++) {
+ MVertex *v = e->getVertex(j);
+ if(j) fprintf(f3,",%g,%g,%g",v->x(),v->y(), v->z());
+ else fprintf(f3,"%g,%g,%g", v->x(),v->y(), v->z());
+ std::map<MVertex*, double*>::iterator it = distance_map.find(v);
+ printf("dist =%g \n",*(it->second) );
+ dist.push_back(*(it->second));
+ }
+ fprintf(f3,"){");
+ for (int i=0; i<dist.size(); i++){
+ if (i) fprintf(f3,",%g", dist[i]);
+ else fprintf(f3,"%g", dist[i]);
+ }
+ fprintf(f3,"};\n");
+ }
+ fprintf(f3,"};\n");
+ fclose(f3);
+
+// #if defined(HAVE_ANN)
+// Msg::Info("Computing Distance function to the boundaries with ANN");
+// std::map<MVertex*,double > distance;
+
+// // add all points from boundaries into kdtree
+// ANNpointArray kdnodes = annAllocPts(numnodes, 4);
+// k = 0;
+// for(unsigned int i = 0; i < entities.size()-1; i++)
+// for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+// MVertex *v = entities[i]->mesh_vertices[j];
+// distance.insert(std::make_pair(v, 0.0));
+// kdnodes[k][0] = v->x();
+// kdnodes[k][1] = v->y();
+// kdnodes[k][2] = v->z();
+// k++;
+// }
+// ANNkd_tree *kdtree = new ANNkd_tree(kdnodes, numnodes, 3);
+
+// // compute the distances with the kdtree
+// ANNidxArray index = new ANNidx[1];
+// ANNdistArray dist = new ANNdist[1];
+// GEntity* ga = entities[entities.size()-1];
+// for(unsigned int j = 0; j < ga->mesh_vertices.size(); j++){
+// MVertex *v = ga->mesh_vertices[j];
+// double xyz[3] = {v->x(), v->y(), v->z()};
+// kdtree->annkSearch(xyz, 1, index, dist);
+// double d = sqrt(dist[0]);
+// distance.insert(std::make_pair(v, d));
+// }
+// delete [] index;
+// delete [] dist;
+// delete kdtree;
+// annDeallocPts(kdnodes);
+
+// //write distance in msh file
+// Msg::Info("Writing distance-ANN.pos");
+// FILE * f = fopen("distance-ANN.pos","w");
+// fprintf(f,"View \"distance ANN\"{\n");
+// for(std::map<MVertex*, double>::iterator it = distance.begin(); it != distance.end(); it++) {
+// fprintf(f,"SP(%g,%g,%g){%g};\n",
+// it->first->x(), it->first->y(), it->first->z(),
+// it->second);
+// }
+// fprintf(f,"};\n");
+// fclose(f);
+
+// #endif
// Msg::Info("Writing distance-bgm.pos");
// FILE * f2 = fopen("distance-bgm.pos","w");
-// fprintf(f,"View \"distance\"{\n");
+// fprintf(f2,"View \"distance\"{\n");
// GEntity* ge = entities[entities.size()-1];
// for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
// MElement *e = ge->getMeshElement(i);
@@ -910,10 +1012,12 @@ int GModel::writeDistanceMSH(const std::string &name, double version, bool binar
// }
// fprintf(f2,"};\n");
// }
-// fprintf(f,"};\n");
+// fprintf(f2,"};\n");
// fclose(f2);
-
+
+
return 1;
+
}
int GModel::writePOS(const std::string &name, bool printElementary,
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 3ead5be996847ebd2a1e863192da613571cffea6..f60a0e713b93b5dfe0ab3828da036c5d6d64fbae 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -723,6 +723,7 @@ distance to segment
*/
void signedDistancesPointsTriangle(std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
const std::vector<SPoint3> &pts,
const SPoint3 &p1,
const SPoint3 &p2,
@@ -742,6 +743,8 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
distances.clear();
distances.resize(pts.size());
+ closePts.clear();
+ closePts.resize(pts.size());
for (unsigned int i = 0; i < pts.size(); i++)
distances[i] = 1.e22;
@@ -763,6 +766,7 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
if (d == 0) sign = 1.e10;
if (u >= 0 && v >= 0 && 1.-u-v >= 0.0){
distances[i] = d;
+ closePts[i] = SPoint3(0.,0.,0.);//TO DO
}
else {
const double t12 = dot(pp1, t1) / n2t1;
@@ -771,21 +775,80 @@ void signedDistancesPointsTriangle(std::vector<double>&distances,
const double t23 = dot(pp2, t3) / n2t3;
d = 1.e10;
bool found = false;
+ SPoint3 closePt;
if (t12 >= 0 && t12 <= 1.){
- d = sign * std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
+ d = sign * std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
+ closePt = p1 + (p2 - p1) * t12;
found = true;
}
if (t13 >= 0 && t13 <= 1.){
+ if (p.distance(p1 + (p3 - p1) * t13) < fabs(d)) closePt = p1 + (p3 - p1) * t13;
d = sign * std::min(fabs(d), p.distance(p1 + (p3 - p1) * t13));
found = true;
}
if (t23 >= 0 && t23 <= 1.){
+ if (p.distance(p2 + (p3 - p2) * t23) < fabs(d)) closePt = p2 + (p3 - p2) * t23;
d = sign * std::min(fabs(d), p.distance(p2 + (p3 - p2) * t23));
found = true;
}
- d = sign * std::min(fabs(d), std::min(std::min(p.distance(p1), p.distance(p2)),
- p.distance(p3)));
+ if (p.distance(p1) < fabs(d)){
+ closePt = p1;
+ d = sign * std::min(fabs(d), p.distance(p1));
+ }
+ if (p.distance(p2) < fabs(d)){
+ closePt = p2;
+ d = sign * std::min(fabs(d), p.distance(p2));
+ }
+ if (p.distance(p3) < fabs(d)){
+ closePt = p3;
+ d = sign * std::min(fabs(d), p.distance(p3));
+ }
+ //d = sign * std::min(fabs(d), std::min(std::min(p.distance(p1), p.distance(p2)),p.distance(p3)));
distances[i] = d;
+ closePts[i] = closePt;
}
}
}
+
+void signedDistancesPointsLine (std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2){
+
+ SVector3 t1 = p2 - p1;
+ const double n2t1 = dot(t1, t1);
+
+ distances.clear();
+ distances.resize(pts.size());
+ closePts.clear();
+ closePts.resize(pts.size());
+
+ double d;
+ for (unsigned int i = 0; i < pts.size();i++){
+ const SPoint3 &p = pts[i];
+ SVector3 pp1 = p - p1;
+ const double t12 = dot(pp1, t1) / n2t1;
+ d = 1.e10;
+ bool found = false;
+ SPoint3 closePt;
+ if (t12 >= 0 && t12 <= 1.){
+ d = std::min(fabs(d), p.distance(p1 + (p2 - p1) * t12));
+ closePt = p1 + (p2 - p1) * t12;
+ found = true;
+ }
+
+ if (p.distance(p1) < fabs(d)){
+ closePt = p1;
+ d = std::min(fabs(d), p.distance(p1));
+ }
+ if (p.distance(p2) < fabs(d)){
+ closePt = p2;
+ d = std::min(fabs(d), p.distance(p2));
+ }
+
+ distances[i] = d;
+ closePts[i] = closePt;
+ }
+
+}
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 9678825457fccf2761f8591f74e6865389ffd83e..726f9c6ed526977ed39e112d765776930482d89b 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -78,8 +78,14 @@ bool newton_fd(void (*func)(fullVector<double> &, fullVector<double> &, void *),
double minimize_grad_fd(double (*func)(fullVector<double> &, void *),
fullVector<double> &x, void *data);
void signedDistancesPointsTriangle (std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
const std::vector<SPoint3> &pts,
const SPoint3 &p1,
const SPoint3 &p2,
const SPoint3 &p3);
+void signedDistancesPointsLine (std::vector<double>&distances,
+ std::vector<SPoint3>&closePts,
+ const std::vector<SPoint3> &pts,
+ const SPoint3 &p1,
+ const SPoint3 &p2);
#endif
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 58e7d028e54e138229c567f795040eb21c5b1b96..6d0f639489ee27a43f010bb4a4470ba2bb4a998d 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -215,8 +215,6 @@ static void recur_compute_centers_ (double R, double a1, double a2,
SPoint2 PL (R*cos(a1),R*sin(a1));
SPoint2 PR (R*cos(a2),R*sin(a2));
- printf("R=%g a1=%g PL=%g %g \n", R, a1*180/M_PI, PL.x(), PL.y());
- printf("R=%g a2=%g PR=%g %g \n", R, a2*180/M_PI, PR.x(), PR.y());
centers.push_back(std::make_pair(PL,zero));
centers.push_back(std::make_pair(PR,zero));
for (int i=0;i<root->children.size();i++){