From a7a669c9d62c9a4f9107baf1c892ad3cfb29d831 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 6 Oct 2010 09:14:30 +0000
Subject: [PATCH]

---
 Fltk/statisticsWindow.cpp  | 143 +++++++++++++++++++++----------------
 Geo/GFaceCompound.cpp      |   4 +-
 Geo/GModelIO_Mesh.cpp      |  17 +++--
 Mesh/meshGEdge.cpp         |   9 +--
 Mesh/meshGFaceOptimize.cpp |  26 +++++--
 5 files changed, 117 insertions(+), 82 deletions(-)

diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 576202d9ad..beeb11201c 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -205,69 +205,90 @@ void statisticsWindow::compute(bool elementQuality)
   // meanAngle  = meanAngle / count;
   // printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
   //hack emi
-  /*
+
+  
   //Emi hack - MESH DEGREE VERTICES
-  std::vector<GEntity*> entities;
-  std::set<MEdge, Less_Edge> edges;
-  GModel::current()->getEntities(entities);
-  std::map<MVertex*, int > vert2Deg;
-  for(unsigned int i = 0; i < entities.size(); i++){
-    if(entities[i]->dim() < 2) continue;
-     for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
-      MElement *e =  entities[i]->getMeshElement(j);
-      for(unsigned int k = 0; k < e->getNumEdges(); k++){
-	edges.insert(e->getEdge(k));
-      }
-      for(unsigned int k = 0; k < e->getNumVertices(); k++){
-  	MVertex *v = e->getVertex(k);
-  	if (v->onWhat()->dim() < 2) continue; 
-  	std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
-  	if (it == vert2Deg.end()) {
-  	  vert2Deg.insert(std::make_pair(v,1));
-  	}
-  	else{
-  	  int nbE = it->second+1;
-  	  it->second = nbE;
-  	}
-      }
-    }
-  }
-  int dMin = 10;
-  int dMax = 0;
-  int d4 = 0;
-  int nbElems = vert2Deg.size();
-  std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
-  for(; itmap !=vert2Deg.end(); itmap++){
-    MVertex *v = itmap->first;
-    int nbE =  itmap->second;
-    dMin = std::min(nbE, dMin);
-    dMax = std::max(nbE, dMax);
-    if (nbE == 4) d4 += 1;
-  }
-  if (nbElems > 0)
-    printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
-  //end emi hack
-
-  FieldManager *fields = GModel::current()->getFields();
-  Field *f = fields->get(fields->background_field);
-  if(fields->background_field > 0){
-    std::set<MEdge, Less_Edge>::iterator it = edges.begin();
-    double sum = 0;
-    for (; it !=edges.end();++it){
-      MVertex *v0 = it->getVertex(0);
-      MVertex *v1 = it->getVertex(1);
-      double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
-		      (v0->y()-v1->y())*(v0->y()-v1->y())+
-		      (v0->z()-v1->z())*(v0->z()-v1->z()));
-      double lf =  (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
-      double e = (l>lf) ? lf/l : l/lf;
-      sum += e - 1.0;
-    }
-    double tau = exp ((1./edges.size()) * sum);
-    printf("nedegs = %d tau = %g\n",edges.size(),tau);
-  }
-  */
+  // std::vector<GEntity*> entities;
+  // std::set<MEdge, Less_Edge> edges;
+  // GModel::current()->getEntities(entities);
+  // std::map<MVertex*, int > vert2Deg;
+  // for(unsigned int i = 0; i < entities.size(); i++){
+  //   printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag());
+  //   if(entities[i]->dim() < 2 ) continue;
+  //   if(entities[i]->tag() != 100) continue;
+  //   printf("continuing \n");
+  //   for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+  //     MElement *e =  entities[i]->getMeshElement(j);
+  //     for(unsigned int k = 0; k < e->getNumEdges(); k++){
+  // 	edges.insert(e->getEdge(k));
+  //     }
+  //     for(unsigned int k = 0; k < e->getNumVertices(); k++){
+  // 	MVertex *v = e->getVertex(k);
+  // 	if (v->onWhat()->dim() < 2) continue; 
+  // 	std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+  // 	if (it == vert2Deg.end()) {
+  // 	  vert2Deg.insert(std::make_pair(v,1));
+  // 	}
+  // 	else{
+  // 	  int nbE = it->second+1;
+  // 	  it->second = nbE;
+  // 	}
+  //     }
+  //   }
+  // }
+  // int dMin = 10;
+  // int dMax = 0;
+  // int d4 = 0;
+  // int nbElems = vert2Deg.size();
+  // std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+  // for(; itmap !=vert2Deg.end(); itmap++){
+  //   MVertex *v = itmap->first;
+  //   int nbE =  itmap->second;
+  //   dMin = std::min(nbE, dMin);
+  //   dMax = std::max(nbE, dMax);
+  //   if (nbE == 4) d4 += 1;
+  // }
+  // if (nbElems > 0)
+  //   printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n", dMin, dMax, (double)d4/nbElems);
+
+  // FieldManager *fields = GModel::current()->getFields();
+  // Field *f = fields->get(fields->background_field);
+  // int nbEdges = edges.size();
+  // printf("nb edges =%d \n", nbEdges);
+  // system("rm qualEdges.txt");
+  // FILE *fp = fopen("qualEdges.txt", "w");
+  // std::vector<int> qualE;
+  // int nbS = 50;
+  // qualE.resize(nbS);
+  // if(fields->background_field > 0){
+  //   printf("found field \n");
+  //   std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+  //   double sum = 0;
+  //   for (; it !=edges.end();++it){
+  //     MVertex *v0 = it->getVertex(0);
+  //     MVertex *v1 = it->getVertex(1);
+  //     double l = sqrt((v0->x()-v1->x())*(v0->x()-v1->x())+
+  // 		      (v0->y()-v1->y())*(v0->y()-v1->y())+
+  // 		      (v0->z()-v1->z())*(v0->z()-v1->z()));
+  //     double lf =  (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),0.5*(v0->z()+v1->z()),v0->onWhat());
+  //     double el = l/lf;
+  //     int index = (int) ceil(el*nbS*0.5);
+  //     qualE[index]+= 1;
+  //     double e = (l>lf) ? lf/l : l/lf;
+  //     sum += e - 1.0;
+  //   }
+  //   double tau = exp ((1./edges.size()) * sum);
+  //   printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
 
