From 17451aa23eda6682b226bd1e746d098b4a04db82 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 20 Mar 2010 10:11:10 +0000
Subject: [PATCH] move distance computation to Plugin(Distance)
---
Common/CommandLine.cpp | 5 -
Common/Context.h | 2 +-
Common/CreateFile.cpp | 6 -
Common/DefaultOptions.h | 2 -
Common/Options.cpp | 6 -
Common/Options.h | 1 -
Geo/GModelIO_Mesh.cpp | 329 -------------------------------------
Plugin/CMakeLists.txt | 1 +
Plugin/Distance.cpp | 343 +++++++++++++++++++++++++++++++++++++++
Plugin/Distance.h | 33 ++++
Plugin/Plugin.h | 6 +-
Plugin/PluginManager.cpp | 3 +
12 files changed, 385 insertions(+), 352 deletions(-)
create mode 100644 Plugin/Distance.cpp
create mode 100644 Plugin/Distance.h
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 8b96247668..75e58c9525 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -66,7 +66,6 @@ void PrintUsage(const char *name)
Msg::Direct(" bdf, p3d, cgns, med, fea)");
Msg::Direct(" -bin Use binary format when available");
Msg::Direct(" -parametric Save vertices with their parametric coordinates");
- Msg::Direct(" -distance Save distance of vertices to the boundaries");
Msg::Direct(" -numsubedges Set the number of subdivisions when displaying high order elements");
Msg::Direct(" -algo string Select mesh algorithm (meshadapt, del2d, front2d, del3d, front3d)");
Msg::Direct(" -smooth int Set number of mesh smoothing steps");
@@ -496,10 +495,6 @@ void GetOptions(int argc, char *argv[])
i++;
CTX::instance()->mesh.saveParametric = 1;
}
- else if(!strcmp(argv[i] + 1, "distance")) {
- i++;
- CTX::instance()->mesh.saveDistance = 1;
- }
else if(!strcmp(argv[i] + 1, "algo")) {
i++;
if(argv[i]) {
diff --git a/Common/Context.h b/Common/Context.h
index 9d22f18bd4..0658b032cb 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -30,7 +30,7 @@ struct contextMeshOptions {
int order, secondOrderLinear, secondOrderIncomplete;
int secondOrderExperimental, meshOnlyVisible;
int smoothInternalEdges, minCircPoints, minCurvPoints;
- int saveAll, saveGroupsOfNodes, binary, bdfFieldFormat, saveParametric, saveDistance;
+ int saveAll, saveGroupsOfNodes, binary, bdfFieldFormat, saveParametric;
int smoothNormals, reverseAllNormals, zoneDefinition, clip;
int saveElementTagType;
};
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index c43b878361..6160fbc619 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -183,12 +183,6 @@ void CreateOutputFile(std::string fileName, int format)
CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
}
- if (CTX::instance()->mesh.saveDistance){
- GModel::current()->writeDistanceMSH
- (fileName, CTX::instance()->mesh.mshFileVersion,
- CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
- CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
- }
break;
case FORMAT_STL:
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 41d2d2c6ac..fe3e73d1b0 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1208,8 +1208,6 @@ StringXNumber MeshOptions_Number[] = {
"physical or partition ids (1=elementary, 2=physical, 3=partition)" },
{ F|O, "SaveParametric" , opt_mesh_save_parametric , 0. ,
"Save parametric coordinates of nodes" },
- { F|O, "SaveDistance" , opt_mesh_save_distance , 0. ,
- "Save distance of nodes to boundaries" },
{ F|O, "SaveGroupsOfNodes" , opt_mesh_save_groups_of_nodes , 0. ,
"Save groups of nodes for each physical line and surface (UNV mesh "
"format only)" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 2cf5913fe1..f1077ee938 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5695,12 +5695,6 @@ double opt_mesh_save_parametric(OPT_ARGS_NUM)
CTX::instance()->mesh.saveParametric = val ? 1 : 0;
return CTX::instance()->mesh.saveParametric;
}
-double opt_mesh_save_distance(OPT_ARGS_NUM)
-{
- if(action & GMSH_SET)
- CTX::instance()->mesh.saveDistance = val ? 1 : 0;
- return CTX::instance()->mesh.saveDistance;
-}
double opt_mesh_save_groups_of_nodes(OPT_ARGS_NUM)
{
diff --git a/Common/Options.h b/Common/Options.h
index 9ac991aafb..5a7cbc189f 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -531,7 +531,6 @@ double opt_mesh_draw_skin_only(OPT_ARGS_NUM);
double opt_mesh_save_all(OPT_ARGS_NUM);
double opt_mesh_save_element_tag_type(OPT_ARGS_NUM);
double opt_mesh_save_parametric(OPT_ARGS_NUM);
-double opt_mesh_save_distance(OPT_ARGS_NUM);
double opt_mesh_save_groups_of_nodes(OPT_ARGS_NUM);
double opt_mesh_color_carousel(OPT_ARGS_NUM);
double opt_mesh_zone_definition(OPT_ARGS_NUM);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index ea0ffc427f..ccb4f58a45 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -28,21 +28,6 @@
#include "discreteEdge.h"
#include "discreteFace.h"
#include "discreteRegion.h"
-#include "MElement.h"
-#include "GEdgeCompound.h"
-#include "GFaceCompound.h"
-#include "simpleFunction.h"
-#include "distanceTerm.h"
-#include "Context.h"
-#include "Numeric.h"
-#include "dofManager.h"
-#include "linearSystemGMM.h"
-#include "linearSystemCSR.h"
-#include "linearSystemFull.h"
-
-#if defined(HAVE_ANN)
-#include "ANN/ANN.h"
-#endif
void GModel::_storePhysicalTagsInEntities(int dim,
std::map<int, std::map<int, std::string> > &map)
@@ -858,320 +843,6 @@ int GModel::writePartitionedMSH(const std::string &baseName, bool binary,
return 1;
}
-int GModel::writeDistanceMSH(const std::string &name, double version, bool binary,
- bool saveAll, bool saveParametric, double scalingFactor)
-{
-
- std::vector<GEntity*> entities;
- getEntities(entities);
-
- 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_map;
- std::map<MVertex*,double > distance_map2;
- 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);
-
- int k=0;
- int maxDim = 0;
- for(unsigned int i = 0; i < entities.size(); i++){
- GEntity* ge = entities[i];
- maxDim = std::max(maxDim, ge->dim());
- 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));
- k++;
- }
- }
-
- //Compute geometrical distance to mesh boundaries
- //-------------------------------------------------
-
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() == (maxDim-1)) {
- 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]);
- }
- }
- }
- }
-
-// 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]);
-// }
-// }
-// }
-
-
- //Compute PDE for distance function
- //-------------------------------------------------
-
-#ifdef HAVE_TAUCS
- linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
-#else
- linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
- lsys->setNoisy(1);
- lsys->setGmres(1);
- lsys->setPrec(5.e-8);
-#endif
-
- 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());
- }
- }
- }
-
- 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;
- printf("L=%g \n", L);
-
- simpleFunction<double> DIFF(mu * mu), MONE(1.0);
- distanceTerm distance(this, 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);
- }
- groupOfElements g((GFace*)ge);
- distance.addToRightHandSide(myAssembler, g);
- }
- }
-
- 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 = std::min(0.9999, myAssembler.getDofValue(v, 0, 1));
- double dist = -mu * log(1. - value);
- itv->second = dist;
- }
-
- lsys->clear();
-
- //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");
- //GEntity* ge = entities[entities.size()-1];
- 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");
- }
- }
- }
- 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);
-
-
- return 1;
-
-}
-
int GModel::writePOS(const std::string &name, bool printElementary,
bool printElementNumber, bool printGamma, bool printEta,
bool printRho, bool printDisto, bool saveAll,
diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index 797ca59a77..12b32aa25e 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -23,6 +23,7 @@ set(SRC
Probe.cpp
HarmonicToTime.cpp ModulusPhase.cpp
HomologyComputation.cpp
+ Distance.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Plugin *.h)
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
new file mode 100644
index 0000000000..d4121bb6d6
--- /dev/null
+++ b/Plugin/Distance.cpp
@@ -0,0 +1,343 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+
+#include <stdlib.h>
+#include "Gmsh.h"
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "Distance.h"
+#include "simpleFunction.h"
+#include "distanceTerm.h"
+#include "Context.h"
+#include "Numeric.h"
+#include "dofManager.h"
+#include "linearSystemGMM.h"
+#include "linearSystemCSR.h"
+#include "linearSystemFull.h"
+
+#if defined(HAVE_ANN)
+#include "ANN/ANN.h"
+#endif
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterDistancePlugin()
+ {
+ return new GMSH_DistancePlugin();
+ }
+}
+
+std::string GMSH_DistancePlugin::getHelp() const
+{
+ return "Plugin(Distance) computes distances to boundaries in "
+ "a mesh.\n\n"
+ "Plugin(Distance) creates a bunch of files.";
+}
+
+PView *GMSH_DistancePlugin::execute(PView *v)
+{
+ std::vector<GEntity*> entities;
+ GModel::current()->getEntities(entities);
+
+ 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_map;
+ std::map<MVertex*,double > distance_map2;
+ 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);
+
+ int k = 0;
+ int maxDim = 0;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity* ge = entities[i];
+ maxDim = std::max(maxDim, ge->dim());
+ 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));
+ k++;
+ }
+ }
+
+ // Compute geometrical distance to mesh boundaries
+
+ for(unsigned int i = 0; i < entities.size(); i++){
+ if(entities[i]->dim() == (maxDim - 1)) {
+ 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]);
+ }
+ }
+ }
+ }
+
+// 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]);
+// }
+// }
+// }
+
+
+ // Compute PDE for distance function
+
+#ifdef HAVE_TAUCS
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#else
+ linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
+ lsys->setNoisy(1);
+ lsys->setGmres(1);
+ lsys->setPrec(5.e-8);
+#endif
+
+ 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());
+ }
+ }
+ }
+
+ 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;
+ printf("L=%g \n", L);
+
+ 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);
+ }
+ groupOfElements g((GFace*)ge);
+ distance.addToRightHandSide(myAssembler, g);
+ }
+ }
+
+ 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 = std::min(0.9999, myAssembler.getDofValue(v, 0, 1));
+ double dist = -mu * log(1. - value);
+ itv->second = dist;
+ }
+
+ lsys->clear();
+
+ //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");
+ }
+ }
+ }
+ 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);
+
+ return 0;
+}
diff --git a/Plugin/Distance.h b/Plugin/Distance.h
new file mode 100644
index 0000000000..35cd24d1d9
--- /dev/null
+++ b/Plugin/Distance.h
@@ -0,0 +1,33 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+
+#ifndef _DISTANCE_H_
+#define _DISTANCE_H_
+
+#include <string>
+#include "Plugin.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterDistancePlugin();
+}
+
+class GMSH_DistancePlugin : public GMSH_PostPlugin
+{
+ public:
+ GMSH_DistancePlugin(){}
+ std::string getName() const { return "Distance"; }
+ std::string getShortHelp() const
+ {
+ return "Compute distance to boundaries";
+ }
+ std::string getHelp() const;
+ std::string getAuthor() const { return "E. Marchandise"; }
+ PView *execute(PView *);
+};
+
+#endif
diff --git a/Plugin/Plugin.h b/Plugin/Plugin.h
index 455c1e204f..7ccf9d7339 100644
--- a/Plugin/Plugin.h
+++ b/Plugin/Plugin.h
@@ -80,6 +80,8 @@ class GMSH_PostPlugin : public GMSH_Plugin
{
public:
inline GMSH_PLUGIN_TYPE getType() const { return GMSH_Plugin::GMSH_POST_PLUGIN; }
+ virtual int getNbOptions() const { return 0; }
+ virtual StringXNumber *getOption(int iopt) { return 0; };
virtual int getNbOptionsStr() const { return 0; }
virtual StringXString *getOptionStr(int iopt) { return NULL; }
// run the plugin
@@ -107,10 +109,10 @@ class GMSH_PostPlugin : public GMSH_Plugin
class GMSH_SolverPlugin : public GMSH_Plugin
{
public:
- virtual int getNbOptionsStr() const { return 0; }
- virtual StringXString *getOptionStr(int iopt) { return 0; }
virtual int getNbOptions() const { return 0; }
virtual StringXNumber *getOption(int iopt) { return 0; };
+ virtual int getNbOptionsStr() const { return 0; }
+ virtual StringXString *getOptionStr(int iopt) { return 0; }
inline GMSH_PLUGIN_TYPE getType() const { return GMSH_Plugin::GMSH_SOLVER_PLUGIN; }
virtual void run() {} // do nothing
// popup dialog box
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index 17c2e423cd..2848acc40a 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -43,6 +43,7 @@
#include "Probe.h"
#include "GSHHS.h"
#include "HomologyComputation.h"
+#include "Distance.h"
// for testing purposes only :-)
#undef HAVE_DLOPEN
@@ -214,6 +215,8 @@ void PluginManager::registerDefaultPlugins()
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("Homology", GMSH_RegisterHomologyComputationPlugin()));
#endif
+ allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+ ("Distance", GMSH_RegisterDistancePlugin()));
}
#if defined(HAVE_FLTK)
--
GitLab