Plugin/Distance: fix + pp

Plugin/Distance.cpp

@@ -54,7 +54,7 @@ GMSH_DistancePlugin::GMSH_DistancePlugin()
 
 std::string GMSH_DistancePlugin::getHelp() const
 {
-  return "Plugin(Distance) computes distances to physical entities in "
+  return "XXX 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 "
@@ -99,8 +99,7 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
   
   double minDist=1.e4;
   double maxDist=0.0;
-  for(std::map<MVertex*,double >::iterator itv = _distance_map.begin(); 
-      itv != _distance_map.end() ; ++itv){
+  for (std::map<MVertex*,double >::iterator itv=_distance_map.begin(); itv!=_distance_map.end(); ++itv){
     double dist = itv->second;
     if (dist>maxDist) maxDist = dist;
     if (dist<minDist) minDist = dist;
@@ -110,16 +109,18 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
   Msg::Info("Writing %s", _fileName.c_str());
   FILE *fName = fopen(_fileName.c_str(),"w");
   fprintf(fName, "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);
-
         int numNodes = e->getNumVertices();
-        if(e->getNumChildren()) numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
+        if (e->getNumChildren())
+          numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
         std::vector<double> x(numNodes), y(numNodes), z(numNodes);
         std::vector<double> *out = _data->incrementList(1, e->getType(), numNodes);
         std::vector<MVertex*> nods;
+        
         if (!e->getNumChildren())
           for(int i = 0; i < numNodes; i++)
             nods.push_back(e->getVertex(i));
@@ -127,20 +128,38 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
           for(int i = 0; i < e->getNumChildren(); i++)
             for(int j = 0; j < e->getChild(i)->getNumVertices(); j++)
               nods.push_back(e->getChild(i)->getVertex(j));
+        
         for (int nod=0; nod<numNodes; nod++) out->push_back((nods[nod])->x());
         for (int nod=0; nod<numNodes; nod++) out->push_back((nods[nod])->y());
         for (int nod=0; nod<numNodes; nod++) out->push_back((nods[nod])->z());
         
-	if (_maxDim == 3) fprintf(fName,"SS(");
-	else if (_maxDim == 2) fprintf(fName,"ST(");
+        if (_maxDim == 2)
+          switch (numNodes) {
+            case 2: fprintf(fName,"SL("); break;
+            case 3: fprintf(fName,"ST("); break;
+            case 4: fprintf(fName,"SQ("); break;
+            default: Msg::Error("Error in Plugin 'Distance' (numNodes=%d).",numNodes); break;
+          }
+        else if (_maxDim == 3)
+          switch (numNodes) {
+            case 4: fprintf(fName,"SS("); break;
+            case 8: fprintf(fName,"SH("); break;
+            case 6: fprintf(fName,"SI("); break;
+            case 5: fprintf(fName,"SY("); break;
+            default: Msg::Error("Error in Plugin 'Distance' (numNodes=%d).",numNodes); break;
+          }
+        
         std::vector<double> dist;
         for (int j=0; j<numNodes; j++) {
           MVertex *v =  nods[j];
-	  if(j) fprintf(fName,",%.16g,%.16g,%.16g",v->x(),v->y(), v->z());
-	  else fprintf(fName,"%.16g,%.16g,%.16g", v->x(),v->y(), v->z());
+          if (j)
+            fprintf(fName, ",%.16g,%.16g,%.16g", v->x(), v->y(), v->z());
+          else
+            fprintf(fName, "%.16g,%.16g,%.16g", v->x(), v->y(), v->z());
           std::map<MVertex*, double>::iterator it = _distance_map.find(v);
           dist.push_back(it->second);
         }
+        
         fprintf(fName,"){");
         for (unsigned int i=0; i<dist.size(); i++) {
           if (_minScale>0 && _maxScale>0)
@@ -148,14 +167,18 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
           else if (_minScale>0 && _maxScale<0)
             dist[i] = _minScale+dist[i];
           out->push_back(dist[i]);
-	  if (i) fprintf(fName,",%.16g", dist[i]);
-	  else fprintf(fName,"%.16g", dist[i]);
+          if (i)
+            fprintf(fName, ",%.16g", dist[i]);
+          else
+            fprintf(fName, "%.16g", dist[i]);
         }
         fprintf(fName,"};\n");
+        
       }
     }
   }
   fprintf(fName,"};\n");
+  
   fclose(fName);
 }
 
@@ -187,10 +210,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
     Msg::Error("This plugin needs a mesh !");
     return view;
   }
+  
   GEntity* ge = _entities[_entities.size()-1];
-  int integrationPointTetra[2];
-  integrationPointTetra[0]=0;
-  integrationPointTetra[1]=0;
+  int integrationPointTetra[2] = {0,0};
   
   int numnodes = 0;
   for (unsigned int i = 0; i < _entities.size()-1; i++)
@@ -213,11 +235,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
   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> isInYarn;
   std::vector<int> isInYarn2;
+  std::vector<SPoint3> closePts;
   std::vector<SPoint3> closePts2;
   
   for (int i=0; i<totNumNodes; i++) {
@@ -232,6 +254,10 @@ 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));
+/* TO DO (by AM)
+      SPoint3 p_empty();
+      _closePts_map.insert(std::make_pair(v, p_empty));
+*/
       pt2Vertex[k] = v;
       k++;
     }
