diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index ec095ad4ecafee29e0ad5b663f9ed3223536efec..e560be4ccb7975efedf6212a8bd1fc23854d00f4 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -2514,8 +2514,8 @@ optionWindow::optionWindow(int deltaFontSize)
       mesh.value[5]->callback(mesh_options_ok_cb);
 
       static Fl_Menu_Item menu_quality_type[] = {
+        {"SICN", 0, 0, 0},
         {"Gamma", 0, 0, 0},
-        {"Eta", 0, 0, 0},
         {"Rho", 0, 0, 0},
         {"Disto", 0, 0, 0},
         {0}
diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 0feddfb744c59b17dea9ac307d526297369281c6..f9bb15bdd94c8406ee57d56cb8e6604e4d3cb949 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -18,6 +18,8 @@
 #include "OS.h"
 #include "Field.h"
 
+enum QM_HISTO {QMH_SICN_XY, QMH_SICN_3D, QMH_GAMMA_XY, QMH_GAMMA_3D, QMH_RHO_XY, QMH_RHO_3D};
+
 void statistics_cb(Fl_Widget *w, void *data)
 {
   FlGui::instance()->stats->show();
@@ -30,64 +32,53 @@ static void statistics_update_cb(Fl_Widget *w, void *data)
 
 static void statistics_histogram_cb(Fl_Widget *w, void *data)
 {
-  std::string name((const char*)data);
+  QM_HISTO qmh = *(QM_HISTO*)data;
 
   std::vector<double> x, y;
 
-  if(name == "Gamma2D"){
+  if (qmh == QMH_SICN_XY) {
     for(int i = 0; i < 100; i++){
-      x.push_back((double)i / 99);
+      x.push_back((double)(2*i-99) / 99);
       y.push_back(FlGui::instance()->stats->quality[0][i]);
     }
-    new PView("Gamma", "# Elements", x, y);
+    new PView("SICN", "# Elements", x, y);
   }
-  else if(name == "Eta2D"){
+  else if (qmh == QMH_GAMMA_XY) {
     for(int i = 0; i < 100; i++){
       x.push_back((double)i / 99);
       y.push_back(FlGui::instance()->stats->quality[1][i]);
     }
-    new PView("Eta", "# Elements", x, y);
+    new PView("Gamma", "# Elements", x, y);
   }
-  else if(name == "Rho2D"){
+  else if (qmh == QMH_RHO_XY) {
     for(int i = 0; i < 100; i++){
       x.push_back((double)i / 99);
       y.push_back(FlGui::instance()->stats->quality[2][i]);
     }
     new PView("Rho", "# Elements", x, y);
   }
-  else if(name == "Disto2D"){
-    for(int i = 0; i < 100; i++){
-      x.push_back((double)i / 99);
-      y.push_back(FlGui::instance()->stats->quality[3][i]);
-    }
-    new PView("Disto", "# Elements", x, y);
-  }
-  else{
+  else {
     std::vector<GEntity*> entities_;
     GModel::current()->getEntities(entities_);
     std::map<int, std::vector<double> > d;
-    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);
-	if(name == "Gamma3D")
-	  d[e->getNum()].push_back(e->gammaShapeMeasure());
-	else if(name == "Eta3D")
-	  d[e->getNum()].push_back(e->etaShapeMeasure());
-	else if(name == "Rho3D")
-#ifdef METRICSHAPEMEASURE
-	  if (e->getDim() == 3)
-	    d[e->getNum()].push_back(e->metricShapeMeasure());
-	  else
-	    d[e->getNum()].push_back(1);
-#else
-    d[e->getNum()].push_back(e->rhoShapeMeasure());
-#endif
-	else
-	  d[e->getNum()].push_back(e->distoShapeMeasure());
+    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);
+        if (qmh == QMH_SICN_3D) {
+          double minSICN, maxSICN;
+          e->invCondNumRange(minSICN, maxSICN);
+          d[e->getNum()].push_back(minSICN);
+        }
+        else if (qmh == QMH_GAMMA_3D)
+          d[e->getNum()].push_back(e->gammaShapeMeasure());
+        else if (qmh == QMH_RHO_3D)
+          d[e->getNum()].push_back(e->rhoShapeMeasure());
       }
     }
-    name.resize(name.size() - 2);
+    std::string name = (qmh == QMH_SICN_3D) ? "SICN" :
+                       (qmh == QMH_GAMMA_3D) ? "Gamma" :
+                       (qmh == QMH_RHO_3D) ? "Rho" : "";
     new PView(name, "ElementData", GModel::current(), d);
   }
 
@@ -101,7 +92,7 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
 
   int num = 0;
   int width = 26 * FL_NORMAL_SIZE;
-  int height = 5 * WB + 18 * BH;
+  int height = 5 * WB + 17 * BH;
 
   win = new paletteWindow
     (width, height, CTX::instance()->nonModalWindows ? true : false, "Statistics");
@@ -135,16 +126,14 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
       value[num++] = new Fl_Output(2 * WB, 2 * WB + 11 * BH, IW, BH, "Time for 2D mesh");
       value[num++] = new Fl_Output(2 * WB, 2 * WB + 12 * BH, IW, BH, "Time for 3D mesh");
 
-      value[num] = new Fl_Output(2 * WB, 2 * WB + 13 * BH, IW, BH, "Gamma");
-      value[num]->tooltip("~ inscribed_radius / circumscribed_radius"); num++;
-      value[num] = new Fl_Output(2 * WB, 2 * WB + 14 * BH, IW, BH, "Eta");
-      value[num]->tooltip("~ volume^(2/3) / sum_edge_length^2"); num++;
+      value[num] = new Fl_Output(2 * WB, 2 * WB + 13 * BH, IW, BH, "SICN");
+      value[num]->tooltip("~ signed inverse condition number"); num++;
+      value[num] = new Fl_Output(2 * WB, 2 * WB + 14 * BH, IW, BH, "Gamma");
+      value[num]->tooltip("~ inscribed_radius / circumscribed_radius (simplices)"); num++;
       value[num] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Rho");
       value[num]->tooltip("~ min_edge_length / max_edge_length"); num++;
-      value[num] = new Fl_Output(2 * WB, 2 * WB + 16 * BH, IW, BH, "Disto");
-      value[num]->tooltip("~ min (J_min/J_0, J_0/J_max)"); num++;
 
-      for(int i = 0; i < 4; i++){
+      for(int i = 0; i < 3; i++){
         int ww = 3 * FL_NORMAL_SIZE;
         new Fl_Box
           (FL_NO_BOX, width - 3 * ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "Plot");
@@ -153,14 +142,14 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
         butt[2 * i + 1] = new Fl_Button
           (width - ww - 2 * WB, 2 * WB + (13 + i) * BH, ww, BH, "3D");
       }
-      butt[0]->callback(statistics_histogram_cb, (void *)"Gamma2D");
-      butt[1]->callback(statistics_histogram_cb, (void *)"Gamma3D");
-      butt[2]->callback(statistics_histogram_cb, (void *)"Eta2D");
-      butt[3]->callback(statistics_histogram_cb, (void *)"Eta3D");
-      butt[4]->callback(statistics_histogram_cb, (void *)"Rho2D");
-      butt[5]->callback(statistics_histogram_cb, (void *)"Rho3D");
-      butt[6]->callback(statistics_histogram_cb, (void *)"Disto2D");
-      butt[7]->callback(statistics_histogram_cb, (void *)"Disto3D");
+      static const QM_HISTO qmh0 = QMH_SICN_XY, qmh1 = QMH_SICN_3D, qmh2 = QMH_GAMMA_XY,
+                            qmh3 = QMH_GAMMA_3D, qmh4 = QMH_RHO_XY, qmh5 = QMH_RHO_3D;
+      butt[0]->callback(statistics_histogram_cb, (void*) &qmh0);
+      butt[1]->callback(statistics_histogram_cb, (void*) &qmh1);
+      butt[2]->callback(statistics_histogram_cb, (void*) &qmh2);
+      butt[3]->callback(statistics_histogram_cb, (void*) &qmh3);
+      butt[4]->callback(statistics_histogram_cb, (void*) &qmh4);
+      butt[5]->callback(statistics_histogram_cb, (void*) &qmh5);
 
       group[1]->end();
     }
@@ -204,188 +193,6 @@ statisticsWindow::statisticsWindow(int deltaFontSize)
 
 void statisticsWindow::compute(bool elementQuality)
 {
-#if 0
-  {
-    // MINIMUM MAXIMUM ANGLES
-    double minAngle = 1.0; //M_PI;
-    double meanAngle = 0.0;
-    int count = 0;
-    std::vector<GEntity*> entities;
-    GModel::current()->getEntities(entities);
-    std::map<int, std::vector<double> > d;
-    for(unsigned int i = 0; i < entities.size(); i++){
-      if(entities[i]->dim() == 3) {// continue;//<3
-	for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
-	  MElement *e = entities[i]->getMeshElement(j);
-	  double angle = e->angleShapeMeasure();
-	  minAngle = std::min(minAngle, angle);
-	  meanAngle += angle;
-	  count++;
-	}
-      }
-    }
-    meanAngle  = meanAngle / count;
-    printf("Angles min =%g av=%g nbhex=%d\n", minAngle, meanAngle, count);
-  }
-
-  // {
-  //   // 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;
-  //     // if(entities[i]->tag() < 100) 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);
-  //   FieldManager *fields = GModel::current()->getFields();
-  //   Field *f = fields->get(fields->background_field);
-  //   int nbEdges = edges.size();
-  //   printf("nb edges =%d \n", nbEdges);
-  //   if(system("rm qualEdges.txt"));
-  //   FILE *fp = Fopen("qualEdges.txt", "w");
-  //   std::vector<int> qualE;
-  //   int nbS = 50;
-  //   qualE.resize(nbS);
-  //   if(fields.getBackgroundField() > 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 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);
-  // }
-  // {
-  //   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;
-  //     if(entities[i]->tag() != 10) 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);
-  //   FieldManager *fields = GModel::current()->getFields();
-  //   Field *f = fields->get(fields.getBackgroundField());
-  //   int nbEdges = edges.size();
-  //   if(system("rm qualEdges.txt"));
-  //   FILE *fp = Fopen("qualEdges.txt", "w");
-  //   std::vector<int> qualE;
-  //   int nbS = 50;
-  //   qualE.resize(nbS);
-  //   if(fields.getBackgroundField() > 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 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);
-  //     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);
-  // }
-#endif
-
   int num = 0;
   static double s[50];
   static char label[50][256];
@@ -418,10 +225,7 @@ void statisticsWindow::compute(bool elementQuality)
   sprintf(label[num], "%g", s[15]); value[num]->value(label[num]); num++;
 
   if(!elementQuality){
-    for(int i = 0; i < 8; i += 2) butt[i]->deactivate();
-    sprintf(label[num], "Press Update");
-    value[num]->deactivate();
-    value[num]->value(label[num]); num++;
+    for(int i = 0; i < 6; i += 2) butt[i]->deactivate();
     sprintf(label[num], "Press Update");
     value[num]->deactivate();
     value[num]->value(label[num]); num++;
@@ -433,7 +237,7 @@ void statisticsWindow::compute(bool elementQuality)
     value[num]->value(label[num]); num++;
   }
   else{
-    for(int i = 0; i < 8; i += 2) butt[i]->activate();
+    for(int i = 0; i < 6; i += 2) butt[i]->activate();
     sprintf(label[num], "%.4g (%.4g->%.4g)", s[17], s[18], s[19]);
     value[num]->activate();
     value[num]->value(label[num]); num++;
@@ -443,9 +247,6 @@ void statisticsWindow::compute(bool elementQuality)
     sprintf(label[num], "%.4g (%.4g->%.4g)", s[23], s[24], s[25]);
     value[num]->activate();
     value[num]->value(label[num]); num++;
-    sprintf(label[num], "%.4g (%.4g->%.4g)", s[46], s[47], s[48]);
-    value[num]->activate();
-    value[num]->value(label[num]); num++;
   }
 
   // post
diff --git a/Fltk/statisticsWindow.h b/Fltk/statisticsWindow.h
index 7d748df7b32af9d2294713e80bc890650123d63f..bc800b72b101a386c53fcfd3ced4584b8c79efcb 100644
--- a/Fltk/statisticsWindow.h
+++ b/Fltk/statisticsWindow.h
@@ -18,7 +18,7 @@ class statisticsWindow{
   Fl_Button *butt[8];
   Fl_Group *group[3];
   Fl_Box *memUsage;
-  double quality[4][100];
+  double quality[3][100];
  public:
   statisticsWindow(int deltaFontSize);
   void compute(bool elementQuality);
diff --git a/Geo/GModelVertexArrays.cpp b/Geo/GModelVertexArrays.cpp
index bae77544602b4d41370fe0d9db7929ca109ecaf4..abc00806fa9ab43a5ddb9e154aa7ef53fa1b3a54 100644
--- a/Geo/GModelVertexArrays.cpp
+++ b/Geo/GModelVertexArrays.cpp
@@ -104,15 +104,13 @@ bool isElementVisible(MElement *ele)
     if(CTX::instance()->mesh.qualityType == 3)
       q = ele->distoShapeMeasure();
     else if(CTX::instance()->mesh.qualityType == 2)
-#ifdef METRICSHAPEMEASURE
-      q = ele->metricShapeMeasure();
-#else
       q = ele->rhoShapeMeasure();
-#endif
     else if(CTX::instance()->mesh.qualityType == 1)
-      q = ele->etaShapeMeasure();
-    else
       q = ele->gammaShapeMeasure();
+    else {
+      double sICNMax;
+      ele->invCondNumRange(q, sICNMax);
+    }
     if(q < CTX::instance()->mesh.qualityInf ||
        q > CTX::instance()->mesh.qualitySup) return false;
   }
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 1caf16832c63064ac8c95dd4b68f1f8c2ad54edb..ad0aa24ba587b637746e5a0c32b2b9998bcc2452 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -49,52 +49,33 @@
 template<class T>
 static void GetQualityMeasure(std::vector<T*> &ele,
                               double &gamma, double &gammaMin, double &gammaMax,
-                              double &eta, double &etaMin, double &etaMax,
+                              double &minSICN, double &minSICNMin, double &minSICNMax,
                               double &rho, double &rhoMin, double &rhoMax,
-                              double &disto, double &distoMin, double &distoMax,
-                              double &scaledJacMin, double &scaledJacMax,
-                              double quality[4][100])
+                              double quality[3][100])
 {
   for(unsigned int i = 0; i < ele.size(); i++){
     double g = ele[i]->gammaShapeMeasure();
     gamma += g;
     gammaMin = std::min(gammaMin, g);
     gammaMax = std::max(gammaMax, g);
-    double e = ele[i]->etaShapeMeasure();
-    eta += e;
-    etaMin = std::min(etaMin, e);
-    etaMax = std::max(etaMax, e);
-#ifdef METRICSHAPEMEASURE
-    static int aa = 0;
-    if (++aa == 1) Msg::Warning("Computing metric shape measure instead of rho shape measure !");
-    double r = ele[i]->metricShapeMeasure();
-    rho += r;
-    rhoMin = std::min(rhoMin, r);
-    rhoMax = std::max(rhoMax, r);
-#else
+    double s, sDum;
+    ele[i]->invCondNumRange(s, sDum);
+    minSICN += s;
+    minSICNMin = std::min(minSICNMin, s);
+    minSICNMax = std::max(minSICNMax, s);
     double r = ele[i]->rhoShapeMeasure();
     rho += r;
     rhoMin = std::min(rhoMin, r);
     rhoMax = std::max(rhoMax, r);
-#endif
-    double jmin,jmax; ele[i]->scaledJacRange(jmin,jmax);
-    double d = jmin;
-    disto += d;
-    distoMin = std::min(distoMin, d);
-    distoMax = std::max(distoMax, d);
-    scaledJacMin = std::min(scaledJacMin, jmin);
-    scaledJacMax = std::max(scaledJacMax, jmax);
-
     for(int j = 0; j < 100; j++){
-      if(g > j / 100. && g <= (j + 1) / 100.) quality[0][j]++;
-      if(e > j / 100. && e <= (j + 1) / 100.) quality[1][j]++;
+      if(s > (2*j-100) / 100. && s <= (2*j-98) / 100.) quality[0][j]++;
+      if(g > j / 100. && g <= (j + 1) / 100.) quality[1][j]++;
       if(r > j / 100. && r <= (j + 1) / 100.) quality[2][j]++;
-      if(d > j / 100. && d <= (j + 1) / 100.) quality[3][j]++;
     }
   }
 }
 
-void GetStatistics(double stat[50], double quality[4][100])
+void GetStatistics(double stat[50], double quality[3][100])
 {
   for(int i = 0; i < 50; i++) stat[i] = 0.;
 
@@ -139,45 +120,36 @@ void GetStatistics(double stat[50], double quality[4][100])
     for(int i = 0; i < 3; i++)
       for(int j = 0; j < 100; j++)
         quality[i][j] = 0.;
+    double minSICN = 0., minSICNMin = 1., minSICNMax = -1.;
     double gamma = 0., gammaMin = 1., gammaMax = 0.;
-    double eta = 0., etaMin = 1., etaMax = 0.;
     double rho = 0., rhoMin = 1., rhoMax = 0.;
-    double disto = 0., distoMin=1., distoMax = 0.;
-    double jmin = 1.e22, jmax = -1.e22;
 
     double N = stat[9] + stat[10] + stat[11] + stat[12];
     if(N){ // if we have 3D elements
       for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
         GetQualityMeasure((*it)->tetrahedra, gamma, gammaMin, gammaMax,
-                          eta, etaMin, etaMax, rho, rhoMin, rhoMax,
-                          disto, distoMin, distoMax, jmin, jmax, quality);
+                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
         GetQualityMeasure((*it)->hexahedra, gamma, gammaMin, gammaMax,
-                          eta, etaMin, etaMax, rho, rhoMin, rhoMax,
-                          disto, distoMin, distoMax, jmin, jmax, quality);
+                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
         GetQualityMeasure((*it)->prisms, gamma, gammaMin, gammaMax,
-                          eta, etaMin, etaMax, rho, rhoMin, rhoMax,
-                          disto, distoMin, distoMax,jmin, jmax, quality);
+                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
         GetQualityMeasure((*it)->pyramids, gamma, gammaMin, gammaMax,
-                          eta, etaMin, etaMax, rho, rhoMin, rhoMax,
-                          disto, distoMin, distoMax, jmin, jmax, quality);
+                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
       }
     }
     else{ // 2D elements
       N = stat[7] + stat[8];
       for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
         GetQualityMeasure((*it)->quadrangles, gamma, gammaMin, gammaMax,
-                          eta, etaMin, etaMax, rho, rhoMin, rhoMax,
-                          disto, distoMin, distoMax, jmin, jmax, quality);
+                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
         GetQualityMeasure((*it)->triangles, gamma, gammaMin, gammaMax,
-                          eta, etaMin, etaMax, rho, rhoMin, rhoMax,
-                          disto, distoMin, distoMax, jmin, jmax, quality);
+                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
       }
     }
     if(N){
-      stat[17] = gamma / N ;  stat[18] = gammaMin;  stat[19] = gammaMax;
-      stat[20] = eta / N ;    stat[21] = etaMin;    stat[22] = etaMax;
-      stat[23] = rho / N ;    stat[24] = rhoMin;    stat[25] = rhoMax;
-      stat[46] = disto / N ;  stat[47] = distoMin;  stat[48] = distoMax;
+      stat[17] = minSICN / N; stat[18] = minSICNMin; stat[19] = minSICNMax;
+      stat[20] = gamma / N;   stat[21] = gammaMin;   stat[22] = gammaMax;
+      stat[23] = rho / N;     stat[24] = rhoMin;     stat[25] = rhoMax;
     }
   }