diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 7025cf611d8671a334a600641363bfebb505c550..a6a0498312de1560caf1d11b31b2758cd59aa685 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -109,7 +109,10 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
while(wire.size()){
ges = nextOne(prevOne, wire);
if(ges.getSign() == 0){ // oops
- Msg::Error("Something wrong in edge loop, no sign !");
+ Msg::Error("Something wrong in edge loop of size=%d, no sign !", wire.size());
+ for (std::list<GEdge* >::iterator it = wire.begin(); it != wire.end(); it++){
+ printf("GEdge=%d begin=%g end =%d \n", (*it)->tag(), (*it)->getBeginVertex()->tag(), (*it)->getEndVertex()->tag() );
+ }
break;
}
prevOne = ⩾
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e0b4a2bbb6b8085eb0b601c0940115b46d3a260a..218e20f1dde1b1846e7f76ed9523b061de0c001c 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -426,13 +426,13 @@ bool GFaceCompound::checkOrientation(int iter) const
int iterMax = 10;
if(!oriented && iter < iterMax){
- if (iter == 0) Msg::Warning("--- Parametrization is NOT 1 to 1 : applying cavity checks.");
+ if (iter == 0) Msg::Warning("--- Parametrization is flipping : applying cavity checks.");
Msg::Debug("--- Cavity Check - iter %d -",iter);
one2OneMap();
return checkOrientation(iter+1);
}
else if (oriented && iter < iterMax){
- Msg::Info("Parametrization is 1 to 1 :-)");
+ Msg::Info("Parametrization has no flips :-)");
//printStuff();
}
@@ -521,10 +521,12 @@ bool GFaceCompound::parametrize() const
bool withoutFolding = parametrize_conformal_spectral() ;
//bool withoutFolding = parametrize_conformal();
if ( withoutFolding == false ){
- //buildOct(); exit(1);
- Msg::Warning("$$$ Parametrization switched to harmonic map");
- parametrize(ITERU,HARMONIC);
- parametrize(ITERV,HARMONIC);
+
+ double alpha = 0.0;
+ Msg::Warning("$$$ Parametrization switched to combination map: A*conf+(1-A)*harm with A=%g", alpha);
+ parametrize(ITERU,HARMONIC, alpha);
+ parametrize(ITERV,HARMONIC, alpha);
+ //buildOct(); exit(1);
}
}
@@ -982,7 +984,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
return SPoint2(0, 0);
}
-void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
+void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom, double alpha) const
{
dofManager<double> myAssembler(_lsys);
@@ -1098,14 +1100,22 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
MVertex *v = *itv;
double value;
myAssembler.getDofValue(v, 0, 1, value);
- std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
+ std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
+
+ //combination convex and harmonic
+ double valNEW ;
+ if (alpha != 0.0)
+ valNEW = alpha*itf->second[step] + (1-alpha)*value ;
+ else
+ valNEW = value;
+
if(itf == coordinates.end()){
SPoint3 p(0, 0, 0);
- p[step] = value;
+ p[step] = valNEW;
coordinates[v] = p;
}
else{
- itf->second[step]= value;
+ itf->second[step]= valNEW;
}
}
@@ -1181,7 +1191,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//mettre max NC contraintes par bord
int NB = ordered.size();
- int NC = std::min(100,NB);
+ int NC = std::min(200,NB);
int jump = (int) NB/NC;
int nbLoop = (int) NB/jump ;
//printf("nb bound nodes=%d jump =%d \n", NB, jump);
@@ -1197,8 +1207,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//-------------------------------
//printf("Solve eigensystem \n");
eigenSolver eig(&myAssembler, "B" , "A", true);
- int nb = 2;
- bool converged = eig.solve(nb, "largest");
+ bool converged = eig.solve(1, "largest");
if(converged) {
double Linfty = 0.0;
@@ -1218,22 +1227,6 @@ bool GFaceCompound::parametrize_conformal_spectral() const
coordinates[v] = SPoint3(paramu,paramv,0.0);
k = k+2;
}
-
-
- //if folding take second smallest eigenvalue
- bool noFolding = checkFolding(ordered);
- if (!noFolding && nb > 1){
- coordinates.clear();
- int k = 0;
- std::vector<std::complex<double> > &ev = eig.getEigenVector(nb-1);
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- double paramu = ev[k].real();
- double paramv = ev[k+1].real();
- coordinates[v] = SPoint3(paramu,paramv,0.0);
- k = k+2;
- }
- }
lsysA->clear();
lsysB->clear();
@@ -1263,9 +1256,9 @@ bool GFaceCompound::parametrize_conformal() const
MVertex *v1 = ordered[0];
MVertex *v2 = ordered[(int)ceil((double)ordered.size()/2.)];
- myAssembler.fixVertex(v1, 0, 1, 0.);
+ myAssembler.fixVertex(v1, 0, 1, 1.);
myAssembler.fixVertex(v1, 0, 2, 0.);
- myAssembler.fixVertex(v2, 0, 1, 1.);
+ myAssembler.fixVertex(v2, 0, 1, -1.);
myAssembler.fixVertex(v2, 0, 2, 0.);
//printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 887545ec052ac544370294d9b15d58614fdee828..43009c6207d9b9c7390b3b15c69e454622a7a034 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -72,7 +72,7 @@ class GFaceCompound : public GFace {
mutable std::list<MTriangle*> fillTris;
void buildOct() const ;
void buildAllNodes() const;
- void parametrize(iterationStep, typeOfMapping) const;
+ void parametrize(iterationStep, typeOfMapping, double alpha=0.) const;
bool parametrize_conformal() const;
bool parametrize_conformal_spectral() const;
void compute_distance() const;
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 30269164ae813e9528da8a53d5d84f01e1d3ba86..2015481ca59bd12a5b3de3710fb7966eaf47a06d 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -88,7 +88,7 @@ void lloydAlgorithm::operator () (GFace *gf)
ENERGY += E;
double d = sqrt((p.x()-cgs(i,0))*(p.x()-cgs(i,0))+
(p.y()-cgs(i,1))*(p.y()-cgs(i,1)));
- criteria = d/A;
+ criteria += d/A;
}// if (v->onWhat() == gf)
else {
}
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index cf67e2addafc2f706a55ca8c0cdcdab321fdb265..8dc433a105d4d023557321ce99b8e20744bd4ba5 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -23,6 +23,21 @@
#include "ANN/ANN.h"
#endif
+StringXNumber DistanceOptions_Number[] = {
+ {GMSH_FULLRC, "Point", NULL, 0.},
+ {GMSH_FULLRC, "Line", NULL, 0.},
+ {GMSH_FULLRC, "Face", NULL, 0.},
+ {GMSH_FULLRC, "Computation", NULL, 0.},
+ //{GMSH_FULLRC, "Scale", NULL, 0.},
+ //{GMSH_FULLRC, "Min Scale", NULL, 1.e-3},
+ //{GMSH_FULLRC, "Max Scale", NULL, 1}
+};
+
+StringXString DistanceOptions_String[] = {
+ {GMSH_FULLRC, "Filename", NULL, "distance.pos"}
+};
+
+
extern "C"
{
GMSH_Plugin *GMSH_RegisterDistancePlugin()
@@ -33,14 +48,48 @@ extern "C"
std::string GMSH_DistancePlugin::getHelp() const
{
- return "Plugin(Distance) computes distances to boundaries in "
+ return "Plugin(Distance) computes distances to elementary entities in "
"a mesh.\n\n"
- "Plugin(Distance) creates a bunch of files.";
+
+ "Define the elementary entities to which the distance is computed \n\n";
+
+ "Computation=0 computes the geometrical distance (Warning: this is an euclidian distance and not the geodesic distance), and Computation=1 solves a PDE on the mesh\n\n";
+
+ "Plugin(Distance) creates a new distance view.";
+}
+
+int GMSH_DistancePlugin::getNbOptions() const
+{
+ return sizeof(DistanceOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_DistancePlugin::getOption(int iopt)
+{
+ return &DistanceOptions_Number[iopt];
+}
+
+int GMSH_DistancePlugin::getNbOptionsStr() const
+{
+ return sizeof(DistanceOptions_String) / sizeof(StringXString);
+}
+
+StringXString *GMSH_DistancePlugin::getOptionStr(int iopt)
+{
+ return &DistanceOptions_String[iopt];
}
PView *GMSH_DistancePlugin::execute(PView *v)
{
-#if defined(HAVE_SOLVER)
+
+ std::string fileName = DistanceOptions_String[0].def;
+ int id_pt = (int) DistanceOptions_Number[0].def;
+ int id_line = (int) DistanceOptions_Number[1].def;
+ int id_face = (int) DistanceOptions_Number[2].def;
+ int type = (int) DistanceOptions_Number[3].def;
+
+ PView *view = new PView();
+ PViewDataList *data = getDataList(view);
+
std::vector<GEntity*> entities;
GModel::current()->getEntities(entities);
@@ -49,8 +98,8 @@ PView *GMSH_DistancePlugin::execute(PView *v)
numnodes += entities[i]->mesh_vertices.size();
int totNumNodes = numnodes + entities[entities.size() - 1]->mesh_vertices.size();
- std::map<MVertex*,double* > distance_map;
- std::map<MVertex*,double > distance_map2;
+ std::map<MVertex*,double* > distance_map_GEOM;
+ std::map<MVertex*,double > distance_map_PDE;
std::vector<SPoint3> pts;
std::vector<double> distances;
pts.clear();
@@ -67,279 +116,227 @@ PView *GMSH_DistancePlugin::execute(PView *v)
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])));
- distance_map2.insert(std::make_pair(v, 0.0));
+ distance_map_GEOM.insert(std::make_pair(v, &(distances[k])));
+ distance_map_PDE.insert(std::make_pair(v, 0.0));
k++;
}
}
-
+
// Compute geometrical distance to mesh boundaries
+ //------------------------------------------------------
+ if (type == 0){
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() == (maxDim - 1)) {
+ for(unsigned int i = 0; i < entities.size(); i++){
GEntity* g2 = entities[i];
- for(unsigned int k = 0; k < g2->getNumMeshElements(); k++){
- std::vector<double> iDistances;
- std::vector<SPoint3> iClosePts;
- MElement *e = g2->getMeshElement(k);
- 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());
- if(e->getNumVertices() == 2){
- signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
- }
- else if(e->getNumVertices() == 3){
- MVertex *v3 = e->getVertex(2);
- SPoint3 p3 (v3->x(),v3->y(),v3->z());
- signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
- }
- for (unsigned int kk = 0; kk< pts.size(); kk++) {
- if (std::abs(iDistances[kk]) < distances[kk])
- distances[kk] = std::abs(iDistances[kk]);
- }
+ int tag = g2->tag();
+ int gDim = g2->dim();
+ bool computeForEntity = false;
+ if (id_pt==0 && id_line==0 && id_face==0 && gDim==maxDim-1 )
+ computeForEntity = true;
+ else if ( (tag==id_pt && gDim==0)|| (tag==id_line && gDim==1) || (tag==id_face && gDim==2) )
+ computeForEntity = true;
+ if (computeForEntity){
+ for(unsigned int k = 0; k < g2->getNumMeshElements(); k++){
+ std::vector<double> iDistances;
+ std::vector<SPoint3> iClosePts;
+ MElement *e = g2->getMeshElement(k);
+ 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());
+ if(e->getNumVertices() == 2){
+ signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+ }
+ else if(e->getNumVertices() == 3){
+ MVertex *v3 = e->getVertex(2);
+ SPoint3 p3 (v3->x(),v3->y(),v3->z());
+ signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
+ }
+ for (unsigned int kk = 0; kk< pts.size(); kk++) {
+ if (std::abs(iDistances[kk]) < distances[kk])
+ distances[kk] = std::abs(iDistances[kk]);
+ }
+ }
}
}
- }
+ // std::map<int, std::vector<GEntity*> > groups[4];
+ // getPhysicalGroups(groups);
+ // std::map<int, std::vector<GEntity*> >::const_iterator it = groups[1].find(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::map<int, std::vector<GEntity*> > groups[4];
-// getPhysicalGroups(groups);
-// std::map<int, std::vector<GEntity*> >::const_iterator it = groups[1].find(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]);
-// }
-// }
-// }
+ data->setName("distance GEOM");
+ Msg::Info("Writing %g", fileName.c_str());
+ FILE * f3 = fopen( fileName.c_str(),"w");
+ fprintf(f3,"View \"distance GEOM\"{\n");
+ for(unsigned int ii = 0; ii < entities.size(); ii++){
+ if(entities[ii]->dim() == maxDim) {
+ for(unsigned int i = 0; i < entities[ii]->getNumMeshElements(); i++){
+ MElement *e = entities[ii]->getMeshElement(i);
+ int numNodes = e->getNumVertices();
+ std::vector<double> x(numNodes), y(numNodes), z(numNodes);
+ std::vector<double> *out = data->incrementList(1, e->getType());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->x());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->y());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->z());
+ if (maxDim == 3) fprintf(f3,"SS(");
+ else if (maxDim == 2) fprintf(f3,"ST(");
+ std::vector<double> dist;
+ for(int j = 0; j < numNodes; 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_GEOM.find(v);
+ dist.push_back(*(it->second));
+ }
+ fprintf(f3,"){");
+ for (unsigned int i = 0; i < dist.size(); i++){
+ out->push_back(dist[i]);
+ if (i) fprintf(f3,",%g", dist[i]);
+ else fprintf(f3,"%g", dist[i]);
+ }
+ fprintf(f3,"};\n");
+ }
+ }
+ }
+ fprintf(f3,"};\n");
+ fclose(f3);
+ }
+
// Compute PDE for distance function
+ //-----------------------------------
+ else if (type == 1){
+
+#if defined(HAVE_SOLVER)
#ifdef HAVE_TAUCS
- linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
- linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
- lsys->setNoisy(1);
- lsys->setGmres(1);
- lsys->setPrec(5.e-8);
+ linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
+ lsys->setNoisy(1);
+ lsys->setGmres(1);
+ lsys->setPrec(5.e-8);
#endif
- dofManager<double> myAssembler(lsys);
+ dofManager<double> myAssembler(lsys);
- SBoundingBox3d bbox;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < maxDim) {
- GEntity *ge = entities[i];
- for(unsigned int j = 0; j < ge->mesh_vertices.size(); j++){
- MVertex *v = ge->mesh_vertices[j];
- myAssembler.fixVertex(v, 0, 1, 0.0);
- bbox += SPoint3(v->x(),v->y(),v->z());
+ SBoundingBox3d bbox;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity* ge = entities[i];
+ int tag = ge->tag();
+ int gDim = ge->dim();
+ bool fixForEntity = false;
+ if(id_pt==0 && id_line==0 && id_face==0 && gDim < maxDim)
+ fixForEntity = true;
+ else if ( (tag==id_pt && gDim==0)|| (tag==id_line && gDim==1) || (tag==id_face && gDim==2) )
+ fixForEntity = true;
+ if (fixForEntity){
+ for(unsigned int j = 0; j < ge->mesh_vertices.size(); j++){
+ MVertex *v = ge->mesh_vertices[j];
+ myAssembler.fixVertex(v, 0, 1, 0.0);
+ bbox += SPoint3(v->x(),v->y(),v->z());
+ }
}
}
- }
- for(unsigned int ii = 0; ii < entities.size(); ii++){
- if(entities[ii]->dim() == maxDim) {
- GEntity *ge = entities[ii];
- for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
- MElement *t = ge->getMeshElement(i);
- for(int k = 0; k < t->getNumVertices(); k++){
- myAssembler.numberVertex(t->getVertex(k), 0, 1);
- }
- }
- }
- }
+ for(unsigned int ii = 0; ii < entities.size(); ii++){
+ if(entities[ii]->dim() == maxDim) {
+ GEntity *ge = entities[ii];
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
+ MElement *t = ge->getMeshElement(i);
+ for(int k = 0; k < t->getNumVertices(); k++){
+ myAssembler.numberVertex(t->getVertex(k), 0, 1);
+ }
+ }
+ }
+ }
- double L = norm(SVector3(bbox.max(), bbox.min()));
- double mu = L/28;
+ double L = norm(SVector3(bbox.max(), bbox.min()));
+ double mu = L/10; //28;
- simpleFunction<double> DIFF(mu * mu), MONE(1.0);
- distanceTerm distance(GModel::current(), 1, &DIFF, &MONE);
+ simpleFunction<double> DIFF(mu * mu), MONE(1.0);
+ distanceTerm distance(GModel::current(), 1, &DIFF, &MONE);
- for(unsigned int ii = 0; ii < entities.size(); ii++){
- if(entities[ii]->dim() == maxDim) {
- GEntity *ge = entities[ii];
- for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
- SElement se(ge->getMeshElement(i));
- distance.addToMatrix(myAssembler, &se);
+ for(unsigned int ii = 0; ii < entities.size(); ii++){
+ if(entities[ii]->dim() == maxDim) {
+ GEntity *ge = entities[ii];
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i){
+ SElement se(ge->getMeshElement(i));
+ distance.addToMatrix(myAssembler, &se);
+ }
+ groupOfElements g((GFace*)ge);
+ distance.addToRightHandSide(myAssembler, g);
}
- groupOfElements g((GFace*)ge);
- distance.addToRightHandSide(myAssembler, g);
}
- }
- Msg::Info("Distance Computation: Assembly done");
- lsys->systemSolve();
- Msg::Info("Distance Computation: System solved");
+ Msg::Info("Distance Computation: Assembly done");
+ lsys->systemSolve();
+ Msg::Info("Distance Computation: System solved");
- for(std::map<MVertex*,double >::iterator itv = distance_map2.begin();
- itv !=distance_map2.end() ; ++itv){
- MVertex *v = itv->first;
- double value;
- myAssembler.getDofValue(v, 0, 1, value);
- value = std::min(0.9999, value);
- double dist = -mu * log(1. - value);
- itv->second = dist;
- }
+ for(std::map<MVertex*,double >::iterator itv = distance_map_PDE.begin();
+ itv !=distance_map_PDE.end() ; ++itv){
+ MVertex *v = itv->first;
+ double value;
+ myAssembler.getDofValue(v, 0, 1, value);
+ value = std::min(0.9999, value);
+ double dist = -mu * log(1. - value);
+ itv->second = dist;
+ }
+ lsys->clear();
- lsys->clear();
+ data->setName("distance GEOM");
- //Writing pos files
-
- 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");
- for(unsigned int ii = 0; ii < entities.size(); ii++){
- if(entities[ii]->dim() == maxDim) {
- for(unsigned int i = 0; i < entities[ii]->getNumMeshElements(); i++){
- MElement *e = entities[ii]->getMeshElement(i);
- if (maxDim == 3) fprintf(f3,"SS(");
- else if (maxDim == 2) fprintf(f3,"ST(");
- 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);
- dist.push_back(*(it->second));
- }
- fprintf(f3,"){");
- for (unsigned 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);
-
- Msg::Info("Writing distance-bgm2.pos");
- FILE * f4 = fopen("distance-bgm2.pos","w");
- fprintf(f4,"View \"distance\"{\n");
- for(unsigned int ii = 0; ii < entities.size(); ii++){
- if(entities[ii]->dim() == maxDim) {
- for(unsigned int i = 0; i < entities[ii]->getNumMeshElements(); i++){
- MElement *e = entities[ii]->getMeshElement(i);
- if (maxDim == 3) fprintf(f4,"SS(");
- else if (maxDim == 2) fprintf(f4,"ST(");
- std::vector<double> dist;
- for(int j = 0; j < e->getNumVertices(); j++) {
- MVertex *v = e->getVertex(j);
- if(j) fprintf(f4,",%g,%g,%g",v->x(),v->y(), v->z());
- else fprintf(f4,"%g,%g,%g", v->x(),v->y(), v->z());
- std::map<MVertex*, double>::iterator it = distance_map2.find(v);
- dist.push_back(it->second);
- }
- fprintf(f4,"){");
- for (unsigned int i = 0; i < dist.size(); i++){
- if (i) fprintf(f4,",%g", dist[i]);
- else fprintf(f4,"%g", dist[i]);
- }
- fprintf(f4,"};\n");
+ Msg::Info("Writing %d", fileName.c_str());
+ FILE * f4 = fopen(fileName.c_str(),"w");
+ fprintf(f4,"View \"distance PDE\"{\n");
+ for(unsigned int ii = 0; ii < entities.size(); ii++){
+ if(entities[ii]->dim() == maxDim) {
+ for(unsigned int i = 0; i < entities[ii]->getNumMeshElements(); i++){
+ MElement *e = entities[ii]->getMeshElement(i);
+
+ int numNodes = e->getNumVertices();
+ std::vector<double> x(numNodes), y(numNodes), z(numNodes);
+ std::vector<double> *out = data->incrementList(1, e->getType());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->x());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->y());
+ for(int nod = 0; nod < numNodes; nod++) out->push_back((e->getVertex(nod))->z());
+
+ if (maxDim == 3) fprintf(f4,"SS(");
+ else if (maxDim == 2) fprintf(f4,"ST(");
+ std::vector<double> dist;
+ for(int j = 0; j < numNodes; j++) {
+ MVertex *v = e->getVertex(j);
+ if(j) fprintf(f4,",%g,%g,%g",v->x(),v->y(), v->z());
+ else fprintf(f4,"%g,%g,%g", v->x(),v->y(), v->z());
+ std::map<MVertex*, double>::iterator it = distance_map_PDE.find(v);
+ dist.push_back(it->second);
+ }
+ fprintf(f4,"){");
+ for (unsigned int i = 0; i < dist.size(); i++){
+ out->push_back(dist[i]);
+ if (i) fprintf(f4,",%g", dist[i]);
+ else fprintf(f4,"%g", dist[i]);
+ }
+ fprintf(f4,"};\n");
+ }
}
}
- }
- fprintf(f4,"};\n");
- fclose(f4);
-
-// #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(f2,"View \"distance\"{\n");
-// GEntity* ge = entities[entities.size()-1];
-// for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
-// MElement *e = ge->getMeshElement(i);
-// fprintf(f2,"ST(");//TO DO CHANGE this
-// std::vector<double> dist;
-// for(int j = 0; j < e->getNumVertices(); j++) {
-// MVertex *v = e->getVertex(j);
-// if(j) fprintf(f2,",%g,%g,%g",v->x(),v->y(), v->z());
-// else fprintf(f2,"%g,%g,%g", v->x(),v->y(), v->z());
-// std::map<MVertex*, double>::iterator it = distance.find(v);
-// dist.push_back(it->second);
-// }
-// fprintf(f2,"){");
-// for (int i=0; i<dist.size(); i++){
-// if (i) fprintf(f2,",%g", dist[i]);
-// else fprintf(f2,"%g", dist[i]);
-// }
-// fprintf(f2,"};\n");
-// }
-// fprintf(f2,"};\n");
-// fclose(f2);
+ fprintf(f4,"};\n");
+ fclose(f4);
#endif
- return 0;
+ }
+
+ //FILE *fp = fopen(fileName.c_str(), "rb");
+ //data->readPOS(fp, 1.0, false);
+
+ data->Time.push_back(0);
+ data->setFileName(fileName.c_str());
+ data->finalize();
+
+ return view;
}
diff --git a/Plugin/Distance.h b/Plugin/Distance.h
index 35cd24d1d9b474d9227bb02f0cda589fc69c44c7..e8f05d50e347e096bf7d6db922b0ec948aeff962 100644
--- a/Plugin/Distance.h
+++ b/Plugin/Distance.h
@@ -27,6 +27,10 @@ class GMSH_DistancePlugin : public GMSH_PostPlugin
}
std::string getHelp() const;
std::string getAuthor() const { return "E. Marchandise"; }
+ int getNbOptions() const;
+ StringXNumber *getOption(int iopt);
+ int getNbOptionsStr() const;
+ StringXString *getOptionStr(int iopt);
PView *execute(PView *);
};
diff --git a/Post/PViewIO.cpp b/Post/PViewIO.cpp
index 58a1b2d8cbbb199db4cb8115e441e48cc69e155d..b1dc23bb90c19069f0b450a5236b6a77cb556a22 100644
--- a/Post/PViewIO.cpp
+++ b/Post/PViewIO.cpp
@@ -51,7 +51,6 @@ bool PView::readPOS(std::string fileName, int fileIndex)
}
else if(!strncmp(&str[1], "View", 4)){
-
index++;
if(fileIndex < 0 || fileIndex == index){
PViewDataList *d = new PViewDataList();