+  //   double ibegin = 2./(2*nbS);
+  //   double inext = 2./nbS;
+  //   for (int i= 0; i< qualE.size(); i++){
+  //     fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+  //   }
+
+  // }
+  // fclose(fp);
+  //emi end hack
 
   int num = 0;
   static double s[50];
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 002f5232ff..c4c965d62d 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -671,8 +671,8 @@ bool GFaceCompound::parametrize() const
        //   Msg::Warning("!!! Overlap: parametrization switched to 'nonLIN conformal' map");
        //   noOverlap = parametrize_conformal_nonLinear() ;
     //}
-     if ( !noOverlap || !checkOrientation(0) ){
-       Msg::Warning("$$$ Overlap: parametrization switched to 'harmonic' map");
+     if (!noOverlap || !checkOrientation(0) ){
+       Msg::Warning("$$$ Flipping: parametrization switched to 'harmonic' map");
        parametrize(ITERU,HARMONIC); 
        parametrize(ITERV,HARMONIC);
      }
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index a7599223d2..e3e1e33b06 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -29,6 +29,7 @@
 #include "discreteFace.h"
 #include "discreteRegion.h"
 #include "MVertexPositionSet.h"
+#include "Os.h"
 
 #if defined(HAVE_POST)
 #include "PView.h"
@@ -59,15 +60,13 @@ void GModel::_storePhysicalTagsInEntities(int dim,
   }
 }
 