@@ -242,6 +268,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
   if (type < 0.0 ) {
     
     bool existEntity = false;
+    
     for (unsigned int i=0; i<_entities.size(); i++) {
       GEntity* g2 = _entities[i];
       int tag = g2->tag();
@@ -280,6 +307,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
               distances[kk] = std::abs(iDistances[kk]);
               MVertex *v = pt2Vertex[kk];
               _distance_map[v] = distances[kk];
+/* TO DO (by AM)
+              _closePts_map[v] = iClosePts[kk];
+*/
             }
           }
         }
@@ -291,7 +321,12 @@ PView *GMSH_DistancePlugin::execute(PView *v)
       if (id_face != 0) Msg::Error("The Physical Surface does not exist !");
       return view;
     }
+    
     printView(_entities, _distance_map);
+    
+/* TO DO (by AM)
+    printView(_entities, _closePts_map);    
+*/
   }
   
   // Compute PDE for distance function
@@ -321,12 +356,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
           MElement *t = ge->getMeshElement(i);
           for (int k=0; k<t->getNumVertices(); k++) {
             MVertex *v = t->getVertex(k);
-	    dofView->fixVertex(v, 0, 1, 0.0);
+            dofView->fixVertex(v, 0, 1, 0.);
             bbox += SPoint3(v->x(), v->y(), v->z());
           }
         }
       }
     }
+    
     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 !");
@@ -341,12 +377,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
         for(unsigned int i=0; i<ge->getNumMeshElements(); ++i) {
           MElement *t = ge->getMeshElement(i);
           allElems.push_back(t);
-	  for(int k = 0; k < t->getNumVertices(); k++){
+          for (int k=0; k<t->getNumVertices(); k++)
             dofView->numberVertex(t->getVertex(k), 0, 1);
         }
       }
     }
-    }
     
     double L = norm(SVector3(bbox.max(), bbox.min())); 
     double mu = type*L;
@@ -422,6 +457,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
         }
       }
     }  
+    
     int mid = (int)floor(allElems.size() / 2.);
     MElement *e = allElems[mid];
     MVertex *vFIX = e->getVertex(0);
@@ -462,10 +498,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
       MElement *e = *it;
     
     int numNodes = e->getNumVertices();
-    if(e->getType() == TYPE_POLYG) numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
+    if (e->getType() == TYPE_POLYG)
+      numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
     std::vector<double> x(numNodes), y(numNodes), z(numNodes);
     std::vector<double> *out2 = data2->incrementList(1, e->getType(), numNodes);
     std::vector<MVertex*> nods;
+    std::vector<double> orth;
+    
     if(!e->getNumChildren())
       for(int i=0; i<numNodes; i++)
         nods.push_back(e->getVertex(i));
@@ -473,17 +512,33 @@ PView *GMSH_DistancePlugin::execute(PView *v)
       for(int i = 0; i < e->getNumChildren(); i++)
         for(int j = 0; j < e->getChild(i)->getNumVertices(); j++)
           nods.push_back(e->getChild(i)->getVertex(j));
+    
     for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->x()); 
     for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->y()); 
     for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->z()); 
     
-    if (_maxDim == 3) fprintf(f5,"SS(");
-    else if (_maxDim == 2) fprintf(f5,"ST(");
-    std::vector<double> orth;
+    if (_maxDim == 2)
+      switch (numNodes) {
+        case 2: fprintf(f5,"SL("); break;
+        case 3: fprintf(f5,"ST("); break;
+        case 4: fprintf(f5,"SQ("); break;
+        default: Msg::Fatal("Error in Plugin 'Distance' (numNodes=%g).",numNodes); break;
+      }
+    else if (_maxDim == 3)
+      switch (numNodes) {
+        case 4: fprintf(f5,"SS("); break;
+        case 8: fprintf(f5,"SH("); break;
+        case 6: fprintf(f5,"SI("); break;
+        case 5: fprintf(f5,"SY("); break;
+        default: Msg::Fatal("Error in Plugin 'Distance' (numNodes=%g).",numNodes); break;
+      }
+    
     for (int j=0; j<numNodes; j++) {
       MVertex *v =  nods[j];
-      if(j) fprintf(f5,",%g,%g,%g",v->x(),v->y(), v->z());
-      else fprintf(f5,"%g,%g,%g", v->x(),v->y(), v->z());
+      if (j)
+        fprintf(f5, ",%g,%g,%g", v->x(), v->y(), v->z());
+      else
+        fprintf(f5, "%g,%g,%g", v->x(), v->y(), v->z());
       double value;
       myAssembler.getDofValue(v, 0, 1, value);
       orth.push_back(value);
@@ -491,10 +546,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
     fprintf(f5,"){");
     for (unsigned int i=0; i<orth.size(); i++) {
       out2->push_back(orth[i]);
-      if (i) fprintf(f5,",%g", orth[i]);
-      else fprintf(f5,"%g", orth[i]);
+      if (i)
+        fprintf(f5,",%g", orth[i]);
+      else
+        fprintf(f5,"%g", orth[i]);
     }   
     fprintf(f5,"};\n");
+    
   }
   fprintf(f5,"};\n");
   fclose(f5);

Plugin/Distance.h

@@ -25,6 +25,7 @@ class GMSH_DistancePlugin : public GMSH_PostPlugin
   PViewDataList *_data;
  public:
   std::map<MVertex*,double > _distance_map;
+  std::map<MVertex*,SPoint3 > _closePts_map;
   GMSH_DistancePlugin(); 
   std::string getName() const { return "Distance"; }
   std::string getShortHelp() const