diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 3b178a9ab89177c357017dde987171749c2ca632..817d0ae4c207f07da3242c6089216a6030d31c8a 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -9,14 +9,11 @@
#include "Gmsh.h"
#include "GmshConfig.h"
#include "GModel.h"
-#include "MElement.h"
#include "Distance.h"
-#include "Context.h"
-#include "Numeric.h"
-
-#if defined(HAVE_SOLVER)
#include "simpleFunction.h"
#include "distanceTerm.h"
+#include "Context.h"
+#include "Numeric.h"
#include "dofManager.h"
#include "linearSystemGMM.h"
#include "linearSystemCSR.h"
@@ -24,7 +21,7 @@
#include "orthogonalTerm.h"
#include "laplaceTerm.h"
#include "crossConfTerm.h"
-#endif
+
StringXNumber DistanceOptions_Number[] = {
{GMSH_FULLRC, "PhysPoint", NULL, 0.},
@@ -40,6 +37,7 @@ StringXString DistanceOptions_String[] = {
{GMSH_FULLRC, "Filename", NULL, "distance.pos"}
};
+
extern "C"
{
GMSH_Plugin *GMSH_RegisterDistancePlugin()
@@ -61,21 +59,13 @@ std::string GMSH_DistancePlugin::getHelp() const
return "Plugin(Distance) computes distances to physical entities in "
"a mesh.\n\n"
- "Define the physical entities to which the distance is computed. "
- "If Point=0, Line=0, and Surface=0, then the distance is computed "
- "to all the boundaries of the mesh (edges in 2D and faces in 3D).\n\n"
+ "Define the physical entities to which the distance is computed. If Point=0, Line=0, and Surface=0, then the distance is computed to all the boundaries of the mesh (edges in 2D and faces in 3D).\n\n"
- "Computation<0. computes the geometrical euclidian distance "
- "(warning: different than the geodesic distance), and Computation=a>0.0 "
- "solves a PDE on the mesh with the diffusion constant mu = a*bbox, with "
- "bbox being the max size of the bounding box of the mesh (see paper "
- "Legrand 2006).\n\n"
+ "Computation<0. computes the geometrical euclidian distance (warning: different than the geodesic distance), and Computation=a>0.0 solves a PDE on the mesh with the diffusion constant mu = a*bbox, with bbox being the max size of the bounding box of the mesh (see paper Legrand 2006).\n\n"
- "Min Scale and max Scale, scale the distance function. If min Scale<0 "
- "and max Scale<0, then no scaling is applied to the distance function.\n\n"
+ "Min Scale and max Scale, scale the distance function. If min Scale<0 and max Scale<0, then no scaling is applied to the distance function.\n\n"
- "Plugin(Distance) creates a new distance view and also saves the view "
- "in the fileName.pos file.";
+ "Plugin(Distance) creates a new distance view and also saves the view in the fileName.pos file.";
}
int GMSH_DistancePlugin::getNbOptions() const
@@ -162,6 +152,7 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
PView *GMSH_DistancePlugin::execute(PView *v)
{
+
int id_pt = (int) DistanceOptions_Number[0].def;
int id_line = (int) DistanceOptions_Number[1].def;
int id_face = (int) DistanceOptions_Number[2].def;
@@ -170,9 +161,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
PView *view = new PView();
_data = getDataList(view);
-
-#if defined(HAVE_SOLVER)
-#if defined(HAVE_TAUCS)
+#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#else
linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
@@ -181,20 +170,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
lsys->setPrec(5.e-8);
#endif
dofManager<double> * dofView = new dofManager<double>(lsys);
-#endif
std::vector<GEntity*> _entities;
GModel::current()->getEntities(_entities);
GEntity* ge = _entities[_entities.size()-1];
int integrationPointTetra[2];
- if (type==-100){
- integrationPointTetra[0]=1;
- integrationPointTetra[1]=4;
- }
- else{
- integrationPointTetra[0]=0;
- integrationPointTetra[1]=0;
- }
+ integrationPointTetra[0]=0;
+ integrationPointTetra[1]=0;
int numnodes = 0;
for(unsigned int i = 0; i < _entities.size()-1; i++)
@@ -219,47 +201,19 @@ PView *GMSH_DistancePlugin::execute(PView *v)
distances.reserve(totNumNodes);
pt2Vertex.reserve(totNumNodes);
- std::map<MVertex*,double> _distanceE_map;
- std::map<MVertex*,int> _isInYarn_map;
- std::vector<int> index;
- std::vector<double> distancesE;
- std::vector<int> isInYarn;
- std::vector<SPoint3> closePts;
- std::vector<double> distances2;
- std::vector<double> distancesE2;
- std::vector<int> isInYarn2;
- std::vector<SPoint3> closePts2;
-
- if (type==-100){
- index.clear();
- distancesE.clear();
- closePts.clear();
- isInYarn.clear();
- isInYarn.reserve(totNumNodes);
- closePts.reserve(totNumNodes);
- distancesE.reserve(totNumNodes);
- index.reserve(totNumNodes);
- distances2.clear();
- distancesE2.clear();
- closePts2.clear();
- isInYarn2.clear();
- distances2.reserve(totNumNodes);
- isInYarn2.reserve(totNumNodes);
- closePts2.reserve(totNumNodes);
- distancesE2.reserve(totNumNodes);
- }
+ std::map<MVertex*,double> _distanceE_map;
+ std::map<MVertex*,int> _isInYarn_map;
+ std::vector<int> index;
+ std::vector<double> distancesE;
+ std::vector<int> isInYarn;
+ std::vector<SPoint3> closePts;
+ std::vector<double> distances2;
+ std::vector<double> distancesE2;
+ std::vector<int> isInYarn2;
+ std::vector<SPoint3> closePts2;
for (int i = 0; i < totNumNodes; i++) {
distances.push_back(1.e22);
- if (type==-100){
- distancesE.push_back(1.e22);
- isInYarn.push_back(0);
- closePts.push_back(SPoint3(0.,0.,0.));
- distances2.push_back(1.e22);
- distancesE2.push_back(1.e22);
- isInYarn2.push_back(0);
- closePts2.push_back(SPoint3(0.,0.,0.));
- }
}
int k = 0;
@@ -270,44 +224,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MVertex *v = ge->mesh_vertices[j];
pts.push_back(SPoint3(v->x(), v->y(),v->z()));
_distance_map.insert(std::make_pair(v, 0.0));
- if (type==-100){
- index.push_back(v->getIndex());
- _isInYarn_map.insert(std::make_pair(v, 0));
- _distanceE_map.insert(std::make_pair(v, 0.0));
- }
pt2Vertex[k] = v;
k++;
}
}
- if (type==-100){
- double jac[3][3];
- for (unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- IntPt *ptsi;
- int nptsi;
- double uvw[4][3];
- e->getIntegrationPoints(e->getPolynomialOrder(),&nptsi, &ptsi);
- for(int j = 0; j < 4; j++) {
- double xyz[3] = {e->getVertex(j)->x(),
- e->getVertex(j)->y(),
- e->getVertex(j)->z()};
- e->xyz2uvw(xyz, uvw[j]);
- }
-
- for(int ip = 0; ip < nptsi; ip++){
- const double u = ptsi[ip].pt[0];
- const double v = ptsi[ip].pt[1];
- const double w = ptsi[ip].pt[2];
- const double weight = ptsi[ip].weight;
- const double detJ = e->getJacobian(u, v, w, jac);
- SPoint3 p;
- e->pnt(u, v, w, p);
- pts.push_back(p);
- }
- }
- }
-
// Compute geometrical distance to mesh boundaries
//------------------------------------------------------
if (type < 0.0 ){
@@ -319,12 +240,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int gDim = g2->dim();
std::vector<int> phys = g2->getPhysicalEntities();
bool computeForEntity = false;
- for(unsigned int k = 0; k< phys.size(); k++){
+ for(int k = 0; k< phys.size(); k++){
int tagp = phys[k];
if (id_pt==0 && id_line==0 && id_face==0 && gDim==_maxDim-1 )
computeForEntity = true;
- else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) ||
- (tagp==id_face && gDim==2) )
+ else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) || (tagp==id_face && gDim==2) )
computeForEntity = true;
}
if (computeForEntity){
@@ -339,160 +259,31 @@ PView *GMSH_DistancePlugin::execute(PView *v)
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 && order==1) ||
- (e->getNumVertices() == 3 && order==2)){
- if (type==-100){
-// if ( !((p1.x()==p2.x()) & (p1.y()==p2.y()) & (p1.z()==p2.z())) ){
- signedDistancesPointsEllipseLine(iDistances, iDistancesE,
- iIsInYarn, iClosePts, pts, p1,p2);
-// }
- }
- else{
- signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
- }
+ if((e->getNumVertices() == 2 && order==1) || (e->getNumVertices() == 3 && order==2)){
+ signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
}
- else if(e->getNumVertices() == 3 && order==1){
+ else if((e->getNumVertices() == 3 and order==1) or (e->getNumVertices() == 6 and order==2)){
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 (type==-100){
- if( !((p1.x()==p2.x()) & (p1.y()==p2.y()) & (p1.z()==p2.z()))){
- if (iIsInYarn[kk]>0){
- if (isInYarn[kk]==0){
- distances[kk] = iDistances[kk];
- distancesE[kk]= iDistancesE[kk];
- isInYarn[kk] = iIsInYarn[kk];
- closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
- iClosePts[kk].z());
- }
- else{
- if (isInYarn[kk]!=iIsInYarn[kk]){
- if (isInYarn2[kk]==0){
- distances2[kk] = iDistances[kk];
- distancesE2[kk]= iDistancesE[kk];
- isInYarn2[kk] = iIsInYarn[kk];
- closePts2[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
- iClosePts[kk].z());
- }
- else{
- if (isInYarn2[kk]==iIsInYarn[kk]){
- if (iDistancesE[kk] < distancesE2[kk]){
- distances2[kk] = iDistances[kk];
- distancesE2[kk]= iDistancesE[kk];
- closePts2[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
- iClosePts[kk].z());
- }
- }
- }
- }
- else{
- if (iDistancesE[kk] < distancesE[kk]){
- distances[kk] = iDistances[kk];
- distancesE[kk]= iDistancesE[kk];
- closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
- iClosePts[kk].z());
- }
- }
- }
- }
- else{
- if (isInYarn[kk]==0){
- if (iDistancesE[kk] < distancesE[kk]){
- distances[kk] = iDistances[kk];
- distancesE[kk]= iDistancesE[kk];
- closePts[kk]= SPoint3(iClosePts[kk].x(),iClosePts[kk].y(),
- iClosePts[kk].z());
- }
- }
- }
+ if (std::abs(iDistances[kk]) < distances[kk]){
+ distances[kk] = std::abs(iDistances[kk]);
+ MVertex *v = pt2Vertex[kk];
+ _distance_map[v] = distances[kk];
}
- }
- else{
- if (fabs(iDistances[kk]) < distances[kk]){
- distances[kk] = fabs(iDistances[kk]);
- MVertex *v = pt2Vertex[kk];
- _distance_map[v] = distances[kk];
- }
- }
}
}
}
}
- if (type==-100){
- for (unsigned int kk = 0; kk< pts.size(); kk++) {
- if (isInYarn2[kk]>0){
- if (distancesE2[kk]>distancesE[kk]){
- distances[kk]=distances2[kk];
- distancesE[kk]=distancesE2[kk];
- isInYarn[kk]=isInYarn2[kk];
- closePts[kk]=closePts2[kk];
- }
- }
- if (kk<totNodes){
- MVertex *v = pt2Vertex[kk];
- _distance_map[v] = distancesE[kk];
- _distanceE_map[v] = distances[kk];
- _isInYarn_map[v] = isInYarn[kk];
- }
- }
- }
if (!existEntity){
if (id_pt != 0) Msg::Error("The Physical Point does not exist !");
if (id_line != 0) Msg::Error("The Physical Line does not exist !");
if (id_face != 0) Msg::Error("The Physical Surface does not exist !");
return view;
}
- printView(_entities, _distance_map);
- if (type==-100){
- Msg::Info("Writing integrationPointInYarn.pos");
- FILE* f5 = fopen("integrationPointInYarn.pos","w");
- FILE* f6 = fopen("integrationPointInYarn.bin","wb");
- FILE* f7 = fopen("integrationPointInYarn.txt","w");
- int j=0;
- fprintf(f5,"View \"integrationPointInYarn\"{\n");
- for (std::vector<int>::iterator it = isInYarn.begin(); it !=isInYarn.end(); it++) {
- if (j>=totNodes){
- int iPIY=*it;
- fwrite(&iPIY,sizeof(int),1,f6);
- fprintf(f7,"%d %lf %lf %lf\n",iPIY,pts[j].x(), pts[j].y(), pts[j].z());
- fprintf(f5,"SP(%g,%g,%g){%d};\n",
- pts[j].x(), pts[j].y(), pts[j].z(),
- *it);
- }
- j++;
- }
- fclose(f6);
- fclose(f7);
- fprintf(f5,"};\n");
- fclose(f5);
-
- Msg::Info("Writing isInYarn.pos");
- FILE * f4 = fopen("isInYarn.pos","w");
- fprintf(f4,"View \"isInYarn\"{\n");
- for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
- MElement *e = ge->getMeshElement(i);
- fprintf(f4,"SS(");
- std::vector<int> inYarn;
- 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*, int>::iterator it = _isInYarn_map.find(v);
- inYarn.push_back(it->second);
- }
- fprintf(f4,"){");
- for (unsigned int i=0; i<inYarn.size(); i++){
- if (i) fprintf(f4,",%d", inYarn[i]);
- else fprintf(f4,"%d", inYarn[i]);
- }
- fprintf(f4,"};\n");
- }
- fprintf(f4,"};\n");
- fclose(f4);
- }
-
+ printView(_entities, _distance_map);
}
// Compute PDE for distance function
@@ -509,12 +300,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int gDim = ge->dim();
bool fixForEntity = false;
std::vector<int> phys = ge->getPhysicalEntities();
- for(unsigned int k = 0; k< phys.size(); k++){
+ for(int k = 0; k< phys.size(); k++){
int tagp = phys[k];
if (id_pt==0 && id_line==0 && id_face==0 && gDim==_maxDim-1 )
fixForEntity = true;
- else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) ||
- (tagp==id_face && gDim==2) )
+ else if ( (tagp==id_pt && gDim==0)|| (tagp==id_line && gDim==1) || (tagp==id_face && gDim==2) )
fixForEntity = true;
}
if (fixForEntity){