-static void replaceCommaByDot(const std::string name){
-  char myCommand[1000], myCommand2[1000];
-  sprintf(myCommand, "sed 's/,/./g' %s > temp.txt", name.c_str());
-  if(system(myCommand))
-    Msg::Error("sed command failed\n");
-  sprintf(myCommand2, "mv temp.txt %s ", name.c_str());
-  if(system(myCommand2))
-    Msg::Error("mv command failed\n");
-}
+ static void replaceCommaByDot(const std::string name){
+   char myCommand[1000], myCommand2[1000];
+   sprintf(myCommand, "sed 's/,/./g' %s > temp.txt", name.c_str());
+   SystemCall(myCommand);
+   sprintf(myCommand2, "mv temp.txt %s ", name.c_str());
+   SystemCall(myCommand2);
+ }
 
 static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map,
                         std::vector<MVertex*> &vertices)
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index d8144024d3..6462f33865 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -95,7 +95,7 @@ static double F_Lc(GEdge *ge, double t)
     lc_here = BGM_MeshSize(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
   else
     lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
-
+ 
   SVector3 der = ge->firstDer(t);
   //const double d = norm(der);
   //return d / lc_here; 
@@ -250,11 +250,12 @@ static double Integration(GEdge *ge, double t1, double t2,
   from.lc = f(ge, from.t);
   from.p = 0.0;
   Points.push_back(from);
-
+  
   to.t = t2;
   to.lc = f(ge, to.t);
+  
   RecursiveIntegration(ge, &from, &to, f, Points, Prec, &depth);
-
+  
   return Points.back().p;
 }
 
@@ -435,7 +436,7 @@ void meshGEdge::operator() (GEdge *ge)
     }
     mesh_vertices.resize(NUMP - 1);
   }
-  
+
   for(unsigned int i = 0; i < mesh_vertices.size() + 1; i++){
     MVertex *v0 = (i == 0) ?
       ge->getBeginVertex()->mesh_vertices[0] : mesh_vertices[i - 1];
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index cea0fc23a2..919b53e09c 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -17,6 +17,7 @@
 #include "GmshMessage.h"
 #include "Generator.h"
 #include "Context.h"
+#include "Os.h"
 
 #ifdef HAVE_MATCH
 extern "C" int FAILED_NODE;
@@ -1586,33 +1587,42 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
 
 void recombineIntoQuads(GFace *gf)
 {
-  gf->model()->writeMSH("before.msh");
+
+  double t1 = Cpu();
+
+  //gf->model()->writeMSH("before.msh");
   removeFourTrianglesNodes(gf, false);
   int success = _recombineIntoQuads(gf,0);
-  gf->model()->writeMSH("raw.msh");
+  //gf->model()->writeMSH("raw.msh");
   for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++) 
     laplaceSmoothing(gf);
+
+  //blossom-quad algo
   if(success && CTX::instance()->mesh.algoRecombine == 1){    
-    gf->model()->writeMSH("smoothed.msh");
+    //gf->model()->writeMSH("smoothed.msh");
     int COUNT = 0;
     char NAME[256];
     while(1){
       int x = removeTwoQuadsNodes(gf);
-      if (x)sprintf(NAME,"iter%dTQ.msh",COUNT++);
+      //if (x)sprintf(NAME,"iter%dTQ.msh",COUNT++);
       if (x)gf->model()->writeMSH(NAME);
       int y= removeDiamonds(gf);
-      if (y)sprintf(NAME,"iter%dD.msh",COUNT++);
+      //if (y)sprintf(NAME,"iter%dD.msh",COUNT++);
       if (y)gf->model()->writeMSH(NAME);
       int  z = _edgeSwapQuadsForBetterQuality ( gf );
-      if (z)sprintf(NAME,"iter%dS.msh",COUNT++);
+      //if (z)sprintf(NAME,"iter%dS.msh",COUNT++);
       if (z)gf->model()->writeMSH(NAME);      
       if (!(x+y+z))break;
     }
     for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++){ 
       laplaceSmoothing(gf);
     }
+    double t2 = Cpu();
+    Msg::Info("Blossom recombination algorithm completed (%g s)", t2 - t1);
     return;
   }
+
+  //simple recombination algo
   _recombineIntoQuads(gf,0);
   for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++) 
     laplaceSmoothing(gf);
@@ -1620,4 +1630,8 @@ void recombineIntoQuads(GFace *gf)
   for(int i = 0; i < CTX::instance()->mesh.nbSmoothing; i++) 
     laplaceSmoothing(gf);
   //  gf->model()->writeMSH("after.msh");
+
+  double t2 = Cpu();
+  Msg::Info("Simple recombination algorithm completed (%g s)", t2 - t1);
+
 }
-- 
GitLab