diff --git a/Common/GmshMessage.cpp b/Common/GmshMessage.cpp
index 2dde6762d6c0bbf30641a104a2572261307504ec..260266ee4fa8b8d0cc54ae593eaa54d9fe3c3d2c 100644
--- a/Common/GmshMessage.cpp
+++ b/Common/GmshMessage.cpp
@@ -360,7 +360,9 @@ void Msg::Info(const char *fmt, ...)
   if(_client) _client->Info(str);
 
 #if defined(HAVE_FLTK)
+#if defined(_OPENMP)
   #pragma omp critical
+#endif
   {
     if(FlGui::available()){
       FlGui::instance()->check();
@@ -394,7 +396,9 @@ void Msg::Direct(const char *fmt, ...)
   if(_client) _client->Info(str);
 
 #if defined(HAVE_FLTK)
+#if defined(_OPENMP)
 #pragma omp master
+#endif
   {
     if(FlGui::available()){
       FlGui::instance()->check();
@@ -428,7 +432,9 @@ void Msg::StatusBar(bool log, const char *fmt, ...)
   if(_client && log) _client->Info(str);
 
 #if defined(HAVE_FLTK)
+#if defined(_OPENMP)
 #pragma omp master
+#endif
   {
     if(FlGui::available()){
       if(log) FlGui::instance()->check();
diff --git a/Common/OS.cpp b/Common/OS.cpp
index 3202ced94920ee9a73694599b7660ce08328f5ae..027d62e54144d95d04b4eda36b1ac7e660c79e2c 100644
--- a/Common/OS.cpp
+++ b/Common/OS.cpp
@@ -269,7 +269,7 @@ int SystemCall(const std::string &command, bool blocking)
   }
   if(!blocking) cmd += " &";
   Msg::Info("Calling '%s'", cmd.c_str());
-  system(cmd.c_str());
+  if(!system(cmd.c_str())) return 1;
   return 0;
 #endif
 }
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index c1011296183c08842bf37eabbe08870bcc3d7a7f..8ae57ddb43a46b611fa34e8813946f60e60007d8 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -297,7 +297,7 @@ class GEntity {
   // get the number of mesh elements (total and by type) in the entity
   virtual unsigned int getNumMeshElements() { return 0; }
   virtual unsigned int getNumMeshParentElements() { return 0; }
-  virtual void getNumMeshElements(unsigned *const c) const { };
+  virtual void getNumMeshElements(unsigned *const c) const { }
 
   // get the start of the array of a type of element
   virtual MElement *const *getStartElementType(int type) const { return 0; }
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index a77663bbab5314173710b1161cf1c0d07dfce0da..c9fa0ed3d9dc5eb7daa8b2bb8bc7cf3252787d0e 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -51,11 +51,11 @@ void writeMSHPeriodicNodes (FILE *fp, std::vector<GEntity*> &entities)
 void readMSHPeriodicNodes (FILE *fp, GModel *gm)
 {
   int count;
-  fscanf(fp, "%d", &count);
+  if(fscanf(fp, "%d", &count) != 1) return;
   for(int i = 0; i < count; i++){
-    int dim,slave,master;
-    fscanf(fp, "%d %d %d", &dim, &slave, &master);
-    GEntity *s=0,*m=0;
+    int dim, slave, master;
+    if(fscanf(fp, "%d %d %d", &dim, &slave, &master) != 3) continue;
+    GEntity *s = 0, *m = 0;
     switch(dim){
     case 0 : s = gm->getVertexByTag(slave); m = gm->getVertexByTag(master); break;
     case 1 : s = gm->getEdgeByTag(slave); m = gm->getEdgeByTag(master); break;
@@ -64,10 +64,10 @@ void readMSHPeriodicNodes (FILE *fp, GModel *gm)
     if (s && m){
       s->setMeshMaster(m->tag());
       int numv;
-      fscanf(fp, "%d", &numv);
+      if(fscanf(fp, "%d", &numv) != 1) numv = 0;
       for(int j = 0; j < numv; j++){
-	int v1,v2;
-	fscanf(fp,"%d %d", &v1, &v2);
+	int v1, v2;
+	if(fscanf(fp,"%d %d", &v1, &v2) != 2) continue;
 	MVertex *mv1 = gm->getMeshVertexByTag(v1);
 	MVertex *mv2 = gm->getMeshVertexByTag(v2);
 	s->correspondingVertices[mv1] = mv2;
@@ -87,7 +87,7 @@ int GModel::readMSH(const std::string &name)
   char str[256] = "";
 
   // detect prehistoric MSH files
-  fgets(str, sizeof(str), fp);
+  if(!fgets(str, sizeof(str), fp)){ fclose(fp); return 0; }
   if(!strncmp(&str[1], "NOD", 3) || !strncmp(&str[1], "NOE", 3)){
     fclose(fp);
     return _readMSH2(name);
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index a615b2a7248ba6726aab6197e44fc9de567001be..e9c6ae75fb7467cf570758d1a3726a3683e05e72 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2287,6 +2287,9 @@ int Extrude_ProtudePoint(int type, int ip,
     }
     c->end = chapeau;
     break;
+//  case ANALYTICAL:
+//    
+//    break;
   default:
     Msg::Error("Unknown extrusion type");
     return pv->Num;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6dc4adb8accd679991140a64a9f80940700d038e..9fa590258b94c8f9745d326f6cecc5145f419642 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -30,7 +30,9 @@ double MElement::_isInsideTolerance = 1.e-6;
 
 MElement::MElement(int num, int part) : _visible(1)
 {
-#pragma omp critical
+#if defined(_OPENMP)
+  #pragma omp critical
+#endif
   {
     // we should make GModel a mandatory argument to the constructor
     GModel *m = GModel::current();
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index bb40e19b81189af31863d1ff28345337bd959673..d06aa37af52c54529f9c461882d4f8cba5ee4d80 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -45,7 +45,9 @@ double angle3Vertices(const MVertex *p1, const MVertex *p2, const MVertex *p3)
 MVertex::MVertex(double x, double y, double z, GEntity *ge, int num)
   : _visible(1), _order(1), _x(x), _y(y), _z(z), _ge(ge)
 {
+#if defined(_OPENMP)
 #pragma omp critical
+#endif
   {
     // we should make GModel a mandatory argument to the constructor
     GModel *m = GModel::current();
@@ -71,7 +73,9 @@ void MVertex::deleteLast()
 
 void MVertex::forceNum(int num)
 {
+#if defined(_OPENMP)
 #pragma omp critical
+#endif
   {
     _num = num;
     GModel::current()->setMaxVertexNumber
diff --git a/Graphics/drawContext.cpp b/Graphics/drawContext.cpp
index b0f3a0a5b80631e375a084fe9038a43a631bb7a3..e456beef83be6e9faae638dccc4884a8b3c89167 100644
--- a/Graphics/drawContext.cpp
+++ b/Graphics/drawContext.cpp
@@ -836,7 +836,7 @@ bool drawContext::select(int type, bool multiple, bool mesh,
     //   for triangle, 4 for quad) and the fourth is the index of the element in
     //   the vertex array
     GLuint names = *ptr++;
-    *ptr++; // mindepth
+    ptr++; // mindepth
     GLuint maxdepth = *ptr++;
     if(names == 2){
       GLuint depth = maxdepth;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index cf18547702e5f814f81f22344a8890bda8927ca0..a07b8d964dd78b12958991cc5259bddcaf2fd2ff 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -480,19 +480,25 @@ static void Mesh2D(GModel *m)
     while(1){
       int nPending = 0;
       std::vector<GFace*> _temp; _temp.insert(_temp.begin(), f.begin(), f.end());
+#if defined(_OPENMP)
 #pragma omp parallel for schedule (dynamic)
+#endif
       for(size_t K = 0 ; K < _temp.size() ; K++){
 	if (_temp[K]->meshStatistics.status == GFace::PENDING){
 	  meshGFace mesher (true, CTX::instance()->mesh.multiplePasses);
 	  mesher(_temp[K]);
+#if defined(_OPENMP)
 #pragma omp critical
+#endif
 	  {
 	    nPending++;
 	  }
 	}
         if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
       }
+#if defined(_OPENMP)
 #pragma omp master
+#endif
       for(std::set<GFace*>::iterator it = cf.begin(); it != cf.end(); ++it){
         if ((*it)->meshStatistics.status == GFace::PENDING){
 	  meshGFace mesher (true, CTX::instance()->mesh.multiplePasses);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index e7c5e29075fb2805f2f36943601cd386fe6c9510..e7ff6d65bf84ce38d6b9f283a2f8b5b3a976a123 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -2257,7 +2257,7 @@ int untangleInvalidQuads(GFace *gf, int niter)
   return N;
 }
 
-static int orientationOK (GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3)
+/*static int orientationOK (GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3)
 {
   SPoint2 p1, p2, p3;
   reparamMeshVertexOnFace(v1, gf, p1);
@@ -2265,7 +2265,7 @@ static int orientationOK (GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3)
   reparamMeshVertexOnFace(v3, gf, p3);
   if (robustPredicates::orient2d(p1, p2, p3) < 0) return true;
   return false;
-}
+}*/
 
 static int allowSwap (GFace *gf, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4)
 {
diff --git a/Mesh/periodical.cpp b/Mesh/periodical.cpp
index ffb29169c0cc4b23705354db2559153a5fa65794..4e26735f5a1eda9552a1c21a212b591a8d2aee52 100644
--- a/Mesh/periodical.cpp
+++ b/Mesh/periodical.cpp
@@ -51,7 +51,7 @@ int geo_cell::search_line(std::pair<int,int> line){
 
   for(i=0;i<lines.size();i++){
     if(lines[i].first==line.first && lines[i].second==line.second) return i;
-	if(lines[i].first==line.second && lines[i].second==line.first) return i;
+    if(lines[i].first==line.second && lines[i].second==line.first) return i;
   }
 
   return -1;
@@ -71,9 +71,9 @@ void voroMetal3D::execute(double h){
   for(it=model->firstRegion();it!=model->lastRegion();it++)
   {
     gr = *it;
-	if(gr->getNumMeshElements()>0){
-	  execute(gr,h);
-	}
+    if(gr->getNumMeshElements()>0){
+      execute(gr,h);
+    }
   }
 }
 
@@ -90,19 +90,19 @@ void voroMetal3D::execute(GRegion* gr,double h){
   vertices2.clear();
   radii.clear();
   vertices.clear();
-	
+
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	
-	for(j=0;j<element->getNumVertices();j++){
-	  vertex = element->getVertex(j);
-	  vertices.insert(vertex);
-	}
+
+    for(j=0;j<element->getNumVertices();j++){
+      vertex = element->getVertex(j);
+      vertices.insert(vertex);
+    }
   }
 
   for(it=vertices.begin();it!=vertices.end();it++){
     vertices2.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
-	radii.push_back(1.0);
+    radii.push_back(1.0);
   }
 
   execute(vertices2,radii,0,h);
@@ -112,15 +112,15 @@ void voroMetal3D::execute(std::vector<double>& properties,int radical,double h){
   unsigned int i;
   std::vector<SPoint3> vertices;
   std::vector<double> radii;
-  
+
   vertices.clear();
   radii.clear();
-	
-  for(i=0;i<properties.size()/4;i++){	  
-	vertices.push_back(SPoint3(properties[4*i],properties[4*i+1],properties[4*i+2]));
-	radii.push_back(properties[4*i+3]);
+
+  for(i=0;i<properties.size()/4;i++){
+    vertices.push_back(SPoint3(properties[4*i],properties[4*i+1],properties[4*i+2]));
+    radii.push_back(properties[4*i+3]);
   }
-  
+
   execute(vertices,radii,radical,h);
 }
 
@@ -164,11 +164,11 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
   max_z = -1000000000.0;
   for(i=0;i<vertices.size();i++){
     min_x = std::min(vertices[i].x(),min_x);
-	max_x = std::max(vertices[i].x(),max_x);
-	min_y = std::min(vertices[i].y(),min_y);
-	max_y = std::max(vertices[i].y(),max_y);
-	min_z = std::min(vertices[i].z(),min_z);
-	max_z = std::max(vertices[i].z(),max_z);
+    max_x = std::max(vertices[i].x(),max_x);
+    min_y = std::min(vertices[i].y(),min_y);
+    max_y = std::max(vertices[i].y(),max_y);
+    min_z = std::min(vertices[i].z(),min_z);
+    max_z = std::max(vertices[i].z(),max_z);
   }
 
   delta = 0.2*(max_x - min_x);
@@ -182,56 +182,56 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
   for(i=0;i<vertices.size();i++){
     if(radical==0){
       contA.put(i,vertices[i].x(),vertices[i].y(),vertices[i].z());
-	}
-	else{
-	  contB.put(i,vertices[i].x(),vertices[i].y(),vertices[i].z(),radii[i]);
-	}
+    }
+    else{
+      contB.put(i,vertices[i].x(),vertices[i].y(),vertices[i].z(),radii[i]);
+    }
   }
 
   c_loop_all loopA(contA);
   c_loop_all loopB(contB);
-	
+
   if(radical==0){
     loopA.start();
     do{
-	  contA.compute_cell(cell,loopA);
-	  loopA.pos(x,y,z);
-	  pointer = new voronoicell_neighbor();
-	  *pointer = cell;
-	  pointers.push_back(pointer);
-	  generators.push_back(SPoint3(x,y,z));
+      contA.compute_cell(cell,loopA);
+      loopA.pos(x,y,z);
+      pointer = new voronoicell_neighbor();
+      *pointer = cell;
+      pointers.push_back(pointer);
+      generators.push_back(SPoint3(x,y,z));
     }while(loopA.inc());
   }
   else{
     loopB.start();
-	do{
-	  contB.compute_cell(cell,loopB);
-	  loopB.pos(x,y,z);
-	  pointer = new voronoicell_neighbor();
-	  *pointer = cell;
-	  pointers.push_back(pointer);
-	  generators.push_back(SPoint3(x,y,z));
-	}while(loopB.inc());
+    do{
+      contB.compute_cell(cell,loopB);
+      loopB.pos(x,y,z);
+      pointer = new voronoicell_neighbor();
+      *pointer = cell;
+      pointers.push_back(pointer);
+      generators.push_back(SPoint3(x,y,z));
+    }while(loopB.inc());
   }
 
   initialize_counter();
 
   min_area = 1000000000.0;
-	
+
   for(i=0;i<pointers.size();i++){
     areas.clear();
-	
+
     pointers[i]->face_areas(areas);
-	  
+
     for(j=0;j<areas.size();j++){
-	  if(areas[j]<min_area){
-	    min_area = areas[j];
-	  }
-	}
+      if(areas[j]<min_area){
+        min_area = areas[j];
+      }
+    }
   }
-	
+
   printf("\nSquared root of smallest face area : %.9f\n\n",sqrt(min_area));
-		
+
   std::ofstream file("MicrostructurePolycrystal3D.pos");
   file << "View \"test\" {\n";
 
@@ -239,121 +239,121 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
   std::ofstream file5("SET.map");
   file2 << "c=" << h << ";\n";
   int countPeriodSurf=0;
-  int countVolume=0;	
+  int countVolume=0;
   for(i=0;i<pointers.size();i++){
-	obj = geo_cell();
+    obj = geo_cell();
 
-	faces.clear();
-	voronoi_vertices.clear();
-	pointers[i]->face_vertices(faces);
-	pointers[i]->vertices(generators[i].x(),generators[i].y(),generators[i].z(),voronoi_vertices);
+    faces.clear();
+    voronoi_vertices.clear();
+    pointers[i]->face_vertices(faces);
+    pointers[i]->vertices(generators[i].x(),generators[i].y(),generators[i].z(),voronoi_vertices);
 
-	obj.line_loops.resize(pointers[i]->number_of_faces());
-	obj.orientations.resize(pointers[i]->number_of_faces());
+    obj.line_loops.resize(pointers[i]->number_of_faces());
+    obj.orientations.resize(pointers[i]->number_of_faces());
 
-	face_number = 0;
-	end = 0;
+    face_number = 0;
+    end = 0;
     while(end<faces.size()){
-	  start = end + 1;
-	  end = start + faces[end];
-	  for(j=start;j<end;j++){
-		if(j<end-1){
-	      index1 = faces[j];
-		  index2 = faces[j+1];
-		}
-		else{
-		  index1 = faces[end-1];
-		  index2 = faces[start];
-		}
-		x1 = voronoi_vertices[3*index1];
-		y1 = voronoi_vertices[3*index1+1];
-		z1 = voronoi_vertices[3*index1+2];
-		x2 = voronoi_vertices[3*index2];
-		y2 = voronoi_vertices[3*index2+1];
-		z2 = voronoi_vertices[3*index2+2];
-		print_segment(SPoint3(x1,y1,z1),SPoint3(x2,y2,z2),file);
-
-		val = obj.search_line(std::pair<int,int>(index1,index2));
-		if(val==-1){
-		  obj.lines.push_back(std::pair<int,int>(index1,index2));
-		  obj.line_loops[face_number].push_back(obj.lines.size()-1);
-		  val = obj.lines.size()-1;
-		}
-		else{
-		  obj.line_loops[face_number].push_back(val);
-		}
-
-		last = obj.line_loops[face_number].size()-1;
-		if(last==0){
-	      obj.orientations[face_number].push_back(0);
-		}
-		else if(obj.lines[obj.line_loops[face_number][last-1]].second==obj.lines[val].first){
-		  obj.orientations[face_number][last-1] = 0;
-		  obj.orientations[face_number].push_back(0);
-		}
-		else if(obj.lines[obj.line_loops[face_number][last-1]].first==obj.lines[val].first){
-		  obj.orientations[face_number][last-1] = 1;
-		  obj.orientations[face_number].push_back(0);
-		}
-		else if(obj.lines[obj.line_loops[face_number][last-1]].second==obj.lines[val].second){
-		  obj.orientations[face_number][last-1] = 0;
-		  obj.orientations[face_number].push_back(1);
-		}
-		else{
-		  obj.orientations[face_number][last-1] = 1;
-		  obj.orientations[face_number].push_back(1);
-		}
-	  }
-
-	  face_number++;
-	}
-
-	for(j=0;j<voronoi_vertices.size()/3;j++){
-	  print_geo_point(get_counter(),voronoi_vertices[3*j],voronoi_vertices[3*j+1],voronoi_vertices[3*j+2],file2);
-	  obj.points2.push_back(get_counter());
-	  increase_counter();
-	}
-
-	for(j=0;j<obj.lines.size();j++){
-	  print_geo_line(get_counter(),obj.points2[obj.lines[j].first],obj.points2[obj.lines[j].second],file2);
-	  obj.lines2.push_back(get_counter());
-	  increase_counter();
-	}
-
-	for(j=0;j<obj.line_loops.size();j++){
-	  temp.clear();
-	  temp2.clear();
-	  for(k=0;k<obj.line_loops[j].size();k++){
-	    temp.push_back(obj.lines2[obj.line_loops[j][k]]);
-		temp2.push_back(obj.orientations[j][k]);
-	  }
-	  print_geo_line_loop(get_counter(),temp,temp2,file2);
-	  obj.line_loops2.push_back(get_counter());
-	  increase_counter();
-	}
-
-	for(j=0;j<obj.line_loops2.size();j++){
-	  print_geo_face(get_counter(),obj.line_loops2[j],file2);
-	  countPeriodSurf++;
-	  file5 <<get_counter()<< "\t" <<"SURFACE"<<get_counter()<<"\t"<<"NSET\n";
-	  obj.faces2.push_back(get_counter());
-	  increase_counter();
-	}
-	  
-	print_geo_face_loop(get_counter(),obj.faces2,file2);
-	obj.face_loops2 = get_counter();
-	increase_counter();
-
-	print_geo_volume(get_counter(),obj.face_loops2,file2);
-	mem = get_counter();
-	countVolume++;
-	file5 <<get_counter()<< "\t" <<"GRAIN"<<countVolume<<"\t"<<"ELSET\n";
-        increase_counter();
-        print_geo_physical_volume(get_counter(),mem,file2);
-	increase_counter();
-  }
-
-  file2 << "Coherence;\n";	
+      start = end + 1;
+      end = start + faces[end];
+      for(j=start;j<end;j++){
+        if(j<end-1){
+          index1 = faces[j];
+          index2 = faces[j+1];
+        }
+        else{
+          index1 = faces[end-1];
+          index2 = faces[start];
+        }
+        x1 = voronoi_vertices[3*index1];
+        y1 = voronoi_vertices[3*index1+1];
+        z1 = voronoi_vertices[3*index1+2];
+        x2 = voronoi_vertices[3*index2];
+        y2 = voronoi_vertices[3*index2+1];
+        z2 = voronoi_vertices[3*index2+2];
+        print_segment(SPoint3(x1,y1,z1),SPoint3(x2,y2,z2),file);
+
+        val = obj.search_line(std::pair<int,int>(index1,index2));
+        if(val==-1){
+          obj.lines.push_back(std::pair<int,int>(index1,index2));
+          obj.line_loops[face_number].push_back(obj.lines.size()-1);
+          val = obj.lines.size()-1;
+        }
+        else{
+          obj.line_loops[face_number].push_back(val);
+        }
+
+        last = obj.line_loops[face_number].size()-1;
+        if(last==0){
+          obj.orientations[face_number].push_back(0);
+        }
+        else if(obj.lines[obj.line_loops[face_number][last-1]].second==obj.lines[val].first){
+          obj.orientations[face_number][last-1] = 0;
+          obj.orientations[face_number].push_back(0);
+        }
+        else if(obj.lines[obj.line_loops[face_number][last-1]].first==obj.lines[val].first){
+          obj.orientations[face_number][last-1] = 1;
+          obj.orientations[face_number].push_back(0);
+        }
+        else if(obj.lines[obj.line_loops[face_number][last-1]].second==obj.lines[val].second){
+          obj.orientations[face_number][last-1] = 0;
+          obj.orientations[face_number].push_back(1);
+        }
+        else{
+          obj.orientations[face_number][last-1] = 1;
+          obj.orientations[face_number].push_back(1);
+        }
+      }
+
+      face_number++;
+    }
+
+    for(j=0;j<voronoi_vertices.size()/3;j++){
+      print_geo_point(get_counter(),voronoi_vertices[3*j],voronoi_vertices[3*j+1],voronoi_vertices[3*j+2],file2);
+      obj.points2.push_back(get_counter());
+      increase_counter();
+    }
+
+    for(j=0;j<obj.lines.size();j++){
+      print_geo_line(get_counter(),obj.points2[obj.lines[j].first],obj.points2[obj.lines[j].second],file2);
+      obj.lines2.push_back(get_counter());
+      increase_counter();
+    }
+
+    for(j=0;j<obj.line_loops.size();j++){
+      temp.clear();
+      temp2.clear();
+      for(k=0;k<obj.line_loops[j].size();k++){
+        temp.push_back(obj.lines2[obj.line_loops[j][k]]);
+        temp2.push_back(obj.orientations[j][k]);
+      }
+      print_geo_line_loop(get_counter(),temp,temp2,file2);
+      obj.line_loops2.push_back(get_counter());
+      increase_counter();
+    }
+
+    for(j=0;j<obj.line_loops2.size();j++){
+      print_geo_face(get_counter(),obj.line_loops2[j],file2);
+      countPeriodSurf++;
+      file5 <<get_counter()<< "\t" <<"SURFACE"<<get_counter()<<"\t"<<"NSET\n";
+      obj.faces2.push_back(get_counter());
+      increase_counter();
+    }
+
+    print_geo_face_loop(get_counter(),obj.faces2,file2);
+    obj.face_loops2 = get_counter();
+    increase_counter();
+
+    print_geo_volume(get_counter(),obj.face_loops2,file2);
+    mem = get_counter();
+    countVolume++;
+    file5 <<get_counter()<< "\t" <<"GRAIN"<<countVolume<<"\t"<<"ELSET\n";
+    increase_counter();
+    print_geo_physical_volume(get_counter(),mem,file2);
+    increase_counter();
+  }
+
+  file2 << "Coherence;\n";
   file << "};\n";
 
   for(i=0;i<pointers.size();i++) delete pointers[i];
@@ -421,9 +421,9 @@ void voroMetal3D::print_geo_line_loop(int index,std::vector<int>& indices,std::v
   file << "Line Loop(" << index << ")={";
 
   for(i=0;i<indices.size();i++){
-	if(orientations[i]==1) file << "-";
+    if(orientations[i]==1) file << "-";
     file << indices[i];
-	if(i<indices.size()-1) file << ",";
+    if(i<indices.size()-1) file << ",";
   }
 
   file << "};\n";
@@ -436,7 +436,7 @@ void voroMetal3D::print_geo_face_loop(int index,std::vector<int>& indices,std::o
 
   for(i=0;i<indices.size();i++){
     file << indices[i];
-	if(i<indices.size()-1) file << ",";
+    if(i<indices.size()-1) file << ",";
   }
 
   file << "};\n";
@@ -447,7 +447,7 @@ void voroMetal3D::correspondance(double e){
   unsigned int j;
   int count;
   int val;
-  int normal;
+  //int normal;
   int phase;
   bool flag;
   bool flag1;
@@ -492,387 +492,387 @@ void voroMetal3D::correspondance(double e){
   std::list<int>::iterator it9;
   std::list<int>::iterator it10;
   std::list<GEdge*>::iterator mem;
-	
-  faces.clear();	
-	
+
+  faces.clear();
+
   for(it=model->firstFace();it!=model->lastFace();it++)
   {
     gf = *it;
-	if(gf->numRegions()==1){
-	  faces.push_back(gf);
-	}
+    if(gf->numRegions()==1){
+      faces.push_back(gf);
+    }
   }
-	
+
   centers.clear();
   markings.clear();
   pairs.clear();
   categories.clear();
-	
+
   for(i=0;i<faces.size();i++){
-	x = 0.0;
-	y = 0.0;
-	z = 0.0;
-	
-	vertices.clear();
-	
-	vertices = faces[i]->vertices();
-	
-	for(it2=vertices.begin();it2!=vertices.end();it2++){
-	  x = x + (*it2)->x();
-	  y = y + (*it2)->y();
-	  z = z + (*it2)->z();
+    x = 0.0;
+    y = 0.0;
+    z = 0.0;
+
+    vertices.clear();
+
+    vertices = faces[i]->vertices();
+
+    for(it2=vertices.begin();it2!=vertices.end();it2++){
+      x = x + (*it2)->x();
+      y = y + (*it2)->y();
+      z = z + (*it2)->z();
     }
-	  
-	x = x/vertices.size();
-	y = y/vertices.size();
-	z = z/vertices.size();
-	  
-	centers.insert(std::pair<GFace*,SPoint3>(faces[i],SPoint3(x,y,z)));
-  }	
-	
+
+    x = x/vertices.size();
+    y = y/vertices.size();
+    z = z/vertices.size();
+
+    centers.insert(std::pair<GFace*,SPoint3>(faces[i],SPoint3(x,y,z)));
+  }
+
   for(i=0;i<faces.size();i++){
     markings.insert(std::pair<GFace*,bool>(faces[i],0));
   }
-	
+
   count = 0;
-	
+
   std::ofstream file("MicrostructurePolycrystal3D.pos");
   file << "View \"test\" {\n";
-  
+
   std::ofstream file2("PERIODIC.map");
-	
+
   for(i=0;i<faces.size();i++){
     for(j=0;j<faces.size();j++){
-	  it3 = centers.find(faces[i]);
-	  it4 = centers.find(faces[j]);
-		
-	  p1 = it3->second;
-	  p2 = it4->second;
-		
-	  delta_x = fabs(p2.x()-p1.x());
-	  delta_y = fabs(p2.y()-p1.y());
-	  delta_z = fabs(p2.z()-p1.z());
-		
-	  flag = correspondance(delta_x,delta_y,delta_z,e,val);
-		
-	  if(flag){
-	    it5 = markings.find(faces[i]);
-	    it6 = markings.find(faces[j]);
-		
-		if(it5->second==0 && it6->second==0){
-		  it5->second = 1;
-		  it6->second = 1;
-		  
-		  pairs.push_back(std::pair<GFace*,GFace*>(faces[i],faces[j]));
-		  categories.push_back(val);
-			
-		  print_segment(p1,p2,file);
-		  
-		  //file2 << "PERIODIC\tSURFACE"<<faces[i]->tag() << "\tSURFACE" << faces[j]->tag() << "\t" << p2.x()-p1.x() << "\t" << p2.y()-p1.y() << "\t" << p2.z()-p1.z() << "\n";	
-		  if (abs((p2.x()-p1.x()-1.0))<0.0001){
-			if (abs((p2.y()-p1.y()))<0.0001){
-				if (abs((p2.z()-p1.z()))<0.0001){
-					file2 << "NSET\tFRONT = FRONT + SURFACE"<<faces[j]->tag()<<"\n";
-					file2 << "NSET\tBACK = BACK + SURFACE"<<faces[i]->tag()<<"\n";
-				}else if(abs((p2.z()-p1.z()-1.0))<0.0001){
-					file2 << "NSET\tFRONTTOP = FRONTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-					file2 << "NSET\tBACKBOTTOM = BACKBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-				}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-					file2 << "NSET\tFRONTBOTTOM = FRONTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKTOP = BACKTOP + SURFACE"<<faces[i]->tag()<<"\n";
-				}
-			}else if (abs((p2.y()-p1.y()-1.0))<0.0001){
-				if (abs((p2.z()-p1.z()))<0.0001){
-                                        file2 << "NSET\tFRONTRIGHT = FRONTRIGHT + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKLEFT = BACKLEFT + SURFACE"<<faces[i]->tag()<<"\n";
-                                }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
-				        file2 << "NSET\tFRONTRIGHTTOP = FRONTRIGHTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKLEFTBOTTOM = BACKLEFTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-                                }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTRIGHTBOTTOM = FRONTRIGHTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKLEFTTOP = BACKLEFTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-				}
-			}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
-				if (abs((p2.z()-p1.z()))<0.0001){
-                                        file2 << "NSET\tFRONTLEFT = FRONTLEFT + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKRIGHT = BACKRIGHT + SURFACE"<<faces[i]->tag()<<"\n";
-                                }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTLEFTTOP = FRONTLEFTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKRIGHTBOTTOM = BACKRIGHTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-                                }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTLEFTBOTTOM = FRONTLEFTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tBACKRIGHTTOP = BACKRIGHTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-                                }
-			}
-		  }else if (abs((p1.x()-p2.x()-1.0))<0.0001){
-			if (abs((p2.y()-p1.y()))<0.0001){
-                                if (abs((p2.z()-p1.z()))<0.0001){
-                                        file2 << "NSET\tFRONT = FRONT + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACK = BACK + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if(abs((p2.z()-p1.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTBOTTOM = FRONTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKTOP = BACKTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTTOP = FRONTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKBOTTOM = BACKBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                }
-                        }else if (abs((p2.y()-p1.y()-1.0))<0.0001){
-                                if (abs((p2.z()-p1.z()))<0.0001){
-                                        file2 << "NSET\tFRONTLEFT = FRONTLEFT + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKRIGHT = BACKRIGHT + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTLEFTBOTTOM = FRONTLEFTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKRIGHTTOP = BACKRIGHTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTLEFTTOP = FRONTLEFTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKRIGHTBOTTOM = BACKRIGHTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                }
-                        }else if (abs((p1.y()-p2.y()-1.0))<0.0001){
-                                if (abs((p2.z()-p1.z()))<0.0001){
-                                        file2 << "NSET\tFRONTRIGHT = FRONTRIGHT + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKLEFT = BACKLEFT + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTRIGHTBOTTOM = FRONTRIGHTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKLEFTTOP = BACKLEFTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-                                        file2 << "NSET\tFRONTRIGHTTOP = FRONTRIGHTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tBACKLEFTBOTTOM = BACKLEFTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                }
-                        }
-		  }else if (abs((p1.x()-p2.x()))<0.0001){
-			  if (abs((p2.y()-p1.y()-1.0))<0.0001){
-				if (abs((p2.z()-p1.z()))<0.0001){
-        	                	file2 << "NSET\tRIGHT = RIGHT + SURFACE"<<faces[j]->tag()<<"\n";
-                	        	file2 << "NSET\tLEFT = LEFT + SURFACE"<<faces[i]->tag()<<"\n";
-				}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
-					file2 << "NSET\tRIGHTTOP = RIGHTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tLEFTBOTTOM = LEFTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-				}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-					file2 << "NSET\tRIGHTBOTTOM = RIGHTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                        file2 << "NSET\tLEFTTOP = LEFTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-				}
-	                  }else if (abs((p1.y()-p2.y()-1.0))<0.0001){
-				if (abs((p2.z()-p1.z()))<0.0001){
-                                        file2 << "NSET\tRIGHT = RIGHT + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tLEFT = LEFT + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
-                                        file2 << "NSET\tRIGHTBOTTOM = RIGHTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tLEFTTOP = LEFTTOP + SURFACE"<<faces[j]->tag()<<"\n";
-                                }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-                                        file2 << "NSET\tRIGHTTOP = RIGHTTOP + SURFACE"<<faces[i]->tag()<<"\n";
-                                        file2 << "NSET\tLEFTBOTTOM = LEFTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                                }
-	                  }else if (abs((p1.y()-p2.y()))<0.0001){
-			  	if (abs((p2.z()-p1.z()-1.0))<0.0001){
-                	        	file2 << "NSET\tTOP = TOP + SURFACE"<<faces[j]->tag()<<"\n";
-                        		file2 << "NSET\tBOTTOM = BOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
-	                  	}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
-        	                	file2 << "NSET\tTOP = TOP + SURFACE"<<faces[i]->tag()<<"\n";
-                	        	file2 << "NSET\tBOTTOM = BOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
-                  	  	}
-			  }
-		  }
-		  count++;
-		}
-	  }
-	}
+      it3 = centers.find(faces[i]);
+      it4 = centers.find(faces[j]);
+
+      p1 = it3->second;
+      p2 = it4->second;
+
+      delta_x = fabs(p2.x()-p1.x());
+      delta_y = fabs(p2.y()-p1.y());
+      delta_z = fabs(p2.z()-p1.z());
+
+      flag = correspondance(delta_x,delta_y,delta_z,e,val);
+
+      if(flag){
+        it5 = markings.find(faces[i]);
+        it6 = markings.find(faces[j]);
+
+        if(it5->second==0 && it6->second==0){
+          it5->second = 1;
+          it6->second = 1;
+
+          pairs.push_back(std::pair<GFace*,GFace*>(faces[i],faces[j]));
+          categories.push_back(val);
+
+          print_segment(p1,p2,file);
+
+          //file2 << "PERIODIC\tSURFACE"<<faces[i]->tag() << "\tSURFACE" << faces[j]->tag() << "\t" << p2.x()-p1.x() << "\t" << p2.y()-p1.y() << "\t" << p2.z()-p1.z() << "\n";
+          if (abs((p2.x()-p1.x()-1.0))<0.0001){
+            if (abs((p2.y()-p1.y()))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tFRONT = FRONT + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACK = BACK + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if(abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTTOP = FRONTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKBOTTOM = BACKBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTBOTTOM = FRONTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKTOP = BACKTOP + SURFACE"<<faces[i]->tag()<<"\n";
+              }
+            }else if (abs((p2.y()-p1.y()-1.0))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tFRONTRIGHT = FRONTRIGHT + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKLEFT = BACKLEFT + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTRIGHTTOP = FRONTRIGHTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKLEFTBOTTOM = BACKLEFTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTRIGHTBOTTOM = FRONTRIGHTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKLEFTTOP = BACKLEFTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+              }
+            }else if (abs((p1.y()-p2.y()-1.0))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tFRONTLEFT = FRONTLEFT + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKRIGHT = BACKRIGHT + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTLEFTTOP = FRONTLEFTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKRIGHTBOTTOM = BACKRIGHTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTLEFTBOTTOM = FRONTLEFTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBACKRIGHTTOP = BACKRIGHTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+              }
+            }
+          }else if (abs((p1.x()-p2.x()-1.0))<0.0001){
+            if (abs((p2.y()-p1.y()))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tFRONT = FRONT + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACK = BACK + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if(abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTBOTTOM = FRONTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKTOP = BACKTOP + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTTOP = FRONTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKBOTTOM = BACKBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+              }
+            }else if (abs((p2.y()-p1.y()-1.0))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tFRONTLEFT = FRONTLEFT + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKRIGHT = BACKRIGHT + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTLEFTBOTTOM = FRONTLEFTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKRIGHTTOP = BACKRIGHTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTLEFTTOP = FRONTLEFTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKRIGHTBOTTOM = BACKRIGHTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+              }
+            }else if (abs((p1.y()-p2.y()-1.0))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tFRONTRIGHT = FRONTRIGHT + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKLEFT = BACKLEFT + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTRIGHTBOTTOM = FRONTRIGHTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKLEFTTOP = BACKLEFTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tFRONTRIGHTTOP = FRONTRIGHTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBACKLEFTBOTTOM = BACKLEFTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+              }
+            }
+          }else if (abs((p1.x()-p2.x()))<0.0001){
+            if (abs((p2.y()-p1.y()-1.0))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tRIGHT = RIGHT + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tLEFT = LEFT + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tRIGHTTOP = RIGHTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tLEFTBOTTOM = LEFTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tRIGHTBOTTOM = RIGHTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tLEFTTOP = LEFTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+              }
+            }else if (abs((p1.y()-p2.y()-1.0))<0.0001){
+              if (abs((p2.z()-p1.z()))<0.0001){
+                file2 << "NSET\tRIGHT = RIGHT + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tLEFT = LEFT + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tRIGHTBOTTOM = RIGHTBOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tLEFTTOP = LEFTTOP + SURFACE"<<faces[j]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tRIGHTTOP = RIGHTTOP + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tLEFTBOTTOM = LEFTBOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+              }
+            }else if (abs((p1.y()-p2.y()))<0.0001){
+              if (abs((p2.z()-p1.z()-1.0))<0.0001){
+                file2 << "NSET\tTOP = TOP + SURFACE"<<faces[j]->tag()<<"\n";
+                file2 << "NSET\tBOTTOM = BOTTOM + SURFACE"<<faces[i]->tag()<<"\n";
+              }else if (abs((p1.z()-p2.z()-1.0))<0.0001){
+                file2 << "NSET\tTOP = TOP + SURFACE"<<faces[i]->tag()<<"\n";
+                file2 << "NSET\tBOTTOM = BOTTOM + SURFACE"<<faces[j]->tag()<<"\n";
+              }
+            }
+          }
+          count++;
+        }
+      }
+    }
   }
 
   file << "};\n";
-	
+
   printf("\nNumber of exterior face periodicities : %d\n",2*count);
   printf("Total number of exterior faces : %zu\n\n",faces.size());
-	
+
   std::ofstream file3;
   file3.open("MicrostructurePolycrystal3D.geo",std::ios::out | std::ios::app);
-  
+
   std::ofstream file4("MicrostructurePolycrystal3D2.pos");
   file4 << "View \"test\" {\n";
-	
+
   file3 << "Physical Surface(11)={";
-	
+
   count = 0;
   for(it=model->firstFace();it!=model->lastFace();it++)
   {
     gf = *it;
-	if(count>0) file3 << ",";
-	file3 << gf->tag();
-	count++;
+    if(count>0) file3 << ",";
+    file3 << gf->tag();
+    count++;
   }
-	
+
   file3 << "};\n";
-	
+
   for(i=0;i<pairs.size();i++){
     gf1 = pairs[i].first;
-	gf2 = pairs[i].second;
-  
-	edges1 = gf1->edges();
-	edges2 = gf2->edges();  
-	 
-	orientations1 = gf1->edgeOrientations();
-	orientations2 = gf2->edgeOrientations();
-	  
-	indices1.clear();
-	indices2.clear();
-	indices3.clear();	  
-	phase = 1;
-	normal = 0;  
-	  
-	it9 = orientations1.begin(); 
-	it10 = orientations2.begin();
-	for(it8=edges2.begin();it8!=edges2.end();it8++,it10++){
-	  if (*it10==1)
-	  	indices3.push_back((*it8)->tag());
-	  else
-		indices3.push_back(-(*it8)->tag());
-	}
-	int countReverseEdge=0;
-	for(it7=edges1.begin();it7!=edges1.end();it7++,it9++){
-	  v1 = (*it7)->getBeginVertex();
-	  v2 = (*it7)->getEndVertex();
-	  
-	  if (*it9==1)	
-	  	indices1.push_back((*it7)->tag());
-	  else
-		indices1.push_back(-(*it7)->tag());
-		
-	  flag1 = 0;
-	  flag2 = 0;
-	  flag3 = 0;
-	  flag4 = 0;
-			
-	  it10 = orientations2.begin();
-	  for(it8=edges2.begin();it8!=edges2.end();it8++,it10++){
-	        v3 = (*it8)->getBeginVertex();
-		v4 = (*it8)->getEndVertex();
-		  
-		correspondance(fabs(v3->x()-v1->x()),fabs(v3->y()-v1->y()),fabs(v3->z()-v1->z()),e,categories[i],flag1);
-		correspondance(fabs(v4->x()-v2->x()),fabs(v4->y()-v2->y()),fabs(v4->z()-v2->z()),e,categories[i],flag2);
-		  
-		correspondance(fabs(v4->x()-v1->x()),fabs(v4->y()-v1->y()),fabs(v4->z()-v1->z()),e,categories[i],flag3);
-		correspondance(fabs(v3->x()-v2->x()),fabs(v3->y()-v2->y()),fabs(v3->z()-v2->z()),e,categories[i],flag4);
-		  
-		if(flag1 && flag2){
-	          if(phase==1){
-		    mem = it8;
-		    phase = 2;
-		  }
-		  else if(phase==2){
-			mem++;
-			  
-			if(it8==mem){
-			  normal = 1;
-			}
-			else{
-			  normal = -1;
-			}
-			  
-			phase = 3;
-		  }
-		  if (*it9==1)	
-		  	indices2.push_back((*it8)->tag());
-		  else
-			indices2.push_back(-(*it8)->tag());
-		  if (*it9!=*it10)
-		  	countReverseEdge++;
-		  print_segment(SPoint3(v3->x(),v3->y(),v3->z()),SPoint3(v1->x(),v1->y(),v1->z()),file4);
-		  print_segment(SPoint3(v4->x(),v4->y(),v4->z()),SPoint3(v2->x(),v2->y(),v2->z()),file4);
-		}
-		else if(flag3 && flag4){
-		  if(phase==1){
-		    mem = it8;
-			phase = 2;
-		  }
-		  else if(phase==2){
-		    mem++;
-			  
-			if(it8==mem){
-			  normal = 1;
-			}
-			else{
-			  normal = -1;
-			}
-			  
-			phase = 3;
-		  }
-		  if (*it9==1)
-		  	indices2.push_back(-(*it8)->tag());
-		  else
-			indices2.push_back((*it8)->tag());
-		  if (*it9!=*it10)
-                        countReverseEdge++;
-		  print_segment(SPoint3(v4->x(),v4->y(),v4->z()),SPoint3(v1->x(),v1->y(),v1->z()),file4);
-		  print_segment(SPoint3(v3->x(),v3->y(),v3->z()),SPoint3(v2->x(),v2->y(),v2->z()),file4);
-		}
-	  }
-	}
-	  
-	if(indices1.size()!=indices2.size()){
-	  printf("Error\n\n");
-	}
-	  
-	file3 << "Periodic Surface " << gf1->tag() << " {";  
-	  
-	for(j=0;j<indices1.size();j++){
-	  if(j>0) file3 << ",";
-	  file3 << indices1[j];
-	}
-	file3 << "} = " << gf2->tag() << " {";
-	  
-	for(j=0;j<indices2.size();j++){
+    gf2 = pairs[i].second;
+
+    edges1 = gf1->edges();
+    edges2 = gf2->edges();
+
+    orientations1 = gf1->edgeOrientations();
+    orientations2 = gf2->edgeOrientations();
+
+    indices1.clear();
+    indices2.clear();
+    indices3.clear();
+    phase = 1;
+    //normal = 0;
+
+    it9 = orientations1.begin();
+    it10 = orientations2.begin();
+    for(it8=edges2.begin();it8!=edges2.end();it8++,it10++){
+      if (*it10==1)
+        indices3.push_back((*it8)->tag());
+      else
+        indices3.push_back(-(*it8)->tag());
+    }
+    int countReverseEdge=0;
+    for(it7=edges1.begin();it7!=edges1.end();it7++,it9++){
+      v1 = (*it7)->getBeginVertex();
+      v2 = (*it7)->getEndVertex();
+
+      if (*it9==1)
+        indices1.push_back((*it7)->tag());
+      else
+        indices1.push_back(-(*it7)->tag());
+
+      flag1 = 0;
+      flag2 = 0;
+      flag3 = 0;
+      flag4 = 0;
+
+      it10 = orientations2.begin();
+      for(it8=edges2.begin();it8!=edges2.end();it8++,it10++){
+        v3 = (*it8)->getBeginVertex();
+        v4 = (*it8)->getEndVertex();
+
+        correspondance(fabs(v3->x()-v1->x()),fabs(v3->y()-v1->y()),fabs(v3->z()-v1->z()),e,categories[i],flag1);
+        correspondance(fabs(v4->x()-v2->x()),fabs(v4->y()-v2->y()),fabs(v4->z()-v2->z()),e,categories[i],flag2);
+
+        correspondance(fabs(v4->x()-v1->x()),fabs(v4->y()-v1->y()),fabs(v4->z()-v1->z()),e,categories[i],flag3);
+        correspondance(fabs(v3->x()-v2->x()),fabs(v3->y()-v2->y()),fabs(v3->z()-v2->z()),e,categories[i],flag4);
+
+        if(flag1 && flag2){
+          if(phase==1){
+            mem = it8;
+            phase = 2;
+          }
+          else if(phase==2){
+            mem++;
+
+            /*if(it8==mem){
+              normal = 1;
+            }
+            else{
+              normal = -1;
+            }*/
+
+            phase = 3;
+          }
+          if (*it9==1)
+            indices2.push_back((*it8)->tag());
+          else
+            indices2.push_back(-(*it8)->tag());
+          if (*it9!=*it10)
+              countReverseEdge++;
+          print_segment(SPoint3(v3->x(),v3->y(),v3->z()),SPoint3(v1->x(),v1->y(),v1->z()),file4);
+          print_segment(SPoint3(v4->x(),v4->y(),v4->z()),SPoint3(v2->x(),v2->y(),v2->z()),file4);
+        }
+        else if(flag3 && flag4){
+          if(phase==1){
+            mem = it8;
+            phase = 2;
+          }
+          else if(phase==2){
+            mem++;
+
+            /*if(it8==mem){
+              normal = 1;
+            }
+            else{
+              normal = -1;
+            }*/
+
+            phase = 3;
+          }
+          if (*it9==1)
+            indices2.push_back(-(*it8)->tag());
+          else
+            indices2.push_back((*it8)->tag());
+          if (*it9!=*it10)
+            countReverseEdge++;
+          print_segment(SPoint3(v4->x(),v4->y(),v4->z()),SPoint3(v1->x(),v1->y(),v1->z()),file4);
+          print_segment(SPoint3(v3->x(),v3->y(),v3->z()),SPoint3(v2->x(),v2->y(),v2->z()),file4);
+        }
+      }
+    }
+
+    if(indices1.size()!=indices2.size()){
+      printf("Error\n\n");
+    }
+
+    file3 << "Periodic Surface " << gf1->tag() << " {";
+
+    for(j=0;j<indices1.size();j++){
       if(j>0) file3 << ",";
-	  file3 << indices2[j];
-	}
-	  
-	file3 << "};\n";
+      file3 << indices1[j];
+    }
+    file3 << "} = " << gf2->tag() << " {";
+
+    for(j=0;j<indices2.size();j++){
+      if(j>0) file3 << ",";
+      file3 << indices2[j];
+    }
+
+    file3 << "};\n";
   }
-	
+
   file4 << "};\n";
 }
 
 bool voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,double e,int& val){
   bool flag;
-	
+
   flag = 0;
   val = 1000;
-	
+
   if(equal(delta_x,1.0,e) && equal(delta_y,0.0,e) && equal(delta_z,0.0,e)){
     flag = 1;
-	val = 1;
+    val = 1;
   }
   if(equal(delta_x,0.0,e) && equal(delta_y,1.0,e) && equal(delta_z,0.0,e)){
     flag = 1;
-	val = 2;
+    val = 2;
   }
   if(equal(delta_x,0.0,e) && equal(delta_y,0.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
-	val = 3;
+    val = 3;
   }
-	
+
   if(equal(delta_x,1.0,e) && equal(delta_y,1.0,e) && equal(delta_z,0.0,e)){
     flag = 1;
-	val = 4;
+    val = 4;
   }
   if(equal(delta_x,0.0,e) && equal(delta_y,1.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
-	val = 5;
+    val = 5;
   }
   if(equal(delta_x,1.0,e) && equal(delta_y,0.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
-	val = 6;
+    val = 6;
   }
-	
+
   if(equal(delta_x,1.0,e) && equal(delta_y,1.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
-	val = 7;
+    val = 7;
   }
-	
+
   return flag;
 }
 
 void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,double e,int val,bool& flag){
   flag = 0;
-	
+
   if(val==1 && equal(delta_x,1.0,e) && equal(delta_y,0.0,e) && equal(delta_z,0.0,e)){
     flag = 1;
   }
@@ -882,7 +882,7 @@ void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,do
   if(val==3 && equal(delta_x,0.0,e) && equal(delta_y,0.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
   }
-	
+
   if(val==4 && equal(delta_x,1.0,e) && equal(delta_y,1.0,e) && equal(delta_z,0.0,e)){
     flag = 1;
   }
@@ -892,7 +892,7 @@ void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,do
   if(val==6 && equal(delta_x,1.0,e) && equal(delta_y,0.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
   }
-	
+
   if(val==7 && equal(delta_x,1.0,e) && equal(delta_y,1.0,e) && equal(delta_z,1.0,e)){
     flag = 1;
   }
@@ -900,7 +900,7 @@ void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,do
 
 bool voroMetal3D::equal(double x,double y,double e){
   bool flag;
-	
+
   if(x>y-e && x<y+e){
     flag = 1;
   }
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index c1c7a170a178a0fbf68177b11f63bc97f6a5a098..a90921e3b2259e2a738f33448eea46e8060d02bc 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -182,8 +182,8 @@ bool compute4neighbors (GFace *gf,   // the surface
 			SMetric3 &metricField, FILE *crossf = 0) // the mesh metric
 {
 
-  Range<double> rangeU = gf->parBounds(0);
-  Range<double> rangeV = gf->parBounds(1);
+  //Range<double> rangeU = gf->parBounds(0);
+  //Range<double> rangeV = gf->parBounds(1);
 
   // we assume that v is on surface gf
 
@@ -476,7 +476,7 @@ void packingOfParallelograms(GFace* gf,  std::vector<MVertex*> &packed, std::vec
   char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
   FILE *f = fopen(ccc,"w");
   fprintf(f,"View \"\"{\n");
-  for (int i=0;i<vertices.size();i++){
+  for (unsigned int i=0;i<vertices.size();i++){
     //    if(vertices[i]->_v->onWhat() != gf) 
       vertices[i]->print(f,i);
     if(vertices[i]->_v->onWhat() == gf) {
diff --git a/Mesh/yamakawa.cpp b/Mesh/yamakawa.cpp
index 2315d061341399c6deb1dec3a316ba136f460772..4e6f7f7027d607b54202c41a3bb6be56f4370350 100755
--- a/Mesh/yamakawa.cpp
+++ b/Mesh/yamakawa.cpp
@@ -328,9 +328,9 @@ void Recombinator::execute(){
   for(it=model->firstRegion();it!=model->lastRegion();it++)
   {
     gr = *it;
-	if(gr->getNumMeshElements()>0){
-	  execute(gr);
-	}
+    if(gr->getNumMeshElements()>0){
+      execute(gr);
+    }
   }
 }
 
@@ -373,7 +373,7 @@ void Recombinator::init_markings(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	markings.insert(std::pair<MElement*,bool>(element,0));
+    markings.insert(std::pair<MElement*,bool>(element,0));
   }
 }
 
@@ -397,55 +397,55 @@ void Recombinator::patern1(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	max_scaled_jacobian(element,index);
-
-	a = element->getVertex(index);
-	b = element->getVertex((index+1)%4);
-	c = element->getVertex((index+2)%4);
-	d = element->getVertex((index+3)%4);
-
-	already.clear();
-	already.push_back(a);
-	already.push_back(b);
-	already.push_back(c);
-	already.push_back(d);
-	bin1.clear();
-	bin2.clear();
-	bin3.clear();
-	find(b,d,already,bin1);
-	find(b,c,already,bin2);
-	find(c,d,already,bin3);
-
-	for(it1=bin1.begin();it1!=bin1.end();it1++){
-	  p = *it1;
-	  for(it2=bin2.begin();it2!=bin2.end();it2++){
-	    q = *it2;
-	    for(it3=bin3.begin();it3!=bin3.end();it3++){
-		  r = *it3;
-		  if(p!=q && p!=r && q!=r){
-		    already.clear();
-		    already.push_back(a);
-		    already.push_back(b);
-		    already.push_back(c);
-		    already.push_back(d);
-		    already.push_back(p);
-		    already.push_back(q);
-		    already.push_back(r);
-		    bin4.clear();
-		    find(p,q,r,already,bin4);
-		    for(it4=bin4.begin();it4!=bin4.end();it4++){
-			  s = *it4;
-			  hex = Hex(a,b,q,c,d,p,s,r);
-			  quality = min_scaled_jacobian(hex);
-			  hex.set_quality(quality);
-			  if(valid(hex)){
-			    potential.push_back(hex);
-			  }
-		    }
-		  }
-		}
-	  }
-	}
+    max_scaled_jacobian(element,index);
+
+    a = element->getVertex(index);
+    b = element->getVertex((index+1)%4);
+    c = element->getVertex((index+2)%4);
+    d = element->getVertex((index+3)%4);
+
+    already.clear();
+    already.push_back(a);
+    already.push_back(b);
+    already.push_back(c);
+    already.push_back(d);
+    bin1.clear();
+    bin2.clear();
+    bin3.clear();
+    find(b,d,already,bin1);
+    find(b,c,already,bin2);
+    find(c,d,already,bin3);
+
+    for(it1=bin1.begin();it1!=bin1.end();it1++){
+      p = *it1;
+      for(it2=bin2.begin();it2!=bin2.end();it2++){
+        q = *it2;
+        for(it3=bin3.begin();it3!=bin3.end();it3++){
+          r = *it3;
+          if(p!=q && p!=r && q!=r){
+            already.clear();
+            already.push_back(a);
+            already.push_back(b);
+            already.push_back(c);
+            already.push_back(d);
+            already.push_back(p);
+            already.push_back(q);
+            already.push_back(r);
+            bin4.clear();
+            find(p,q,r,already,bin4);
+            for(it4=bin4.begin();it4!=bin4.end();it4++){
+              s = *it4;
+              hex = Hex(a,b,q,c,d,p,s,r);
+              quality = min_scaled_jacobian(hex);
+              hex.set_quality(quality);
+              if(valid(hex)){
+                potential.push_back(hex);
+              }
+            }
+          }
+        }
+      }
+    }
   }
 }
 
@@ -464,36 +464,36 @@ void Recombinator::patern2(GRegion* gr){
     diagonal(element,index1,index2);
     two_others(index1,index2,index3,index4);
 
-	b = element->getVertex(index1);
+    b = element->getVertex(index1);
     d = element->getVertex(index2);
     a = element->getVertex(index3);
     c = element->getVertex(index4);
 
-	verif.clear();
-	find(b,d,verif);
-	if(verif.size()==6){
-	  s = find(a,b,d,c,verif);
-	  p = find(b,c,d,a,verif);
-	  if(s!=0 && p!=0){
-	    r = find(s,b,d,a,verif);
-	    q = find(p,b,d,c,verif);
-		if(r!=0 && q!=0){
-	      hex = Hex(a,s,b,c,d,r,q,p);
-	      quality = min_scaled_jacobian(hex);
-	      hex.set_quality(quality);
-	      if(valid(hex)){
-	        potential.push_back(hex);
-	      }
-
-	      hex = Hex(a,c,d,s,b,p,q,r);
-	      quality = min_scaled_jacobian(hex);
-	      hex.set_quality(quality);
-	      if(valid(hex)){
-	        potential.push_back(hex);
-	      }
-		}
-	  }
-	}
+    verif.clear();
+    find(b,d,verif);
+    if(verif.size()==6){
+      s = find(a,b,d,c,verif);
+      p = find(b,c,d,a,verif);
+      if(s!=0 && p!=0){
+        r = find(s,b,d,a,verif);
+        q = find(p,b,d,c,verif);
+        if(r!=0 && q!=0){
+          hex = Hex(a,s,b,c,d,r,q,p);
+          quality = min_scaled_jacobian(hex);
+          hex.set_quality(quality);
+          if(valid(hex)){
+            potential.push_back(hex);
+          }
+
+          hex = Hex(a,c,d,s,b,p,q,r);
+          quality = min_scaled_jacobian(hex);
+          hex.set_quality(quality);
+          if(valid(hex)){
+            potential.push_back(hex);
+          }
+        }
+      }
+    }
   }
 }
 
@@ -513,108 +513,108 @@ void Recombinator::patern3(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	diagonal(element,index1,index2);
-	two_others(index1,index2,index3,index4);
-
-	b = element->getVertex(index1);
-	d = element->getVertex(index2);
-	a = element->getVertex(index3);
-	c = element->getVertex(index4);
-
-	verif1.clear();
-	verif2.clear();
-	find(b,d,verif1);
-	find(a,c,verif2);
-
-	if(verif1.size()==4 && verif2.size()==4){
-	  fA = find(b,d,a,c,verif1);
-	  fB = find(b,d,c,a,verif1);
-	  bA = find(a,c,b,d,verif2);
-	  bB = find(a,c,d,b,verif2);
-
-	  if(fA!=0 && fB!=0 && bA!=0 && bB!=0 && fA!=fB && bA!=bB){
-		if(scalar(fA,fB,a,b)>scalar(fA,fB,b,c) && scalar(bA,bB,a,b)>scalar(bA,bB,b,c)){
-		  if(distance(fA,b,c)<distance(fB,b,c)){
-		    p = fA;
-			q = fB;
-		  }
-		  else{
-		    p = fB;
-			q = fA;
-		  }
-
-		  if(distance(bA,b,c)<distance(bB,b,c)){
-		    r = bA;
-			s = bB;
-		  }
-		  else{
-		    r = bB;
-			s = bA;
-		  }
-
-		  c1 = linked(b,p);
-		  c2 = linked(c,p);
-		  c3 = linked(p,q);
-		  c4 = linked(a,q);
-		  c5 = linked(d,q);
-
-		  c6 = linked(b,r);
-		  c7 = linked(c,r);
-		  c8 = linked(r,s);
-		  c9 = linked(a,s);
-		  c10 = linked(d,s);
-
-		  if(c1 && c2 && c3 && c4 && c5 && c6 && c7 && c8 && c9 && c10){
-		    hex = Hex(p,c,r,b,q,d,s,a);
-			quality = min_scaled_jacobian(hex);
-			hex.set_quality(quality);
-			if(valid(hex)){
-			  potential.push_back(hex);
-			}
-		  }
-		}
-		else if(scalar(fA,fB,a,b)<=scalar(fA,fB,b,c) && scalar(bA,bB,a,b)<=scalar(bA,bB,b,c)){
-		  if(distance(fA,a,b)<distance(fB,a,b)){
-		    p = fA;
-			q = fB;
-		  }
-		  else{
-		    p = fB;
-			q = fA;
-		  }
-
-		  if(distance(bA,a,b)<distance(bB,a,b)){
-		    r = bA;
-			s = bB;
-		  }
-		  else{
-		    r = bB;
-			s = bA;
-		  }
-
-		  c1 = linked(b,p);
-		  c2 = linked(a,p);
-		  c3 = linked(p,q);
-		  c4 = linked(c,q);
-		  c5 = linked(d,q);
-
-		  c6 = linked(b,r);
-		  c7 = linked(a,r);
-		  c8 = linked(r,s);
-		  c9 = linked(c,s);
-		  c10 = linked(d,s);
-
-		  if(c1 && c2 && c3 && c4 && c5 && c6 && c7 && c8 && c9 && c10){
-		    hex = Hex(p,b,r,a,q,c,s,d);
-			quality = min_scaled_jacobian(hex);
-			hex.set_quality(quality);
-			if(valid(hex)){
-		      potential.push_back(hex);
-			}
-		  }
-		}
-	  }
-	}
+    diagonal(element,index1,index2);
+    two_others(index1,index2,index3,index4);
+
+    b = element->getVertex(index1);
+    d = element->getVertex(index2);
+    a = element->getVertex(index3);
+    c = element->getVertex(index4);
+
+    verif1.clear();
+    verif2.clear();
+    find(b,d,verif1);
+    find(a,c,verif2);
+
+    if(verif1.size()==4 && verif2.size()==4){
+      fA = find(b,d,a,c,verif1);
+      fB = find(b,d,c,a,verif1);
+      bA = find(a,c,b,d,verif2);
+      bB = find(a,c,d,b,verif2);
+
+      if(fA!=0 && fB!=0 && bA!=0 && bB!=0 && fA!=fB && bA!=bB){
+        if(scalar(fA,fB,a,b)>scalar(fA,fB,b,c) && scalar(bA,bB,a,b)>scalar(bA,bB,b,c)){
+          if(distance(fA,b,c)<distance(fB,b,c)){
+            p = fA;
+            q = fB;
+          }
+          else{
+            p = fB;
+            q = fA;
+          }
+
+          if(distance(bA,b,c)<distance(bB,b,c)){
+            r = bA;
+            s = bB;
+          }
+          else{
+            r = bB;
+            s = bA;
+          }
+
+          c1 = linked(b,p);
+          c2 = linked(c,p);
+          c3 = linked(p,q);
+          c4 = linked(a,q);
+          c5 = linked(d,q);
+
+          c6 = linked(b,r);
+          c7 = linked(c,r);
+          c8 = linked(r,s);
+          c9 = linked(a,s);
+          c10 = linked(d,s);
+
+          if(c1 && c2 && c3 && c4 && c5 && c6 && c7 && c8 && c9 && c10){
+            hex = Hex(p,c,r,b,q,d,s,a);
+            quality = min_scaled_jacobian(hex);
+            hex.set_quality(quality);
+            if(valid(hex)){
+              potential.push_back(hex);
+            }
+          }
+        }
+        else if(scalar(fA,fB,a,b)<=scalar(fA,fB,b,c) && scalar(bA,bB,a,b)<=scalar(bA,bB,b,c)){
+          if(distance(fA,a,b)<distance(fB,a,b)){
+            p = fA;
+            q = fB;
+          }
+          else{
+            p = fB;
+            q = fA;
+          }
+
+          if(distance(bA,a,b)<distance(bB,a,b)){
+            r = bA;
+            s = bB;
+          }
+          else{
+            r = bB;
+            s = bA;
+          }
+
+          c1 = linked(b,p);
+          c2 = linked(a,p);
+          c3 = linked(p,q);
+          c4 = linked(c,q);
+          c5 = linked(d,q);
+
+          c6 = linked(b,r);
+          c7 = linked(a,r);
+          c8 = linked(r,s);
+          c9 = linked(c,s);
+          c10 = linked(d,s);
+
+          if(c1 && c2 && c3 && c4 && c5 && c6 && c7 && c8 && c9 && c10){
+            hex = Hex(p,b,r,a,q,c,s,d);
+            quality = min_scaled_jacobian(hex);
+            hex.set_quality(quality);
+            if(valid(hex)){
+              potential.push_back(hex);
+            }
+          }
+        }
+      }
+    }
   }
 }
 
@@ -639,73 +639,73 @@ void Recombinator::merge(GRegion* gr){
   for(i=0;i<potential.size();i++){
     hex = potential[i];
 
-	threshold = 0.25;
-	if(hex.get_quality()<threshold){
-	  break;
-	}
-
-	a = hex.get_a();
-	b = hex.get_b();
-	c = hex.get_c();
-	d = hex.get_d();
-	e = hex.get_e();
-	f = hex.get_f();
-	g = hex.get_g();
-	h = hex.get_h();
-
-	parts.clear();
-	find(a,hex,parts);
-	find(b,hex,parts);
-	find(c,hex,parts);
-	find(d,hex,parts);
-	find(e,hex,parts);
-	find(f,hex,parts);
-	find(g,hex,parts);
-	find(h,hex,parts);
-
-	flag = 1;
-	for(it=parts.begin();it!=parts.end();it++){
-	  element = *it;
-	  it2 = markings.find(element);
-	  if(it2->second==1 && !sliver(element,hex)){
-	    flag = 0;
-		break;
-	  }
-	}
-	if(!flag) continue;
-
-	if(!valid(hex,parts)){
-	  continue;
-	}
-
-	if(!conformityA(hex)){
-	  continue;
-	}
-
-	if(!conformityB(hex)){
-	  continue;
-	}
-
-	if(!conformityC(hex)){
-	  continue;
-	}
-
-	if(!faces_statuquo(hex)){
-	  continue;
-	}
-
-	//printf("%d - %d/%d - %f\n",count,i,(int)potential.size(),hex.get_quality());
-	quality = quality + hex.get_quality();
-	for(it=parts.begin();it!=parts.end();it++){
-	  element = *it;
-	  it2 = markings.find(element);
-	  it2->second = 1;
-	}
-	gr->addHexahedron(new MHexahedron(a,b,c,d,e,f,g,h));
-	build_hash_tableA(hex);
-	build_hash_tableB(hex);
-	build_hash_tableC(hex);
-	count++;
+    threshold = 0.25;
+    if(hex.get_quality()<threshold){
+      break;
+    }
+
+    a = hex.get_a();
+    b = hex.get_b();
+    c = hex.get_c();
+    d = hex.get_d();
+    e = hex.get_e();
+    f = hex.get_f();
+    g = hex.get_g();
+    h = hex.get_h();
+
+    parts.clear();
+    find(a,hex,parts);
+    find(b,hex,parts);
+    find(c,hex,parts);
+    find(d,hex,parts);
+    find(e,hex,parts);
+    find(f,hex,parts);
+    find(g,hex,parts);
+    find(h,hex,parts);
+
+    flag = 1;
+    for(it=parts.begin();it!=parts.end();it++){
+      element = *it;
+      it2 = markings.find(element);
+      if(it2->second==1 && !sliver(element,hex)){
+        flag = 0;
+        break;
+      }
+    }
+    if(!flag) continue;
+
+    if(!valid(hex,parts)){
+      continue;
+    }
+
+    if(!conformityA(hex)){
+      continue;
+    }
+
+    if(!conformityB(hex)){
+      continue;
+    }
+
+    if(!conformityC(hex)){
+      continue;
+    }
+
+    if(!faces_statuquo(hex)){
+      continue;
+    }
+
+    //printf("%d - %d/%d - %f\n",count,i,(int)potential.size(),hex.get_quality());
+    quality = quality + hex.get_quality();
+    for(it=parts.begin();it!=parts.end();it++){
+      element = *it;
+      it2 = markings.find(element);
+      it2->second = 1;
+    }
+    gr->addHexahedron(new MHexahedron(a,b,c,d,e,f,g,h));
+    build_hash_tableA(hex);
+    build_hash_tableB(hex);
+    build_hash_tableC(hex);
+    count++;
   }
 
   opt.clear();
@@ -715,10 +715,10 @@ void Recombinator::merge(GRegion* gr){
 
   for(i=0;i<opt.size();i++){
     element = (MElement*)(opt[i]);
-	it2 = markings.find(element);
-	if(it2->second==0){
-	  gr->tetrahedra.push_back(opt[i]);
-	}
+    it2 = markings.find(element);
+    if(it2->second==0){
+      gr->tetrahedra.push_back(opt[i]);
+    }
   }
 
   printf("hexahedra average quality (0->1) : %f\n",quality/count);
@@ -746,69 +746,69 @@ void Recombinator::improved_merge(GRegion* gr){
   for(i=0;i<potential.size();i++){
     hex = potential[i];
 
-	threshold = 0.25;
-	if(hex.get_quality()<threshold){
-	  break;
-	}
-
-	a = hex.get_a();
-	b = hex.get_b();
-	c = hex.get_c();
-	d = hex.get_d();
-	e = hex.get_e();
-	f = hex.get_f();
-	g = hex.get_g();
-	h = hex.get_h();
-
-	parts.clear();
-	find(a,hex,parts);
-	find(b,hex,parts);
-	find(c,hex,parts);
-	find(d,hex,parts);
-	find(e,hex,parts);
-	find(f,hex,parts);
-	find(g,hex,parts);
-	find(h,hex,parts);
-
-	flag = 1;
-	for(it=parts.begin();it!=parts.end();it++){
-	  element = *it;
-	  it2 = markings.find(element);
-	  if(it2->second==1 && !sliver(element,hex)){
-	    flag = 0;
-		break;
-	  }
-	}
-	if(!flag) continue;
-
-	if(!valid(hex,parts)){
-	  continue;
-	}
-
-	if(!conformityA(hex)){
-	  continue;
-	}
-
-	if(!conformityB(hex)){
-	  continue;
-	}
-
-	if(!conformityC(hex)){
-	  continue;
-	}
-
-	//printf("%d - %d/%d - %f\n",count,i,(int)potential.size(),hex.get_quality());
-	quality = quality + hex.get_quality();
-	for(it=parts.begin();it!=parts.end();it++){
-	  element = *it;
-	  it2 = markings.find(element);
-	  it2->second = 1;
-	}
-	gr->addHexahedron(new MHexahedron(a,b,c,d,e,f,g,h));
-	build_hash_tableA(hex);
-	build_hash_tableB(hex);
-	build_hash_tableC(hex);
-	count++;
+    threshold = 0.25;
+    if(hex.get_quality()<threshold){
+      break;
+    }
+
+    a = hex.get_a();
+    b = hex.get_b();
+    c = hex.get_c();
+    d = hex.get_d();
+    e = hex.get_e();
+    f = hex.get_f();
+    g = hex.get_g();
+    h = hex.get_h();
+
+    parts.clear();
+    find(a,hex,parts);
+    find(b,hex,parts);
+    find(c,hex,parts);
+    find(d,hex,parts);
+    find(e,hex,parts);
+    find(f,hex,parts);
+    find(g,hex,parts);
+    find(h,hex,parts);
+
+    flag = 1;
+    for(it=parts.begin();it!=parts.end();it++){
+      element = *it;
+      it2 = markings.find(element);
+      if(it2->second==1 && !sliver(element,hex)){
+        flag = 0;
+        break;
+      }
+    }
+    if(!flag) continue;
+
+    if(!valid(hex,parts)){
+      continue;
+    }
+
+    if(!conformityA(hex)){
+      continue;
+    }
+
+    if(!conformityB(hex)){
+      continue;
+    }
+
+    if(!conformityC(hex)){
+      continue;
+    }
+
+    //printf("%d - %d/%d - %f\n",count,i,(int)potential.size(),hex.get_quality());
+    quality = quality + hex.get_quality();
+    for(it=parts.begin();it!=parts.end();it++){
+      element = *it;
+      it2 = markings.find(element);
+      it2->second = 1;
+    }
+    gr->addHexahedron(new MHexahedron(a,b,c,d,e,f,g,h));
+    build_hash_tableA(hex);
+    build_hash_tableB(hex);
+    build_hash_tableC(hex);
+    count++;
   }
 
   opt.clear();
@@ -818,10 +818,10 @@ void Recombinator::improved_merge(GRegion* gr){
 
   for(i=0;i<opt.size();i++){
     element = (MElement*)(opt[i]);
-	it2 = markings.find(element);
-	if(it2->second==0){
-	  gr->tetrahedra.push_back(opt[i]);
-	}
+    it2 = markings.find(element);
+    if(it2->second==0){
+      gr->tetrahedra.push_back(opt[i]);
+    }
   }
 
   printf("hexahedra average quality (0->1) : %f\n",quality/count);
@@ -833,7 +833,7 @@ void Recombinator::rearrange(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	element->setVolumePositive();
+    element->setVolumePositive();
   }
 }
 
@@ -850,15 +850,15 @@ void Recombinator::statistics(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	volume = element->getVolume();
+    volume = element->getVolume();
 
-	if(element->getNumVertices()==8){
-	  hex_nbr = hex_nbr + 1;
-	  hex_volume = hex_volume + volume;
-	}
+    if(element->getNumVertices()==8){
+      hex_nbr = hex_nbr + 1;
+      hex_volume = hex_volume + volume;
+    }
 
-	all_nbr = all_nbr + 1;
-	all_volume = all_volume + volume;
+    all_nbr = all_nbr + 1;
+    all_volume = all_volume + volume;
   }
 
   printf("percentage of hexahedra (number) : %.2f\n",hex_nbr*100.0/all_nbr);
@@ -883,17 +883,17 @@ void Recombinator::build_tuples(GRegion* gr){
   {
     gf = *it;
 
-	for(i=0;i<gf->getNumMeshElements();i++){
-	  element = gf->getMeshElement(i);
+    for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);
 
-	  if(element->getNumVertices()==3){
-		a = element->getVertex(0);
-		b = element->getVertex(1);
-		c = element->getVertex(2);
+      if(element->getNumVertices()==3){
+        a = element->getVertex(0);
+        b = element->getVertex(1);
+        c = element->getVertex(2);
 
-		tuples.insert(Tuple(a,b,c,element,gf));
-	  }
-	}
+        tuples.insert(Tuple(a,b,c,element,gf));
+      }
+    }
   }
 }
 
@@ -911,23 +911,23 @@ void Recombinator::modify_surfaces(GRegion* gr){
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
 
-	if(element->getNumVertices()==8){
-	  a = element->getVertex(0);
-	  b = element->getVertex(1);
-	  c = element->getVertex(2);
-	  d = element->getVertex(3);
-	  e = element->getVertex(4);
-	  f = element->getVertex(5);
-	  g = element->getVertex(6);
-	  h = element->getVertex(7);
-
-	  modify_surfaces(a,b,c,d);
-	  modify_surfaces(e,f,g,h);
-	  modify_surfaces(a,e,h,d);
-	  modify_surfaces(b,f,g,c);
-	  modify_surfaces(a,e,f,b);
-	  modify_surfaces(d,h,g,c);
-	}
+    if(element->getNumVertices()==8){
+      a = element->getVertex(0);
+      b = element->getVertex(1);
+      c = element->getVertex(2);
+      d = element->getVertex(3);
+      e = element->getVertex(4);
+      f = element->getVertex(5);
+      g = element->getVertex(6);
+      h = element->getVertex(7);
+
+      modify_surfaces(a,b,c,d);
+      modify_surfaces(e,f,g,h);
+      modify_surfaces(a,e,h,d);
+      modify_surfaces(b,f,g,c);
+      modify_surfaces(a,e,f,b);
+      modify_surfaces(d,h,g,c);
+    }
   }
 
   faces = gr->faces();
@@ -936,37 +936,37 @@ void Recombinator::modify_surfaces(GRegion* gr){
   {
     gf = *it;
 
-	opt.clear();
+    opt.clear();
 
-	for(i=0;i<gf->getNumMeshElements();i++){
-	  element = gf->getMeshElement(i);
+    for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);
 
-	  if(element->getNumVertices()==3){
-	    it2 = triangles.find(element);
-		if(it2==triangles.end()){
-		  opt.push_back(element);
-		}
-	  }
-	}
+      if(element->getNumVertices()==3){
+        it2 = triangles.find(element);
+        if(it2==triangles.end()){
+          opt.push_back(element);
+        }
+      }
+    }
 
-	gf->triangles.clear();
+    gf->triangles.clear();
 
-	for(i=0;i<opt.size();i++){
-	  gf->triangles.push_back((MTriangle*)opt[i]);
-	}
+    for(i=0;i<opt.size();i++){
+      gf->triangles.push_back((MTriangle*)opt[i]);
+    }
   }
 }
 
 void Recombinator::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
   bool flag1,flag2;
   MElement *element1,*element2;
-  GFace *gf1,*gf2;
+  GFace *gf1;//,*gf2;
   Tuple tuple1,tuple2;
   std::multiset<Tuple>::iterator it1;
   std::multiset<Tuple>::iterator it2;
 
   gf1 = NULL;
-  gf2 = NULL;
+  //gf2 = NULL;
 
   tuple1 = Tuple(a,b,c);
   tuple2 = Tuple(c,d,a);
@@ -979,37 +979,37 @@ void Recombinator::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  element1 = it1->get_element();
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      element1 = it1->get_element();
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  element2 = it2->get_element();
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      element2 = it2->get_element();
+      //gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     triangles.insert(element1);
-	triangles.insert(element2);
+    triangles.insert(element2);
 
-	gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
+    gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
   }
 
   tuple1 = Tuple(a,b,d);
@@ -1024,36 +1024,36 @@ void Recombinator::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
       break;
-	}
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  element1 = it1->get_element();
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      element1 = it1->get_element();
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  element2 = it2->get_element();
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      element2 = it2->get_element();
+      //gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     triangles.insert(element1);
-	triangles.insert(element2);
+    triangles.insert(element2);
 
-	gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
+    gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
   }
 }
 
@@ -1127,33 +1127,33 @@ double Recombinator::diagonal(MElement* element,int& index1,int& index2){
 
   if(l1>=l2 && l1>=l3 && l1>=l4 && l1>=l5 && l1>=l6){
     index1 = 0;
-	index2 = 1;
-	max = l1;
+    index2 = 1;
+    max = l1;
   }
   else if(l2>=l1 && l2>=l3 && l2>=l4 && l2>=l5 && l2>=l6){
     index1 = 0;
-	index2 = 2;
-	max = l2;
+    index2 = 2;
+    max = l2;
   }
   else if(l3>=l1 && l3>=l2 && l3>=l4 && l3>=l5 && l3>=l6){
     index1 = 0;
-	index2 = 3;
-	max = l3;
+    index2 = 3;
+    max = l3;
   }
   else if(l4>=l1 && l4>=l2 && l4>=l3 && l4>=l5 && l4>=l6){
     index1 = 1;
-	index2 = 2;
-	max = l4;
+    index2 = 2;
+    max = l4;
   }
   else if(l5>=l1 && l5>=l2 && l5>=l3 && l5>=l4 && l5>=l6){
     index1 = 2;
-	index2 = 3;
-	max = l5;
+    index2 = 3;
+    max = l5;
   }
   else if(l6>=l1 && l6>=l2 && l6>=l3 && l6>=l4 && l6>=l5){
     index1 = 3;
-	index2 = 1;
-	max = l6;
+    index2 = 1;
+    max = l6;
   }
 
   return max;
@@ -1203,16 +1203,16 @@ void Recombinator::two_others(int index1,int index2,int& index3,int& index4){
 
   for(i=0;i<4;i++){
     if(i!=index1 && i!=index2){
-	  index3 = i;
-	  break;
-	}
+      index3 = i;
+      break;
+    }
   }
 
   for(i=0;i<4;i++){
     if(i!=index1 && i!=index2 && i!=index3){
-	  index4 = i;
-	  break;
-	}
+      index4 = i;
+      break;
+    }
   }
 }
 
@@ -1333,10 +1333,10 @@ bool Recombinator::linked(MVertex* v1,MVertex* v2){
   flag = 0;
 
   for(it2=(it->second).begin();it2!=(it->second).end();it2++){
-	if(*it2==v2){
-	  flag = 1;
-	  break;
-	}
+    if(*it2==v2){
+      flag = 1;
+      break;
+    }
   }
 
   return flag;
@@ -1384,18 +1384,18 @@ void Recombinator::find(MVertex* vertex,Hex hex,std::set<MElement*>& final){
 
   for(it2=(it->second).begin();it2!=(it->second).end();it2++){
     a = (*it2)->getVertex(0);
-	b = (*it2)->getVertex(1);
-	c = (*it2)->getVertex(2);
-	d = (*it2)->getVertex(3);
+    b = (*it2)->getVertex(1);
+    c = (*it2)->getVertex(2);
+    d = (*it2)->getVertex(3);
 
-	flag1 = inclusion(a,hex);
-	flag2 = inclusion(b,hex);
-	flag3 = inclusion(c,hex);
-	flag4 = inclusion(d,hex);
+    flag1 = inclusion(a,hex);
+    flag2 = inclusion(b,hex);
+    flag3 = inclusion(c,hex);
+    flag4 = inclusion(d,hex);
 
-	if(flag1 && flag2 && flag3 && flag4){
-	  final.insert(*it2);
-	}
+    if(flag1 && flag2 && flag3 && flag4){
+      final.insert(*it2);
+    }
   }
 }
 
@@ -1411,38 +1411,38 @@ MVertex* Recombinator::find(MVertex* v1,MVertex* v2,MVertex* v3,MVertex* already
   for(it=bin.begin();it!=bin.end();it++){
     element = *it;
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
-
-	flag1 = inclusion(v1,a,b,c,d);
-	flag2 = inclusion(v2,a,b,c,d);
-	flag3 = inclusion(v3,a,b,c,d);
-	flag4 = inclusion(already,a,b,c,d);
-
-	if(flag1 && flag2 && flag3 && !flag4){
-	  if(a!=v1 && a!=v2 && a!=v3){
-	    pointer = a;
-	  }
-	  else if(b!=v1 && b!=v2 && b!=v3){
-	    pointer = b;
-	  }
-	  else if(c!=v1 && c!=v2 && c!=v3){
-	    pointer = c;
-	  }
-	  else{
-	    pointer = d;
-	  }
-	  break;
-	}
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
+
+    flag1 = inclusion(v1,a,b,c,d);
+    flag2 = inclusion(v2,a,b,c,d);
+    flag3 = inclusion(v3,a,b,c,d);
+    flag4 = inclusion(already,a,b,c,d);
+
+    if(flag1 && flag2 && flag3 && !flag4){
+      if(a!=v1 && a!=v2 && a!=v3){
+        pointer = a;
+      }
+      else if(b!=v1 && b!=v2 && b!=v3){
+        pointer = b;
+      }
+      else if(c!=v1 && c!=v2 && c!=v3){
+        pointer = c;
+      }
+      else{
+        pointer = d;
+      }
+      break;
+    }
   }
 
   return pointer;
 }
 
 void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MVertex*>& bin2,
-								const std::vector<MVertex*>& already,std::set<MVertex*>& final){
+                                const std::vector<MVertex*>& already,std::set<MVertex*>& final){
   size_t i;
   bool ok;
   std::set<MVertex*> temp;
@@ -1451,23 +1451,23 @@ void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MV
   std::set_intersection(bin1.begin(),bin1.end(),bin2.begin(),bin2.end(),std::inserter(temp,temp.end()));
 
   for(it=temp.begin();it!=temp.end();it++){
-	ok = 1;
+    ok = 1;
 
-	for(i=0;i<already.size();i++){
-	  if((*it)==already[i]){
-	    ok = 0;
-		break;
-	  }
-	}
+    for(i=0;i<already.size();i++){
+      if((*it)==already[i]){
+        ok = 0;
+        break;
+      }
+    }
 
-	if(ok){
-	  final.insert(*it);
-	}
+    if(ok){
+      final.insert(*it);
+    }
   }
 }
 
 void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MVertex*>& bin2,const std::set<MVertex*>& bin3,
-								const std::vector<MVertex*>& already,std::set<MVertex*>& final){
+                                const std::vector<MVertex*>& already,std::set<MVertex*>& final){
   size_t i;
   bool ok;
   std::set<MVertex*> temp;
@@ -1480,16 +1480,16 @@ void Recombinator::intersection(const std::set<MVertex*>& bin1,const std::set<MV
   for(it=temp2.begin();it!=temp2.end();it++){
     ok = 1;
 
-	for(i=0;i<already.size();i++){
-	  if((*it)==already[i]){
-	    ok = 0;
-		break;
-	  }
+    for(i=0;i<already.size();i++){
+      if((*it)==already[i]){
+        ok = 0;
+        break;
+      }
     }
 
-	if(ok){
-	  final.insert(*it);
-	}
+    if(ok){
+      final.insert(*it);
+    }
   }
 }
 
@@ -1539,19 +1539,19 @@ bool Recombinator::inclusion(MVertex* v1,MVertex* v2,MVertex* v3,const std::set<
   for(it=bin.begin();it!=bin.end();it++){
     element = *it;
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
 
-	flag1 = inclusion(v1,a,b,c,d);
-	flag2 = inclusion(v2,a,b,c,d);
-	flag3 = inclusion(v3,a,b,c,d);
+    flag1 = inclusion(v1,a,b,c,d);
+    flag2 = inclusion(v2,a,b,c,d);
+    flag3 = inclusion(v3,a,b,c,d);
 
-	if(flag1 && flag2 && flag3){
-	  ok = 1;
-	  break;
-	}
+    if(flag1 && flag2 && flag3){
+      ok = 1;
+      break;
+    }
   }
 
   return ok;
@@ -1566,15 +1566,15 @@ bool Recombinator::inclusion(Facet facet){
 
   while(it!=hash_tableA.end()){
     if(facet.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(facet.same_vertices(*it)){
-	  flag = 1;
-	  break;
-	}
+    if(facet.same_vertices(*it)){
+      flag = 1;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   return flag;
@@ -1589,15 +1589,15 @@ bool Recombinator::inclusion(Diagonal diagonal){
 
   while(it!=hash_tableB.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 1;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 1;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   return flag;
@@ -1612,15 +1612,15 @@ bool Recombinator::duplicate(Diagonal diagonal){
 
   while(it!=hash_tableC.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 1;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 1;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   return flag;
@@ -1802,34 +1802,34 @@ bool Recombinator::faces_statuquo(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     if(gf1!=gf2){
-	  ok = 0;
-	}
+      ok = 0;
+    }
   }
 
   tuple1 = Tuple(a,b,d);
@@ -1843,34 +1843,34 @@ bool Recombinator::faces_statuquo(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     if(gf1!=gf2){
-	  ok = 0;
-	}
+      ok = 0;
+    }
   }
 
   return ok;
@@ -1896,16 +1896,16 @@ void Recombinator::build_vertex_to_vertices(GRegion* gr){
 
       it = vertex_to_vertices.find(a);
       if(it!=vertex_to_vertices.end()){
-	it->second.insert(b);
-	it->second.insert(c);
-	it->second.insert(d);
+        it->second.insert(b);
+        it->second.insert(c);
+        it->second.insert(d);
       }
       else{
-	bin.clear();
-	bin.insert(b);
-	bin.insert(c);
-	bin.insert(d);
-	vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(a,bin));
+        bin.clear();
+        bin.insert(b);
+        bin.insert(c);
+        bin.insert(d);
+        vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(a,bin));
       }
     }
   }
@@ -1923,19 +1923,19 @@ void Recombinator::build_vertex_to_elements(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	for(j=0;j<element->getNumVertices();j++){
-	  vertex = element->getVertex(j);
+    for(j=0;j<element->getNumVertices();j++){
+      vertex = element->getVertex(j);
 
-	  it = vertex_to_elements.find(vertex);
-	  if(it!=vertex_to_elements.end()){
-	    it->second.insert(element);
-	  }
-	  else{
-	    bin.clear();
-		bin.insert(element);
-		vertex_to_elements.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
-	  }
-	}
+      it = vertex_to_elements.find(vertex);
+      if(it!=vertex_to_elements.end()){
+        it->second.insert(element);
+      }
+      else{
+        bin.clear();
+        bin.insert(element);
+        vertex_to_elements.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
+      }
+    }
   }
 }
 
@@ -1976,15 +1976,15 @@ void Recombinator::build_hash_tableA(Facet facet){
 
   while(it!=hash_tableA.end()){
     if(facet.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(facet.same_vertices(*it)){
-	  flag = 0;
-	  break;
-	}
+    if(facet.same_vertices(*it)){
+      flag = 0;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   if(flag){
@@ -2027,15 +2027,15 @@ void Recombinator::build_hash_tableB(Diagonal diagonal){
 
   while(it!=hash_tableB.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 0;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 0;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   if(flag){
@@ -2079,15 +2079,15 @@ void Recombinator::build_hash_tableC(Diagonal diagonal){
 
   while(it!=hash_tableC.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 0;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 0;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   if(flag){
@@ -2109,14 +2109,14 @@ void Recombinator::print_vertex_to_vertices(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	for(j=0;j<element->getNumVertices();j++){
-	  vertex = element->getVertex(j);
-	  p1 = SPoint3(vertex->x(),vertex->y(),vertex->z());
-	  it = vertex_to_vertices.find(vertex);
-	  for(it2=(it->second).begin();it2!=(it->second).end();it2++){
-	    p2 = SPoint3((*it2)->x(),(*it2)->y(),(*it2)->z());
-		print_segment(p1,p2,file);
-	  }
+    for(j=0;j<element->getNumVertices();j++){
+      vertex = element->getVertex(j);
+      p1 = SPoint3(vertex->x(),vertex->y(),vertex->z());
+      it = vertex_to_vertices.find(vertex);
+      for(it2=(it->second).begin();it2!=(it->second).end();it2++){
+        p2 = SPoint3((*it2)->x(),(*it2)->y(),(*it2)->z());
+        print_segment(p1,p2,file);
+      }
     }
   }
   file << "};\n";
@@ -2132,12 +2132,12 @@ void Recombinator::print_vertex_to_elements(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	for(j=0;j<element->getNumVertices();j++){
-	  vertex = element->getVertex(j);
-	  it = vertex_to_elements.find(vertex);
-	  it2 = vertex_to_vertices.find(vertex);
-	  printf("%d %d\n",(int)(it->second).size(),(int)(it2->second).size());
-	}
+    for(j=0;j<element->getNumVertices();j++){
+      vertex = element->getVertex(j);
+      it = vertex_to_elements.find(vertex);
+      it2 = vertex_to_vertices.find(vertex);
+      printf("%d %d\n",(int)(it->second).size(),(int)(it2->second).size());
+    }
   }
 }
 
@@ -2190,19 +2190,19 @@ double Recombinator::max_scaled_jacobian(MElement* element,int& index){
 
   if(j1>=j2 && j1>=j3 && j1>=j4){
     index = 0;
-	val = j1;
+    val = j1;
   }
   else if(j2>=j3 && j2>=j4 && j2>=j1){
     index = 1;
-	val = j2;
+    val = j2;
   }
   else if(j3>=j4 && j3>=j1 && j3>=j2){
     index = 2;
-	val = j3;
+    val = j3;
   }
   else{
     index = 3;
-	val = j4;
+    val = j4;
   }
 
   return val;
@@ -2246,8 +2246,8 @@ double Recombinator::min_scaled_jacobian(Hex hex){
   min = 1000000000.0;
   for(i=0;i<8;i++){
     if(jacobians[i]<=min){
-	  min = jacobians[i];
-	}
+      min = jacobians[i];
+    }
   }
 
   return min;
@@ -2331,9 +2331,9 @@ void Supplementary::execute(){
   for(it=model->firstRegion();it!=model->lastRegion();it++)
   {
     gr = *it;
-	if(gr->getNumMeshElements()>0){
-	  execute(gr);
-	}
+    if(gr->getNumMeshElements()>0){
+      execute(gr);
+    }
   }
 }
 
@@ -2360,43 +2360,43 @@ void Supplementary::execute(GRegion* gr){
   hash_tableC.clear();
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(eight(element)){
-	  a = element->getVertex(0);
-	  b = element->getVertex(1);
-	  c = element->getVertex(2);
-	  d = element->getVertex(3);
-	  e = element->getVertex(4);
-	  f = element->getVertex(5);
-	  g = element->getVertex(6);
-	  h = element->getVertex(7);
-
-	  build_hash_tableA(a,b,c,d);
-	  build_hash_tableA(e,f,g,h);
-	  build_hash_tableA(a,b,f,e);
-	  build_hash_tableA(b,c,g,f);
-	  build_hash_tableA(d,c,g,h);
-	  build_hash_tableA(d,a,e,h);
-
-	  build_hash_tableB(a,b,c,d);
-	  build_hash_tableB(e,f,g,h);
-	  build_hash_tableB(a,b,f,e);
-	  build_hash_tableB(b,c,g,f);
-	  build_hash_tableB(d,c,g,h);
-	  build_hash_tableB(d,a,e,h);
-
-	  build_hash_tableC(Diagonal(a,b));
-	  build_hash_tableC(Diagonal(b,c));
-	  build_hash_tableC(Diagonal(c,d));
-	  build_hash_tableC(Diagonal(d,a));
-	  build_hash_tableC(Diagonal(e,f));
-	  build_hash_tableC(Diagonal(f,g));
-	  build_hash_tableC(Diagonal(g,h));
-	  build_hash_tableC(Diagonal(h,e));
-	  build_hash_tableC(Diagonal(a,e));
-	  build_hash_tableC(Diagonal(b,f));
-	  build_hash_tableC(Diagonal(c,g));
-	  build_hash_tableC(Diagonal(d,h));
-	}
+    if(eight(element)){
+      a = element->getVertex(0);
+      b = element->getVertex(1);
+      c = element->getVertex(2);
+      d = element->getVertex(3);
+      e = element->getVertex(4);
+      f = element->getVertex(5);
+      g = element->getVertex(6);
+      h = element->getVertex(7);
+
+      build_hash_tableA(a,b,c,d);
+      build_hash_tableA(e,f,g,h);
+      build_hash_tableA(a,b,f,e);
+      build_hash_tableA(b,c,g,f);
+      build_hash_tableA(d,c,g,h);
+      build_hash_tableA(d,a,e,h);
+
+      build_hash_tableB(a,b,c,d);
+      build_hash_tableB(e,f,g,h);
+      build_hash_tableB(a,b,f,e);
+      build_hash_tableB(b,c,g,f);
+      build_hash_tableB(d,c,g,h);
+      build_hash_tableB(d,a,e,h);
+
+      build_hash_tableC(Diagonal(a,b));
+      build_hash_tableC(Diagonal(b,c));
+      build_hash_tableC(Diagonal(c,d));
+      build_hash_tableC(Diagonal(d,a));
+      build_hash_tableC(Diagonal(e,f));
+      build_hash_tableC(Diagonal(f,g));
+      build_hash_tableC(Diagonal(g,h));
+      build_hash_tableC(Diagonal(h,e));
+      build_hash_tableC(Diagonal(a,e));
+      build_hash_tableC(Diagonal(b,f));
+      build_hash_tableC(Diagonal(c,g));
+      build_hash_tableC(Diagonal(d,h));
+    }
   }
 
   std::sort(potential.begin(),potential.end());
@@ -2418,9 +2418,9 @@ void Supplementary::init_markings(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(four(element)){
-	  markings.insert(std::pair<MElement*,bool>(element,0));
-	}
+    if(four(element)){
+      markings.insert(std::pair<MElement*,bool>(element,0));
+    }
   }
 }
 
@@ -2443,42 +2443,42 @@ void Supplementary::patern(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(four(element)){
-	  for(j=0;j<4;j++){
-	    a = element->getVertex(j);
-		vertices[0] = element->getVertex((j+1)%4);
-		vertices[1] = element->getVertex((j+2)%4);
-		vertices[2] = element->getVertex((j+3)%4);
-		for(k=0;k<3;k++){
-		  b = vertices[k%3];
-		  c = vertices[(k+1)%3];
-		  d = vertices[(k+2)%3];
-		  already.clear();
-		  already.push_back(a);
-		  already.push_back(b);
-		  already.push_back(c);
-		  already.push_back(d);
-		  bin1.clear();
-		  bin2.clear();
-		  find(b,d,already,bin1);
-		  find(c,d,already,bin2);
-		  for(it1=bin1.begin();it1!=bin1.end();it1++){
-			p = *it1;
-			for(it2=bin2.begin();it2!=bin2.end();it2++){
-			  q = *it2;
-			  if(p!=q && linked(p,q)){
-			    prism = Prism(a,b,c,d,p,q);
-				quality = min_scaled_jacobian(prism);
-				prism.set_quality(quality);
-				if(valid(prism)){
-				  potential.push_back(prism);
-				}
-			  }
-			}
-		  }
-		}
-	  }
-	}
+    if(four(element)){
+      for(j=0;j<4;j++){
+        a = element->getVertex(j);
+        vertices[0] = element->getVertex((j+1)%4);
+        vertices[1] = element->getVertex((j+2)%4);
+        vertices[2] = element->getVertex((j+3)%4);
+        for(k=0;k<3;k++){
+          b = vertices[k%3];
+          c = vertices[(k+1)%3];
+          d = vertices[(k+2)%3];
+          already.clear();
+          already.push_back(a);
+          already.push_back(b);
+          already.push_back(c);
+          already.push_back(d);
+          bin1.clear();
+          bin2.clear();
+          find(b,d,already,bin1);
+          find(c,d,already,bin2);
+          for(it1=bin1.begin();it1!=bin1.end();it1++){
+            p = *it1;
+            for(it2=bin2.begin();it2!=bin2.end();it2++){
+              q = *it2;
+              if(p!=q && linked(p,q)){
+                prism = Prism(a,b,c,d,p,q);
+                quality = min_scaled_jacobian(prism);
+                prism.set_quality(quality);
+                if(valid(prism)){
+                  potential.push_back(prism);
+                }
+              }
+            }
+          }
+        }
+      }
+    }
   }
 }
 
@@ -2503,69 +2503,69 @@ void Supplementary::merge(GRegion* gr){
   for(i=0;i<potential.size();i++){
     prism = potential[i];
 
-	threshold = 0.15;
-	if(prism.get_quality()<threshold){
-	  break;
-	}
-
-	a = prism.get_a();
-	b = prism.get_b();
-	c = prism.get_c();
-	d = prism.get_d();
-	e = prism.get_e();
-	f = prism.get_f();
-
-	parts.clear();
-	find(a,prism,parts);
-	find(b,prism,parts);
-	find(c,prism,parts);
-	find(d,prism,parts);
-	find(e,prism,parts);
-	find(f,prism,parts);
-
-	flag = 1;
-	for(it=parts.begin();it!=parts.end();it++){
-	  element = *it;
-	  it2 = markings.find(element);
-	  if(it2->second==1 && !sliver(element,prism)){
-	    flag = 0;
-		break;
-	  }
-	}
-	if(!flag) continue;
-
-	if(!valid(prism,parts)){
-	  continue;
-	}
-
-	if(!conformityA(prism)){
-	  continue;
-	}
-
-	if(!conformityB(prism)){
-	  continue;
-	}
-
-	if(!conformityC(prism)){
-	  continue;
-	}
-
-	if(!faces_statuquo(prism)){
-	  continue;
-	}
-
-	//printf("%d - %d/%d - %f\n",count,i,(int)potential.size(),prism.get_quality());
-	quality = quality + prism.get_quality();
-	for(it=parts.begin();it!=parts.end();it++){
+    threshold = 0.15;
+    if(prism.get_quality()<threshold){
+      break;
+    }
+
+    a = prism.get_a();
+    b = prism.get_b();
+    c = prism.get_c();
+    d = prism.get_d();
+    e = prism.get_e();
+    f = prism.get_f();
+
+    parts.clear();
+    find(a,prism,parts);
+    find(b,prism,parts);
+    find(c,prism,parts);
+    find(d,prism,parts);
+    find(e,prism,parts);
+    find(f,prism,parts);
+
+    flag = 1;
+    for(it=parts.begin();it!=parts.end();it++){
+      element = *it;
+      it2 = markings.find(element);
+      if(it2->second==1 && !sliver(element,prism)){
+        flag = 0;
+        break;
+      }
+    }
+    if(!flag) continue;
+
+    if(!valid(prism,parts)){
+      continue;
+    }
+
+    if(!conformityA(prism)){
+      continue;
+    }
+
+    if(!conformityB(prism)){
+      continue;
+    }
+
+    if(!conformityC(prism)){
+      continue;
+    }
+
+    if(!faces_statuquo(prism)){
+      continue;
+    }
+
+    //printf("%d - %d/%d - %f\n",count,i,(int)potential.size(),prism.get_quality());
+    quality = quality + prism.get_quality();
+    for(it=parts.begin();it!=parts.end();it++){
       element = *it;
-	  it2 = markings.find(element);
-	  it2->second = 1;
-	}
-	gr->addPrism(new MPrism(a,b,c,d,e,f));
-	build_hash_tableA(prism);
-	build_hash_tableB(prism);
-	build_hash_tableC(prism);
-	count++;
+      it2 = markings.find(element);
+      it2->second = 1;
+    }
+    gr->addPrism(new MPrism(a,b,c,d,e,f));
+    build_hash_tableA(prism);
+    build_hash_tableB(prism);
+    build_hash_tableC(prism);
+    count++;
   }
 
   opt.clear();
@@ -2575,10 +2575,10 @@ void Supplementary::merge(GRegion* gr){
 
   for(i=0;i<opt.size();i++){
     element = (MElement*)(opt[i]);
-	it2 = markings.find(element);
-	if(it2->second==0){
-	  gr->tetrahedra.push_back(opt[i]);
-	}
+    it2 = markings.find(element);
+    if(it2->second==0){
+      gr->tetrahedra.push_back(opt[i]);
+    }
   }
 
   printf("prisms average quality (0->1) : %f\n",quality/count);
@@ -2590,7 +2590,7 @@ void Supplementary::rearrange(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	element->setVolumePositive();
+    element->setVolumePositive();
   }
 }
 
@@ -2607,15 +2607,15 @@ void Supplementary::statistics(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	volume = element->getVolume();
+    volume = element->getVolume();
 
-	if(six(element)){
-	  prism_nbr = prism_nbr + 1;
-	  prism_volume = prism_volume + volume;
-	}
+    if(six(element)){
+      prism_nbr = prism_nbr + 1;
+      prism_volume = prism_volume + volume;
+    }
 
-	all_nbr = all_nbr + 1;
-	all_volume = all_volume + volume;
+    all_nbr = all_nbr + 1;
+    all_volume = all_volume + volume;
   }
 
   printf("percentage of prisms (number) : %.2f\n",prism_nbr*100.0/all_nbr);
@@ -2640,17 +2640,17 @@ void Supplementary::build_tuples(GRegion* gr){
   {
     gf = *it;
 
-	for(i=0;i<gf->getNumMeshElements();i++){
-	  element = gf->getMeshElement(i);
+    for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);
 
-	  if(element->getNumVertices()==3){
-	    a = element->getVertex(0);
-		b = element->getVertex(1);
-		c = element->getVertex(2);
+      if(element->getNumVertices()==3){
+        a = element->getVertex(0);
+        b = element->getVertex(1);
+        c = element->getVertex(2);
 
-		tuples.insert(Tuple(a,b,c,element,gf));
-	  }
-	}
+        tuples.insert(Tuple(a,b,c,element,gf));
+      }
+    }
   }
 }
 
@@ -2668,18 +2668,18 @@ void Supplementary::modify_surfaces(GRegion* gr){
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
 
-	if(element->getNumVertices()==6){
-	  a = element->getVertex(0);
-	  b = element->getVertex(1);
-	  c = element->getVertex(2);
-	  d = element->getVertex(3);
-	  e = element->getVertex(4);
-	  f = element->getVertex(5);
-
-	  modify_surfaces(a,d,e,b);
-	  modify_surfaces(a,d,f,c);
-	  modify_surfaces(b,e,f,c);
-	}
+    if(element->getNumVertices()==6){
+      a = element->getVertex(0);
+      b = element->getVertex(1);
+      c = element->getVertex(2);
+      d = element->getVertex(3);
+      e = element->getVertex(4);
+      f = element->getVertex(5);
+
+      modify_surfaces(a,d,e,b);
+      modify_surfaces(a,d,f,c);
+      modify_surfaces(b,e,f,c);
+    }
   }
 
   faces = gr->faces();
@@ -2688,37 +2688,37 @@ void Supplementary::modify_surfaces(GRegion* gr){
   {
     gf = *it;
 
-	opt.clear();
+    opt.clear();
 
-	for(i=0;i<gf->getNumMeshElements();i++){
-	  element = gf->getMeshElement(i);
+    for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);
 
-	  if(element->getNumVertices()==3){
-	    it2 = triangles.find(element);
-		if(it2==triangles.end()){
-		  opt.push_back(element);
-		}
-	  }
-	}
+      if(element->getNumVertices()==3){
+        it2 = triangles.find(element);
+        if(it2==triangles.end()){
+          opt.push_back(element);
+        }
+      }
+    }
 
-	gf->triangles.clear();
+    gf->triangles.clear();
 
-	for(i=0;i<opt.size();i++){
-	  gf->triangles.push_back((MTriangle*)opt[i]);
-	}
+    for(i=0;i<opt.size();i++){
+      gf->triangles.push_back((MTriangle*)opt[i]);
+    }
   }
 }
 
 void Supplementary::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
   bool flag1,flag2;
   MElement *element1,*element2;
-  GFace *gf1,*gf2;
+  GFace *gf1;//,*gf2;
   Tuple tuple1,tuple2;
   std::multiset<Tuple>::iterator it1;
   std::multiset<Tuple>::iterator it2;
 
   gf1 = NULL;
-  gf2 = NULL;
+  //gf2 = NULL;
 
   tuple1 = Tuple(a,b,c);
   tuple2 = Tuple(c,d,a);
@@ -2731,37 +2731,37 @@ void Supplementary::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d)
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
     if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  element1 = it1->get_element();
-	  gf1 = it1->get_gf();
-	}
+      flag1 = 1;
+      element1 = it1->get_element();
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  element2 = it2->get_element();
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      element2 = it2->get_element();
+      //gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     triangles.insert(element1);
-	triangles.insert(element2);
+    triangles.insert(element2);
 
-	gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
+    gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
   }
 
   tuple1 = Tuple(a,b,d);
@@ -2775,37 +2775,37 @@ void Supplementary::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d)
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  element1 = it1->get_element();
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      element1 = it1->get_element();
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  element2 = it2->get_element();
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      element2 = it2->get_element();
+      //gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     triangles.insert(element1);
-	triangles.insert(element2);
+    triangles.insert(element2);
 
-	gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
+    gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
   }
 }
 
@@ -2955,9 +2955,9 @@ bool Supplementary::linked(MVertex* v1,MVertex* v2){
   if(it!=vertex_to_vertices.end()){
     for(it2=(it->second).begin();it2!=(it->second).end();it2++){
       if(*it2==v2){
-	    flag = 1;
-	    break;
-	  }
+        flag = 1;
+        break;
+      }
     }
   }
 
@@ -2987,24 +2987,24 @@ void Supplementary::find(MVertex* vertex,Prism prism,std::set<MElement*>& final)
   if(it!=vertex_to_tetrahedra.end()){
     for(it2=(it->second).begin();it2!=(it->second).end();it2++){
       a = (*it2)->getVertex(0);
-	  b = (*it2)->getVertex(1);
-	  c = (*it2)->getVertex(2);
-	  d = (*it2)->getVertex(3);
+      b = (*it2)->getVertex(1);
+      c = (*it2)->getVertex(2);
+      d = (*it2)->getVertex(3);
 
-	  flag1 = inclusion(a,prism);
-	  flag2 = inclusion(b,prism);
-	  flag3 = inclusion(c,prism);
-	  flag4 = inclusion(d,prism);
+      flag1 = inclusion(a,prism);
+      flag2 = inclusion(b,prism);
+      flag3 = inclusion(c,prism);
+      flag4 = inclusion(d,prism);
 
-	  if(flag1 && flag2 && flag3 && flag4){
-	    final.insert(*it2);
-	  }
+      if(flag1 && flag2 && flag3 && flag4){
+        final.insert(*it2);
+      }
     }
   }
 }
 
 void Supplementary::intersection(const std::set<MVertex*>& bin1,const std::set<MVertex*>& bin2,
-								const std::vector<MVertex*>& already,std::set<MVertex*>& final){
+                                const std::vector<MVertex*>& already,std::set<MVertex*>& final){
   size_t i;
   bool ok;
   std::set<MVertex*> temp;
@@ -3015,16 +3015,16 @@ void Supplementary::intersection(const std::set<MVertex*>& bin1,const std::set<M
   for(it=temp.begin();it!=temp.end();it++){
     ok = 1;
 
-	for(i=0;i<already.size();i++){
-	  if((*it)==already[i]){
-	    ok = 0;
-		break;
-	  }
-	}
+    for(i=0;i<already.size();i++){
+      if((*it)==already[i]){
+        ok = 0;
+        break;
+      }
+    }
 
-	if(ok){
-	  final.insert(*it);
-	}
+    if(ok){
+      final.insert(*it);
+    }
   }
 }
 
@@ -3068,19 +3068,19 @@ bool Supplementary::inclusion(MVertex* v1,MVertex* v2,MVertex* v3,const std::set
   for(it=bin.begin();it!=bin.end();it++){
     element = *it;
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
 
-	flag1 = inclusion(v1,a,b,c,d);
-	flag2 = inclusion(v2,a,b,c,d);
-	flag3 = inclusion(v3,a,b,c,d);
+    flag1 = inclusion(v1,a,b,c,d);
+    flag2 = inclusion(v2,a,b,c,d);
+    flag3 = inclusion(v3,a,b,c,d);
 
-	if(flag1 && flag2 && flag3){
-	  ok = 1;
-	  break;
-	}
+    if(flag1 && flag2 && flag3){
+      ok = 1;
+      break;
+    }
   }
 
   return ok;
@@ -3095,15 +3095,15 @@ bool Supplementary::inclusion(Facet facet){
 
   while(it!=hash_tableA.end()){
     if(facet.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(facet.same_vertices(*it)){
-	  flag = 1;
-	  break;
-	}
+    if(facet.same_vertices(*it)){
+      flag = 1;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   return flag;
@@ -3118,15 +3118,15 @@ bool Supplementary::inclusion(Diagonal diagonal){
 
   while(it!=hash_tableB.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 1;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 1;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   return flag;
@@ -3141,15 +3141,15 @@ bool Supplementary::duplicate(Diagonal diagonal){
 
   while(it!=hash_tableC.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 1;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 1;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   return flag;
@@ -3298,34 +3298,34 @@ bool Supplementary::faces_statuquo(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     if(gf1!=gf2){
-	  ok = 0;
-	}
+      ok = 0;
+    }
   }
 
   tuple1 = Tuple(a,b,d);
@@ -3339,34 +3339,34 @@ bool Supplementary::faces_statuquo(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     if(gf1!=gf2){
-	  ok = 0;
-	}
+      ok = 0;
+    }
   }
 
   return ok;
@@ -3384,27 +3384,27 @@ void Supplementary::build_vertex_to_vertices(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(four(element)){
-	  for(j=0;j<element->getNumVertices();j++){
-	    a = element->getVertex(j);
-	    b = element->getVertex((j+1)%4);
-	    c = element->getVertex((j+2)%4);
-	    d = element->getVertex((j+3)%4);
-
-	    it = vertex_to_vertices.find(a);
-	    if(it!=vertex_to_vertices.end()){
-	      it->second.insert(b);
-		  it->second.insert(c);
-		  it->second.insert(d);
-	    }
-	    else{
-	      bin.clear();
-		  bin.insert(b);
-		  bin.insert(c);
-		  bin.insert(d);
-		  vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(a,bin));
-	    }
-	  }
+    if(four(element)){
+      for(j=0;j<element->getNumVertices();j++){
+        a = element->getVertex(j);
+        b = element->getVertex((j+1)%4);
+        c = element->getVertex((j+2)%4);
+        d = element->getVertex((j+3)%4);
+
+        it = vertex_to_vertices.find(a);
+        if(it!=vertex_to_vertices.end()){
+          it->second.insert(b);
+          it->second.insert(c);
+          it->second.insert(d);
+        }
+        else{
+          bin.clear();
+          bin.insert(b);
+          bin.insert(c);
+          bin.insert(d);
+          vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(a,bin));
+        }
+      }
     }
   }
 }
@@ -3421,21 +3421,21 @@ void Supplementary::build_vertex_to_tetrahedra(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(four(element)){
-	  for(j=0;j<element->getNumVertices();j++){
-	    vertex = element->getVertex(j);
-
-		it = vertex_to_tetrahedra.find(vertex);
-		if(it!=vertex_to_tetrahedra.end()){
-		  it->second.insert(element);
-		}
-		else{
-		  bin.clear();
-		  bin.insert(element);
-		  vertex_to_tetrahedra.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
-		}
-	  }
-	}
+    if(four(element)){
+      for(j=0;j<element->getNumVertices();j++){
+        vertex = element->getVertex(j);
+
+        it = vertex_to_tetrahedra.find(vertex);
+        if(it!=vertex_to_tetrahedra.end()){
+          it->second.insert(element);
+        }
+        else{
+          bin.clear();
+          bin.insert(element);
+          vertex_to_tetrahedra.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
+        }
+      }
+    }
   }
 }
 
@@ -3471,15 +3471,15 @@ void Supplementary::build_hash_tableA(Facet facet){
 
   while(it!=hash_tableA.end()){
     if(facet.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(facet.same_vertices(*it)){
-	  flag = 0;
-	  break;
-	}
+    if(facet.same_vertices(*it)){
+      flag = 0;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   if(flag){
@@ -3517,15 +3517,15 @@ void Supplementary::build_hash_tableB(Diagonal diagonal){
 
   while(it!=hash_tableB.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 0;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 0;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   if(flag){
@@ -3564,15 +3564,15 @@ void Supplementary::build_hash_tableC(Diagonal diagonal){
 
   while(it!=hash_tableC.end()){
     if(diagonal.get_hash()!=it->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(diagonal.same_vertices(*it)){
-	  flag = 0;
-	  break;
-	}
+    if(diagonal.same_vertices(*it)){
+      flag = 0;
+      break;
+    }
 
-	it++;
+    it++;
   }
 
   if(flag){
@@ -3629,8 +3629,8 @@ double Supplementary::min_scaled_jacobian(Prism prism){
   min = 1000000000.0;
   for(i=0;i<6;i++){
     if(jacobians[i]<=min){
-	  min = jacobians[i];
-	}
+      min = jacobians[i];
+    }
   }
 
   return min;
@@ -3652,9 +3652,9 @@ void PostOp::execute(bool flag){
   for(it=model->firstRegion();it!=model->lastRegion();it++)
   {
     gr = *it;
-	if(gr->getNumMeshElements()>0){
-	  execute(gr,flag);
-	}
+    if(gr->getNumMeshElements()>0){
+      execute(gr,flag);
+    }
   }
 }
 
@@ -3691,9 +3691,9 @@ void PostOp::init_markings(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(four(element) || five(element)){
-	  markings.insert(std::pair<MElement*,bool>(element,0));
-	}
+    if(four(element) || five(element)){
+      markings.insert(std::pair<MElement*,bool>(element,0));
+    }
   }
 }
 
@@ -3712,47 +3712,47 @@ void PostOp::pyramids1(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(eight(element)){
-	  hexahedra.push_back(element);
-	}
-	else if(six(element)){
-	  prisms.push_back(element);
-	}
+    if(eight(element)){
+      hexahedra.push_back(element);
+    }
+    else if(six(element)){
+      prisms.push_back(element);
+    }
   }
 
   for(i=0;i<hexahedra.size();i++){
     element = hexahedra[i];
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
-	e = element->getVertex(4);
-	f = element->getVertex(5);
-	g = element->getVertex(6);
-	h = element->getVertex(7);
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
+    e = element->getVertex(4);
+    f = element->getVertex(5);
+    g = element->getVertex(6);
+    h = element->getVertex(7);
 
-	pyramids1(a,b,c,d,gr);
-	pyramids1(e,f,g,h,gr);
-	pyramids1(a,b,f,e,gr);
-	pyramids1(b,c,g,f,gr);
-	pyramids1(d,c,g,h,gr);
-	pyramids1(d,a,e,h,gr);
+    pyramids1(a,b,c,d,gr);
+    pyramids1(e,f,g,h,gr);
+    pyramids1(a,b,f,e,gr);
+    pyramids1(b,c,g,f,gr);
+    pyramids1(d,c,g,h,gr);
+    pyramids1(d,a,e,h,gr);
   }
 
   for(i=0;i<prisms.size();i++){
     element = prisms[i];
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
-	e = element->getVertex(4);
-	f = element->getVertex(5);
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
+    e = element->getVertex(4);
+    f = element->getVertex(5);
 
-	pyramids1(a,d,f,c,gr);
-	pyramids1(a,b,e,d,gr);
-	pyramids1(b,c,f,e,gr);
+    pyramids1(a,d,f,c,gr);
+    pyramids1(a,b,e,d,gr);
+    pyramids1(b,c,f,e,gr);
   }
 
   opt.clear();
@@ -3762,10 +3762,10 @@ void PostOp::pyramids1(GRegion* gr){
 
   for(i=0;i<opt.size();i++){
     element = (MElement*)(opt[i]);
-	it = markings.find(element);
-	if(it->second==0){
-	  gr->tetrahedra.push_back(opt[i]);
-	}
+    it = markings.find(element);
+    if(it->second==0){
+      gr->tetrahedra.push_back(opt[i]);
+    }
   }
 }
 
@@ -3785,47 +3785,47 @@ void PostOp::pyramids2(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(eight(element)){
-	  hexahedra.push_back(element);
-	}
-	else if(six(element)){
-	  prisms.push_back(element);
-	}
+    if(eight(element)){
+      hexahedra.push_back(element);
+    }
+    else if(six(element)){
+      prisms.push_back(element);
+    }
   }
 
   for(i=0;i<hexahedra.size();i++){
     element = hexahedra[i];
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
-	e = element->getVertex(4);
-	f = element->getVertex(5);
-	g = element->getVertex(6);
-	h = element->getVertex(7);
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
+    e = element->getVertex(4);
+    f = element->getVertex(5);
+    g = element->getVertex(6);
+    h = element->getVertex(7);
 
-	pyramids2(a,b,c,d,gr);
-	pyramids2(e,f,g,h,gr);
-	pyramids2(a,b,f,e,gr);
-	pyramids2(b,c,g,f,gr);
-	pyramids2(d,c,g,h,gr);
-	pyramids2(d,a,e,h,gr);
+    pyramids2(a,b,c,d,gr);
+    pyramids2(e,f,g,h,gr);
+    pyramids2(a,b,f,e,gr);
+    pyramids2(b,c,g,f,gr);
+    pyramids2(d,c,g,h,gr);
+    pyramids2(d,a,e,h,gr);
   }
 
   for(i=0;i<prisms.size();i++){
     element = prisms[i];
 
-	a = element->getVertex(0);
-	b = element->getVertex(1);
-	c = element->getVertex(2);
-	d = element->getVertex(3);
-	e = element->getVertex(4);
-	f = element->getVertex(5);
+    a = element->getVertex(0);
+    b = element->getVertex(1);
+    c = element->getVertex(2);
+    d = element->getVertex(3);
+    e = element->getVertex(4);
+    f = element->getVertex(5);
 
-	pyramids2(a,d,f,c,gr);
-	pyramids2(a,b,e,d,gr);
-	pyramids2(b,c,f,e,gr);
+    pyramids2(a,d,f,c,gr);
+    pyramids2(a,b,e,d,gr);
+    pyramids2(b,c,f,e,gr);
   }
 
   opt1.clear();
@@ -3835,10 +3835,10 @@ void PostOp::pyramids2(GRegion* gr){
 
   for(i=0;i<opt1.size();i++){
     element = (MElement*)(opt1[i]);
-	it = markings.find(element);
-	if(it->second==0){
-	  gr->tetrahedra.push_back(opt1[i]);
-	}
+    it = markings.find(element);
+    if(it->second==0){
+      gr->tetrahedra.push_back(opt1[i]);
+    }
   }
 
   opt2.clear();
@@ -3848,10 +3848,10 @@ void PostOp::pyramids2(GRegion* gr){
 
   for(i=0;i<opt2.size();i++){
     element = (MElement*)(opt2[i]);
-	it = markings.find(element);
-	if(it->second==0){
-	  gr->pyramids.push_back(opt2[i]);
-	}
+    it = markings.find(element);
+    if(it->second==0){
+      gr->pyramids.push_back(opt2[i]);
+    }
   }
 }
 
@@ -3879,19 +3879,19 @@ void PostOp::pyramids1(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
 
   if(bin.size()==2){
     it = bin.begin();
-	it1 = markings.find(*it);
-	it++;
-	it2 = markings.find(*it);
+    it1 = markings.find(*it);
+    it++;
+    it2 = markings.find(*it);
 
-	if(it1->second==0 && it2->second==0){
-	  vertex = find(a,b,c,d,*it);
+    if(it1->second==0 && it2->second==0){
+      vertex = find(a,b,c,d,*it);
 
-	  if(vertex!=0){
-		gr->addPyramid(new MPyramid(a,b,c,d,vertex));
-		it1->second = 1;
-		it2->second = 1;
-	  }
-	}
+      if(vertex!=0){
+        gr->addPyramid(new MPyramid(a,b,c,d,vertex));
+        it1->second = 1;
+        it2->second = 1;
+      }
+    }
   }
 }
 
@@ -3946,19 +3946,19 @@ void PostOp::pyramids2(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
   }
 
   /*if(pyramids.size()==0 && tetrahedra.size()==1){
-	printf("tetrahedron deleted\n");
+    printf("tetrahedron deleted\n");
     it2 = markings.find(*tetrahedra.begin());
-	it2->second = 1;
-	return;
+    it2->second = 1;
+    return;
   }*/
 
   if(flag){
     diagA = a;
-	diagB = c;
+    diagB = c;
   }
   else{
     diagA = b;
-	diagB = d;
+    diagB = d;
   }
 
   Ns.clear();
@@ -3973,109 +3973,109 @@ void PostOp::pyramids2(MVertex* a,MVertex* b,MVertex* c,MVertex* d,GRegion* gr){
   movables.clear();
 
   if(tetrahedra.size()>0 || pyramids.size()>0){
-	estimate1 = estimate1 + tetrahedra.size() + 2*pyramids.size();
-	estimate2 = estimate2 + 1;
-
-	mid = new MVertex(x,y,z,gr);
-	gr->addMeshVertex(mid);
-
-	temp2 = new MPyramid(a,b,c,d,mid);
-	gr->addPyramid(temp2);
-	markings.insert(std::pair<MElement*,bool>(temp2,0));
-	build_vertex_to_pyramids(temp2);
-
-	for(it=tetrahedra.begin();it!=tetrahedra.end();it++){
-	  N1 = other(*it,diagA,diagB);
-	  N2 = other(*it,diagA,diagB,N1);
-
-	  if(N1!=0 && N2!=0){
-		Ns.insert(N1);
-		Ns.insert(N2);
-
-	    temp = new MTetrahedron(N1,N2,diagA,mid);
-		gr->addTetrahedron(temp);
-		markings.insert(std::pair<MElement*,bool>(temp,0));
-		build_vertex_to_tetrahedra(temp);
-		movables.push_back(temp);
-
-		temp = new MTetrahedron(N1,N2,diagB,mid);
-		gr->addTetrahedron(temp);
-		markings.insert(std::pair<MElement*,bool>(temp,0));
-		build_vertex_to_tetrahedra(temp);
-		movables.push_back(temp);
-
-		it2 = markings.find(*it);
-		it2->second = 1;
-		erase_vertex_to_tetrahedra(*it);
-	  }
-	}
-
-	for(it=pyramids.begin();it!=pyramids.end();it++){
-	  v1 = (*it)->getVertex(0);
-	  v2 = (*it)->getVertex(1);
-	  v3 = (*it)->getVertex(2);
-	  v4 = (*it)->getVertex(3);
-	  v5 = (*it)->getVertex(4);
-
-	  if(v1!=diagA && v1!=diagB){
-	    Ns.insert(v1);
-	  }
-	  if(v2!=diagA && v2!=diagB){
-	    Ns.insert(v2);
-	  }
-	  if(v3!=diagA && v3!=diagB){
-	    Ns.insert(v3);
-	  }
-	  if(v4!=diagA && v4!=diagB){
-	    Ns.insert(v4);
-	  }
-	  if(v5!=diagA && v5!=diagB){
-	    Ns.insert(v5);
-	  }
-
-	  temp2 = new MPyramid(v1,v2,v3,v4,mid);
-	  gr->addPyramid(temp2);
-	  markings.insert(std::pair<MElement*,bool>(temp2,0));
-	  build_vertex_to_pyramids(temp2);
-
-	  if(different(v1,v2,diagA,diagB)){
-	    temp = new MTetrahedron(v1,v2,mid,v5);
-		gr->addTetrahedron(temp);
-		markings.insert(std::pair<MElement*,bool>(temp,0));
-		build_vertex_to_tetrahedra(temp);
-		movables.push_back(temp);
-	  }
-
-	  if(different(v2,v3,diagA,diagB)){
-	    temp = new MTetrahedron(v2,v3,mid,v5);
-		gr->addTetrahedron(temp);
-		markings.insert(std::pair<MElement*,bool>(temp,0));
-		build_vertex_to_tetrahedra(temp);
-		movables.push_back(temp);
-	  }
-
-	  if(different(v3,v4,diagA,diagB)){
-	    temp = new MTetrahedron(v3,v4,mid,v5);
-		gr->addTetrahedron(temp);
-		markings.insert(std::pair<MElement*,bool>(temp,0));
-		build_vertex_to_tetrahedra(temp);
-		movables.push_back(temp);
-	  }
-
-	  if(different(v4,v1,diagA,diagB)){
-	    temp = new MTetrahedron(v4,v1,mid,v5);
-		gr->addTetrahedron(temp);
-		markings.insert(std::pair<MElement*,bool>(temp,0));
-		build_vertex_to_tetrahedra(temp);
-		movables.push_back(temp);
-	  }
-
-	  it2 = markings.find(*it);
-	  it2->second = 1;
-	  erase_vertex_to_pyramids(*it);
-	}
-
-	mean(Ns,mid,movables);
+    estimate1 = estimate1 + tetrahedra.size() + 2*pyramids.size();
+    estimate2 = estimate2 + 1;
+
+    mid = new MVertex(x,y,z,gr);
+    gr->addMeshVertex(mid);
+
+    temp2 = new MPyramid(a,b,c,d,mid);
+    gr->addPyramid(temp2);
+    markings.insert(std::pair<MElement*,bool>(temp2,0));
+    build_vertex_to_pyramids(temp2);
+
+    for(it=tetrahedra.begin();it!=tetrahedra.end();it++){
+      N1 = other(*it,diagA,diagB);
+      N2 = other(*it,diagA,diagB,N1);
+
+      if(N1!=0 && N2!=0){
+        Ns.insert(N1);
+        Ns.insert(N2);
+
+        temp = new MTetrahedron(N1,N2,diagA,mid);
+        gr->addTetrahedron(temp);
+        markings.insert(std::pair<MElement*,bool>(temp,0));
+        build_vertex_to_tetrahedra(temp);
+        movables.push_back(temp);
+
+        temp = new MTetrahedron(N1,N2,diagB,mid);
+        gr->addTetrahedron(temp);
+        markings.insert(std::pair<MElement*,bool>(temp,0));
+        build_vertex_to_tetrahedra(temp);
+        movables.push_back(temp);
+
+        it2 = markings.find(*it);
+        it2->second = 1;
+        erase_vertex_to_tetrahedra(*it);
+      }
+    }
+
+    for(it=pyramids.begin();it!=pyramids.end();it++){
+      v1 = (*it)->getVertex(0);
+      v2 = (*it)->getVertex(1);
+      v3 = (*it)->getVertex(2);
+      v4 = (*it)->getVertex(3);
+      v5 = (*it)->getVertex(4);
+
+      if(v1!=diagA && v1!=diagB){
+        Ns.insert(v1);
+      }
+      if(v2!=diagA && v2!=diagB){
+        Ns.insert(v2);
+      }
+      if(v3!=diagA && v3!=diagB){
+        Ns.insert(v3);
+      }
+      if(v4!=diagA && v4!=diagB){
+        Ns.insert(v4);
+      }
+      if(v5!=diagA && v5!=diagB){
+        Ns.insert(v5);
+      }
+
+      temp2 = new MPyramid(v1,v2,v3,v4,mid);
+      gr->addPyramid(temp2);
+      markings.insert(std::pair<MElement*,bool>(temp2,0));
+      build_vertex_to_pyramids(temp2);
+
+      if(different(v1,v2,diagA,diagB)){
+        temp = new MTetrahedron(v1,v2,mid,v5);
+        gr->addTetrahedron(temp);
+        markings.insert(std::pair<MElement*,bool>(temp,0));
+        build_vertex_to_tetrahedra(temp);
+        movables.push_back(temp);
+      }
+
+      if(different(v2,v3,diagA,diagB)){
+        temp = new MTetrahedron(v2,v3,mid,v5);
+        gr->addTetrahedron(temp);
+        markings.insert(std::pair<MElement*,bool>(temp,0));
+        build_vertex_to_tetrahedra(temp);
+        movables.push_back(temp);
+      }
+
+      if(different(v3,v4,diagA,diagB)){
+        temp = new MTetrahedron(v3,v4,mid,v5);
+        gr->addTetrahedron(temp);
+        markings.insert(std::pair<MElement*,bool>(temp,0));
+        build_vertex_to_tetrahedra(temp);
+        movables.push_back(temp);
+      }
+
+      if(different(v4,v1,diagA,diagB)){
+        temp = new MTetrahedron(v4,v1,mid,v5);
+        gr->addTetrahedron(temp);
+        markings.insert(std::pair<MElement*,bool>(temp,0));
+        build_vertex_to_tetrahedra(temp);
+        movables.push_back(temp);
+      }
+
+      it2 = markings.find(*it);
+      it2->second = 1;
+      erase_vertex_to_pyramids(*it);
+    }
+
+    mean(Ns,mid,movables);
   }
 }
 
@@ -4085,7 +4085,7 @@ void PostOp::rearrange(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	element->setVolumePositive();
+    element->setVolumePositive();
   }
 }
 
@@ -4109,34 +4109,34 @@ void PostOp::statistics(GRegion* gr){
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
 
-	if(eight(element)){
-	  nbr8 = nbr8 + 1;
-	  vol8 = vol8 + element->getVolume();
-	}
+    if(eight(element)){
+      nbr8 = nbr8 + 1;
+      vol8 = vol8 + element->getVolume();
+    }
 
-	if(six(element)){
-	  nbr6 = nbr6 + 1;
-	  vol6 = vol6 + element->getVolume();
-	}
+    if(six(element)){
+      nbr6 = nbr6 + 1;
+      vol6 = vol6 + element->getVolume();
+    }
 
     if(five(element)){
-	  nbr5 = nbr5 + 1;
-	  vol5 = vol5 + workaround(element);
-	}
+      nbr5 = nbr5 + 1;
+      vol5 = vol5 + workaround(element);
+    }
 
-	if(four(element)){
-	  nbr4 = nbr4 + 1;
-	  vol4 = vol4 + element->getVolume();
-	}
+    if(four(element)){
+      nbr4 = nbr4 + 1;
+      vol4 = vol4 + element->getVolume();
+    }
 
-	nbr = nbr + 1;
+    nbr = nbr + 1;
 
-	if(!five(element)){
-	  vol = vol + element->getVolume();
-	}
-	else{
-	  vol = vol + workaround(element);
-	}
+    if(!five(element)){
+      vol = vol + element->getVolume();
+    }
+    else{
+      vol = vol + workaround(element);
+    }
   }
 
   printf("Number :\n");
@@ -4172,23 +4172,23 @@ void PostOp::build_tuples(GRegion* gr){
   {
     gf = *it;
 
-	for(i=0;i<gf->getNumMeshElements();i++){
-	  element = gf->getMeshElement(i);
+    for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);
 
-	  if(element->getNumVertices()==3){
-	    a = element->getVertex(0);
-		b = element->getVertex(1);
-		c = element->getVertex(2);
+      if(element->getNumVertices()==3){
+        a = element->getVertex(0);
+        b = element->getVertex(1);
+        c = element->getVertex(2);
 
-		tuples.insert(Tuple(a,b,c,element,gf));
-	  }
-	}
+        tuples.insert(Tuple(a,b,c,element,gf));
+      }
+    }
   }
 }
 
 void PostOp::modify_surfaces(GRegion* gr){
   unsigned int i;
-  MVertex *a,*b,*c,*d,*e;
+  MVertex *a,*b,*c,*d;//,*e;
   MElement* element;
   GFace* gf;
   std::list<GFace*> faces;
@@ -4199,15 +4199,15 @@ void PostOp::modify_surfaces(GRegion* gr){
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
 
-	if(element->getNumVertices()==5){
-	  a = element->getVertex(0);
-	  b = element->getVertex(1);
-	  c = element->getVertex(2);
-	  d = element->getVertex(3);
-	  e = element->getVertex(4);
+    if(element->getNumVertices()==5){
+      a = element->getVertex(0);
+      b = element->getVertex(1);
+      c = element->getVertex(2);
+      d = element->getVertex(3);
+      //e = element->getVertex(4);
 
-	  modify_surfaces(a,b,c,d);
-	}
+      modify_surfaces(a,b,c,d);
+    }
   }
 
   faces = gr->faces();
@@ -4216,37 +4216,37 @@ void PostOp::modify_surfaces(GRegion* gr){
   {
     gf = *it;
 
-	opt.clear();
+    opt.clear();
 
-	for(i=0;i<gf->getNumMeshElements();i++){
-	  element = gf->getMeshElement(i);
+    for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);
 
-	  if(element->getNumVertices()==3){
-	    it2 = triangles.find(element);
-		if(it2==triangles.end()){
-		  opt.push_back(element);
-		}
-	  }
-	}
+      if(element->getNumVertices()==3){
+        it2 = triangles.find(element);
+        if(it2==triangles.end()){
+          opt.push_back(element);
+        }
+      }
+    }
 
-	gf->triangles.clear();
+    gf->triangles.clear();
 
-	for(i=0;i<opt.size();i++){
-	  gf->triangles.push_back((MTriangle*)opt[i]);
-	}
+    for(i=0;i<opt.size();i++){
+      gf->triangles.push_back((MTriangle*)opt[i]);
+    }
   }
 }
 
 void PostOp::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
   bool flag1,flag2;
   MElement *element1,*element2;
-  GFace *gf1,*gf2;
+  GFace *gf1;//,*gf2;
   Tuple tuple1,tuple2;
   std::multiset<Tuple>::iterator it1;
   std::multiset<Tuple>::iterator it2;
 
   gf1 = NULL;
-  gf2 = NULL;
+  //gf2 = NULL;
 
   tuple1 = Tuple(a,b,c);
   tuple2 = Tuple(c,d,a);
@@ -4259,37 +4259,37 @@ void PostOp::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  element1 = it1->get_element();
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      element1 = it1->get_element();
+      gf1 = it1->get_gf();
+    }
 
-	it1++;
+    it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
-	  flag2 = 1;
-	  element2 = it2->get_element();
-	  gf2 = it2->get_gf();
-	}
+    if(tuple2.same_vertices(*it2)){
+      flag2 = 1;
+      element2 = it2->get_element();
+      //gf2 = it2->get_gf();
+    }
 
-	it2++;
+    it2++;
   }
 
   if(flag1 && flag2){
     triangles.insert(element1);
-	triangles.insert(element2);
+    triangles.insert(element2);
 
-	gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
+    gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
   }
 
   tuple1 = Tuple(a,b,d);
@@ -4303,37 +4303,37 @@ void PostOp::modify_surfaces(MVertex* a,MVertex* b,MVertex* c,MVertex* d){
 
   while(it1!=tuples.end()){
     if(tuple1.get_hash()!=it1->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple1.same_vertices(*it1)){
-	  flag1 = 1;
-	  element1 = it1->get_element();
-	  gf1 = it1->get_gf();
-	}
+    if(tuple1.same_vertices(*it1)){
+      flag1 = 1;
+      element1 = it1->get_element();
+      gf1 = it1->get_gf();
+    }
 
     it1++;
   }
 
   while(it2!=tuples.end()){
     if(tuple2.get_hash()!=it2->get_hash()){
-	  break;
-	}
+      break;
+    }
 
-	if(tuple2.same_vertices(*it2)){
+    if(tuple2.same_vertices(*it2)){
       flag2 = 1;
-	  element2 = it2->get_element();
-	  gf2 = it2->get_gf();
-	}
+      element2 = it2->get_element();
+      //gf2 = it2->get_gf();
+    }
 
     it2++;
   }
 
   if(flag1 && flag2){
     triangles.insert(element1);
-	triangles.insert(element2);
+    triangles.insert(element2);
 
-	gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
+    gf1->addQuadrangle(new MQuadrangle(a,b,c,d));
   }
 }
 
@@ -4384,10 +4384,10 @@ MVertex* PostOp::other(MElement* element,MVertex* v1,MVertex* v2){
 
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
-	if(vertex!=v1 && vertex!=v2){
-	  pointer = vertex;
-	  break;
-	}
+    if(vertex!=v1 && vertex!=v2){
+      pointer = vertex;
+      break;
+    }
   }
 
   return pointer;
@@ -4402,10 +4402,10 @@ MVertex* PostOp::other(MElement* element,MVertex* v1,MVertex* v2,MVertex* v3){
 
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
-	if(vertex!=v1 && vertex!=v2 && vertex!=v3){
-	  pointer = vertex;
-	  break;
-	}
+    if(vertex!=v1 && vertex!=v2 && vertex!=v3){
+      pointer = vertex;
+      break;
+    }
   }
 
   return pointer;
@@ -4429,8 +4429,8 @@ void PostOp::mean(const std::set<MVertex*>& Ns,MVertex* mid,const std::vector<ME
 
   for(it=Ns.begin();it!=Ns.end();it++){
     x = x + (*it)->x();
-	y = y + (*it)->y();
-	z = z + (*it)->z();
+    y = y + (*it)->y();
+    z = z + (*it)->z();
   }
 
   x = x/Ns.size();
@@ -4444,23 +4444,23 @@ void PostOp::mean(const std::set<MVertex*>& Ns,MVertex* mid,const std::vector<ME
   mid->setXYZ(x,y,z);
 
   for(j=0;j<100;j++){
-	flag = 0;
+    flag = 0;
 
-	for(i=0;i<movables.size();i++){
-	  if(movables[i]->getVolume()<0.0){
-	    flag = 1;
-	  }
-	}
+    for(i=0;i<movables.size();i++){
+      if(movables[i]->getVolume()<0.0){
+        flag = 1;
+      }
+    }
 
-	if(!flag){
-	  break;
-	}
+    if(!flag){
+      break;
+    }
 
-	x = 0.1*init_x + 0.9*mid->x();
-	y = 0.1*init_y + 0.9*mid->y();
-	z = 0.1*init_z + 0.9*mid->z();
+    x = 0.1*init_x + 0.9*mid->x();
+    y = 0.1*init_y + 0.9*mid->y();
+    z = 0.1*init_z + 0.9*mid->z();
 
-	mid->setXYZ(x,y,z);
+    mid->setXYZ(x,y,z);
   }
 
   iterations = iterations + j;
@@ -4468,21 +4468,21 @@ void PostOp::mean(const std::set<MVertex*>& Ns,MVertex* mid,const std::vector<ME
   for(j=0;j<6;j++){
     flag = 0;
 
-	for(i=0;i<movables.size();i++){
-	  if(movables[i]->gammaShapeMeasure()<0.2){
-	    flag = 1;
-	  }
-	}
+    for(i=0;i<movables.size();i++){
+      if(movables[i]->gammaShapeMeasure()<0.2){
+        flag = 1;
+      }
+    }
 
-	if(!flag){
-	  break;
-	}
+    if(!flag){
+      break;
+    }
 
-	x = 0.1*init_x + 0.9*mid->x();
-	y = 0.1*init_y + 0.9*mid->y();
-	z = 0.1*init_z + 0.9*mid->z();
+    x = 0.1*init_x + 0.9*mid->x();
+    y = 0.1*init_y + 0.9*mid->y();
+    z = 0.1*init_z + 0.9*mid->z();
 
-	mid->setXYZ(x,y,z);
+    mid->setXYZ(x,y,z);
   }
 
   iterations = iterations + j;
@@ -4497,10 +4497,10 @@ double PostOp::workaround(MElement* element){
 
   if(five(element)){
     temp1 = new MTetrahedron(element->getVertex(0),element->getVertex(1),element->getVertex(2),element->getVertex(4));
-	temp2 = new MTetrahedron(element->getVertex(2),element->getVertex(3),element->getVertex(0),element->getVertex(4));
-	volume = fabs(temp1->getVolume()) + fabs(temp2->getVolume());
-	delete temp1;
-	delete temp2;
+    temp2 = new MTetrahedron(element->getVertex(2),element->getVertex(3),element->getVertex(0),element->getVertex(4));
+    volume = fabs(temp1->getVolume()) + fabs(temp2->getVolume());
+    delete temp1;
+    delete temp2;
   }
 
   return volume;
@@ -4515,10 +4515,10 @@ MVertex* PostOp::find(MVertex* v1,MVertex* v2,MVertex* v3,MVertex* v4,MElement*
 
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
-	if(vertex!=v1 && vertex!=v2 && vertex!=v3 && vertex!=v4){
-	  pointer = vertex;
-	  break;
-	}
+    if(vertex!=v1 && vertex!=v2 && vertex!=v3 && vertex!=v4){
+      pointer = vertex;
+      break;
+    }
   }
 
   return pointer;
@@ -4555,16 +4555,16 @@ void PostOp::find_pyramids(MVertex* v1,MVertex* v2,std::set<MElement*>& final){
 
   for(it=temp.begin();it!=temp.end();it++){
     flag1 = equal(v1,v2,(*it)->getVertex(0),(*it)->getVertex(1));
-	flag2 = equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(2));
-	flag3 = equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(3));
-	flag4 = equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(0));
-	flag5 = equal(v1,v2,(*it)->getVertex(0),(*it)->getVertex(4));
-	flag6 = equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(4));
-	flag7 = equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(4));
-	flag8 = equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(4));
-	if(flag1 || flag2 || flag3 || flag4 || flag5 || flag6 || flag7 || flag8){
-	  final.insert(*it);
-	}
+    flag2 = equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(2));
+    flag3 = equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(3));
+    flag4 = equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(0));
+    flag5 = equal(v1,v2,(*it)->getVertex(0),(*it)->getVertex(4));
+    flag6 = equal(v1,v2,(*it)->getVertex(1),(*it)->getVertex(4));
+    flag7 = equal(v1,v2,(*it)->getVertex(2),(*it)->getVertex(4));
+    flag8 = equal(v1,v2,(*it)->getVertex(3),(*it)->getVertex(4));
+    if(flag1 || flag2 || flag3 || flag4 || flag5 || flag6 || flag7 || flag8){
+      final.insert(*it);
+    }
   }
 }
 
@@ -4580,9 +4580,9 @@ void PostOp::build_vertex_to_tetrahedra(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(four(element)){
-	  build_vertex_to_tetrahedra(element);
-	}
+    if(four(element)){
+      build_vertex_to_tetrahedra(element);
+    }
   }
 }
 
@@ -4595,15 +4595,15 @@ void PostOp::build_vertex_to_tetrahedra(MElement* element){
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
 
-	it = vertex_to_tetrahedra.find(vertex);
-	if(it!=vertex_to_tetrahedra.end()){
-	  it->second.insert(element);
-	}
-	else{
-	  bin.clear();
-	  bin.insert(element);
-	  vertex_to_tetrahedra.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
-	}
+    it = vertex_to_tetrahedra.find(vertex);
+    if(it!=vertex_to_tetrahedra.end()){
+      it->second.insert(element);
+    }
+    else{
+      bin.clear();
+      bin.insert(element);
+      vertex_to_tetrahedra.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
+    }
   }
 }
 
@@ -4615,10 +4615,10 @@ void PostOp::erase_vertex_to_tetrahedra(MElement* element){
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
 
-	it = vertex_to_tetrahedra.find(vertex);
-	if(it!=vertex_to_tetrahedra.end()){
-	  it->second.erase(element);
-	}
+    it = vertex_to_tetrahedra.find(vertex);
+    if(it!=vertex_to_tetrahedra.end()){
+      it->second.erase(element);
+    }
   }
 }
 
@@ -4630,9 +4630,9 @@ void PostOp::build_vertex_to_pyramids(GRegion* gr){
 
   for(i=0;i<gr->getNumMeshElements();i++){
     element = gr->getMeshElement(i);
-	if(five(element)){
-	  build_vertex_to_pyramids(element);
-	}
+    if(five(element)){
+      build_vertex_to_pyramids(element);
+    }
   }
 }
 
@@ -4645,15 +4645,15 @@ void PostOp::build_vertex_to_pyramids(MElement* element){
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
 
-	it = vertex_to_pyramids.find(vertex);
-	if(it!=vertex_to_pyramids.end()){
-	  it->second.insert(element);
-	}
-	else{
-	  bin.clear();
-	  bin.insert(element);
-	  vertex_to_pyramids.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
-	}
+    it = vertex_to_pyramids.find(vertex);
+    if(it!=vertex_to_pyramids.end()){
+      it->second.insert(element);
+    }
+    else{
+      bin.clear();
+      bin.insert(element);
+      vertex_to_pyramids.insert(std::pair<MVertex*,std::set<MElement*> >(vertex,bin));
+    }
   }
 }
 
@@ -4665,9 +4665,9 @@ void PostOp::erase_vertex_to_pyramids(MElement* element){
   for(i=0;i<element->getNumVertices();i++){
     vertex = element->getVertex(i);
 
-	it = vertex_to_pyramids.find(vertex);
-	if(it!=vertex_to_pyramids.end()){
-	  it->second.erase(element);
-	}
+    it = vertex_to_pyramids.find(vertex);
+    if(it!=vertex_to_pyramids.end()){
+      it->second.erase(element);
+    }
   }
 }
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index bb948ae2d33263dbecab55be02a97a58318a72eb..35c7f4fe0cb1665156602d3189010d7e81fc94f6 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -963,13 +963,14 @@ void   skipline(void);
 
 #define YY_INPUT(buf,result,max_size)					\
      {									\
-       int c = '*', n;							\
+       int c = '*';							\
+       unsigned int n;                                                  \
        for ( n = 0; n < max_size &&					\
-	       (c = getc( gmsh_yyin )) != EOF && c != '\n'; ++n )		\
+	       (c = getc( gmsh_yyin )) != EOF && c != '\n'; ++n )       \
 	 buf[n] = (char) c;						\
        if ( c == '\n' ){						\
 	 buf[n++] = (char) c;						\
-	 gmsh_yylineno++;							\
+	 gmsh_yylineno++;						\
        }								\
        if ( c == EOF && ferror( gmsh_yyin ) )				\
 	 Msg::Fatal("Input in flex scanner failed");			\
diff --git a/contrib/Chaco/input/input.c b/contrib/Chaco/input/input.c
index 23db3260fb3f7a1adc14ba63b637abc15e4c8702..abab4eb607f486acd51e203b83e3a7094386a71a 100644
--- a/contrib/Chaco/input/input.c
+++ b/contrib/Chaco/input/input.c
@@ -53,7 +53,7 @@ int      *ndims;		/* number of divisions at each stage */
     while (*fin == NULL) {
 	if (PROMPT)
 	    printf("Graph input file: ");
-	scanf("%s", inname);
+	if(scanf("%s", inname) != 1) return;
 
 	*fin = fopen(inname, "r");
 	if (*fin == NULL) {
@@ -65,14 +65,14 @@ int      *ndims;		/* number of divisions at each stage */
     if (OUTPUT_ASSIGN && !SEQUENCE) {
 	if (PROMPT)
 	    printf("Assignment output file: ");
-	scanf("%s", outassignname);
+	if(scanf("%s", outassignname) != 1) return;
     }
 
     /* Name output results file. */
     if (ECHO < 0) {
 	if (PROMPT)
 	    printf("File name for saving run results: ");
-	scanf("%s", outfilename);
+	if(scanf("%s", outfilename) != 1) return;
     }
 
     /* Initialize the method flags */
@@ -105,7 +105,7 @@ int      *ndims;		/* number of divisions at each stage */
 	while (*finassign == NULL) {
 	    if (PROMPT)
 		printf("Assignment input file: ");
-	    scanf("%s", inassignname);
+	    if(scanf("%s", inassignname) != 1) return;
 
 	    *finassign = fopen(inassignname, "r");
 	    if (*finassign == NULL) {
@@ -118,7 +118,7 @@ int      *ndims;		/* number of divisions at each stage */
 	while (*fingeom == NULL) {
 	    if (PROMPT)
 		printf("Geometry input file name: ");
-	    scanf("%s", geomname);
+	    if(scanf("%s", geomname) != 1) return;
 
 	    *fingeom = fopen(geomname, "r");
 	    if (*fingeom == NULL) {
@@ -141,7 +141,7 @@ int      *ndims;		/* number of divisions at each stage */
                 while (*fingeom == NULL) {
                     if (PROMPT)
                         printf("Geometry input file name: ");
-                    scanf("%s", geomname);
+                    if(scanf("%s", geomname) != 1) return;
  
                     *fingeom = fopen(geomname, "r");
                     if (*fingeom == NULL) {
@@ -163,7 +163,7 @@ int      *ndims;		/* number of divisions at each stage */
             while (*fingeom == NULL) {
                 if (PROMPT)
                     printf("Geometry input file name: ");
-                scanf("%s", geomname);
+                if(scanf("%s", geomname) != 1) return;
  
                 *fingeom = fopen(geomname, "r");
                 if (*fingeom == NULL) {
diff --git a/contrib/blossom/concorde97/TSP/bcontrol.c b/contrib/blossom/concorde97/TSP/bcontrol.c
index 6fce5b7edb878ba864fad713da07241a29190536..c268019b698884b2b6437c5a3b885c771712091c 100644
--- a/contrib/blossom/concorde97/TSP/bcontrol.c
+++ b/contrib/blossom/concorde97/TSP/bcontrol.c
@@ -642,16 +642,16 @@ CCtsp_lp *lp;
 
     printf ("Enter the (integer) id's for the two child nodes: ");
     fflush (stdout);
-    scanf ("%d %d", &ch0, &ch1);
+    if (scanf ("%d %d", &ch0, &ch1) != 2) { rval = 1; goto CLEANUP; }
 
     printf ("Enter 0 if edge-branch and 1 if clique-branch: ");
     fflush (stdout);
-    scanf ("%d", &tbran);
+    if (scanf ("%d", &tbran) != 1) { rval = 1; goto CLEANUP; }
 
     if (!tbran) {
         printf ("Enter ends of branching edge (use neg if original): ");
         fflush (stdout);
-        scanf ("%d %d", &bend0, &bend1);
+        if (scanf ("%d %d", &bend0, &bend1) != 2) { rval = 1; goto CLEANUP; }
         if (bend0 < 0) {
             if (bend1 >= 0) {
                 fprintf (stderr, "both ends must be from the same order\n");
@@ -670,7 +670,7 @@ CCtsp_lp *lp;
     } else {
         printf ("Enter the number of segments in clique: ");
         fflush (stdout);
-        scanf ("%d", &nseg);
+        if (scanf ("%d", &nseg) != 1) { rval = 1; goto CLEANUP; }
         slist = CC_SAFE_MALLOC (2*nseg, int);
         if (!slist) {
             fprintf (stderr, "out of memory\n");
@@ -679,7 +679,7 @@ CCtsp_lp *lp;
         printf ("Enter the ends of the segments: ");
         fflush (stdout);
         for (i = 0; i < nseg; i++) {
-            scanf ("%d %d", &slist[2*i], &slist[2*i+1]);
+            if (scanf ("%d %d", &slist[2*i], &slist[2*i+1]) != 2) { rval = 1; goto CLEANUP; }
         }
         c = CC_SAFE_MALLOC (1, CCtsp_lpclique);
         if (!c) {
diff --git a/contrib/blossom/concorde97/TSP/cutcall.c b/contrib/blossom/concorde97/TSP/cutcall.c
index 63f8b25f8a21bdda7f597fd3bd82b238e275d576..24d265a50da1e76f84d0b0af2fdd6db2672cfb0d 100644
--- a/contrib/blossom/concorde97/TSP/cutcall.c
+++ b/contrib/blossom/concorde97/TSP/cutcall.c
@@ -1375,17 +1375,23 @@ int *tour;
             fprintf (stderr, "out of memory in CCtsp_file_cuts\n");
             rval = 1; goto CLEANUP;
         }
-        fscanf (in, "%d", &nhandles);
+        if (fscanf (in, "%d", &nhandles) != 1) {
+            rval = 1; goto CLEANUP;
+        }
         c->handlecount = nhandles;
         for (i = 0; i < ncliques; i++) {
-            fscanf (in, "%d", &size);
+            if (fscanf (in, "%d", &size) != 1) {
+                rval = 1; goto CLEANUP;
+            }
             icliq = CC_SAFE_MALLOC (size, int);
             if (!icliq) {
                 fprintf (stderr, "out of memory in CCtsp_file_cuts\n");
                 rval = 1; goto CLEANUP;
             }
             for (j = 0; j < size; j++) {
-                fscanf (in, "%d", &k);
+                if (fscanf (in, "%d", &k) != 1) {
+                    rval = 1; goto CLEANUP;
+                }
                 icliq[j] = inv[k];
             }
             rval = CCtsp_array_to_lpclique (icliq, size, &(c->cliques[i]));
@@ -1395,7 +1401,9 @@ int *tour;
             }
             CC_FREE (icliq, int);
         }
-        fscanf (in, "%d", &(c->rhs));
+        if (fscanf (in, "%d", &(c->rhs)) != 1) {
+            rval = 1; goto CLEANUP;
+        }
         c->sense = 'G';
         c->branch = 0;
         c->next = *cuts;
diff --git a/contrib/blossom/concorde97/UTIL/getdata.c b/contrib/blossom/concorde97/UTIL/getdata.c
index cb98d098d926c2b0e5e2ae04b17ae2a464d360be..577f56fc2c175f0e7ce6fa53f8fdf31ac4d5fee0 100644
--- a/contrib/blossom/concorde97/UTIL/getdata.c
+++ b/contrib/blossom/concorde97/UTIL/getdata.c
@@ -270,42 +270,49 @@ CCdatagroup *dat;
                     return 1;
                 }
             } else {
-                FILE *datin = fopen (datname, "r");
-                if (datin == (FILE *) NULL) {
-                    perror (datname);
-                    fprintf (stderr, "Unable to open %s for input\n", datname);
-                    return 1;
+              FILE *datin = fopen (datname, "r");
+              if (datin == (FILE *) NULL) {
+                perror (datname);
+                fprintf (stderr, "Unable to open %s for input\n", datname);
+                return 1;
+              }
+              if (fscanf (datin, "%d", ncount) != 1) { fclose(datin); return 1; }
+              printf ("nnodes = %d\n", *ncount);
+              dat->x = CC_SAFE_MALLOC (*ncount, double);
+              if (!dat->x) {
+                fclose (datin);
+                return 1;
+              }
+              dat->y = CC_SAFE_MALLOC (*ncount, double);
+              if (!dat->y) {
+                fclose (datin);
+                CCutil_freedatagroup (*ncount, dat);
+                return 1;
+              }
+              if ((innorm & CC_NORM_SIZE_BITS) == CC_D3_NORM_SIZE) {
+                dat->z = CC_SAFE_MALLOC (*ncount, double);
+                if (!dat->z) {
+                  fclose (datin);
+                  CCutil_freedatagroup (*ncount, dat);
+                  return 1;
                 }
-                fscanf (datin, "%d", ncount);
-                printf ("nnodes = %d\n", *ncount);
-                dat->x = CC_SAFE_MALLOC (*ncount, double);
-                if (!dat->x) {
-                    fclose (datin);
+                for (i = 0; i < *ncount; i++) {
+                  if (fscanf (datin, "%lf %lf %lf", &(dat->x[i]), &(dat->y[i]), &(dat->z[i])) != 1) {
+                    fclose(datin);
+                    CCutil_freedatagroup (*ncount, dat);
                     return 1;
+                  }
                 }
-                dat->y = CC_SAFE_MALLOC (*ncount, double);
-                if (!dat->y) {
-                    fclose (datin);
+              } else {
+                for (i = 0; i < *ncount; i++) {
+                  if (fscanf (datin, "%lf %lf", &(dat->x[i]), &(dat->y[i])) != 1) {
+                    fclose(datin);
                     CCutil_freedatagroup (*ncount, dat);
                     return 1;
+                  }
                 }
-                if ((innorm & CC_NORM_SIZE_BITS) == CC_D3_NORM_SIZE) {
-                    dat->z = CC_SAFE_MALLOC (*ncount, double);
-                    if (!dat->z) {
-                        fclose (datin);
-                        CCutil_freedatagroup (*ncount, dat);
-                        return 1;
-                    }
-                    for (i = 0; i < *ncount; i++) {
-                        fscanf (datin, "%lf %lf %lf",
-                             &(dat->x[i]), &(dat->y[i]), &(dat->z[i]));
-                    }
-                } else {
-                    for (i = 0; i < *ncount; i++) {
-                        fscanf (datin, "%lf %lf", &(dat->x[i]), &(dat->y[i]));
-                    }
-                }
-                fclose (datin);
+              }
+              fclose (datin);
             }
         } else {
             printf ("Random %d point set\n", *ncount);
@@ -505,8 +512,8 @@ CCdatagroup *dat;
                 fprintf (stderr, "Unable to open %s for input\n", datname);
                 return 1;
             }
-            fscanf (datin, "%d", &seed);
-            fscanf (datin, "%d", &maxdist);
+            if (fscanf (datin, "%d", &seed) != 1) { fclose(datin); return 1; }
+            if (fscanf (datin, "%d", &maxdist) != 1) { fclose(datin); return 1; }
             fclose (datin);
         } else {
             seed = 1;
@@ -811,20 +818,28 @@ double **wcoord;
     if (weightname != (char *) NULL) {
         FILE *weightin = fopen (weightname, "r");
         if (weightin == (FILE *) NULL) {
-            perror (weightname);
-            fprintf (stderr, "Unable to open %s for input\n", weightname);
-            CC_FREE (*wcoord, double);
-            return 1;
+          perror (weightname);
+          fprintf (stderr, "Unable to open %s for input\n", weightname);
+          CC_FREE (*wcoord, double);
+          return 1;
+        }
+        if (fscanf (weightin, "%d", &k) != 1) {
+          fclose(weightin);
+          CC_FREE (*wcoord, double);
+          return 1;
         }
-        fscanf (weightin, "%d", &k);
         if (k != ncount) {
-            fprintf (stderr, "Weight file does not match node file\n");
-            fclose (weightin);
-            CC_FREE (*wcoord, double);
-            return 1;
+          fprintf (stderr, "Weight file does not match node file\n");
+          fclose (weightin);
+          CC_FREE (*wcoord, double);
+          return 1;
         }
         for (i = 0; i < k; i++) {
-            fscanf (weightin, "%lf", &((*wcoord)[i]));
+          if (fscanf (weightin, "%lf", &((*wcoord)[i])) != 1) {
+            fclose(weightin);
+            CC_FREE (*wcoord, double);
+            return 1;
+          }
         }
         make_weights_nonnegative (ncount, *wcoord);
         fclose (weightin);
@@ -991,7 +1006,10 @@ CCdatagroup *dat;
                         return 1;
                     }
                     for (i = 0; i < *ncount; i++) {
-                        fscanf (in, "%*d %lf %lf", &(dat->x[i]), &(dat->y[i]));
+                      if (fscanf (in, "%*d %lf %lf", &(dat->x[i]), &(dat->y[i])) != 2) {
+                        CCutil_freedatagroup (*ncount, dat);
+                        return 1;
+                      }
                     }
                 } else if (((dat->norm) & CC_NORM_SIZE_BITS) ==
                                             CC_D3_NORM_SIZE) {
@@ -1011,8 +1029,11 @@ CCdatagroup *dat;
                         return 1;
                     }
                     for (i = 0; i < *ncount; i++) {
-                        fscanf (in, "%*d %lf %lf %lf",
-                               &(dat->x[i]), &(dat->y[i]), &(dat->z[i]));
+                      if (fscanf (in, "%*d %lf %lf %lf",
+                               &(dat->x[i]), &(dat->y[i]), &(dat->z[i])) != 3) {
+                        CCutil_freedatagroup (*ncount, dat);
+                        return 1;
+                      }
                     }
                 } else {
                     fprintf (stderr, "ERROR: Node coordinates with norm %d?\n",
@@ -1043,8 +1064,12 @@ CCdatagroup *dat;
                                 CCutil_freedatagroup (*ncount, dat);
                                 return 1;
                             }
-                            for (j = 0; j <= i; j++)
-                                fscanf (in, "%d", &(dat->adj[i][j]));
+                            for (j = 0; j <= i; j++){
+                              if (fscanf (in, "%d", &(dat->adj[i][j])) != 1) {
+                                CCutil_freedatagroup (*ncount, dat);
+                                return 1;
+                              }
+                            }
                         }
                     } else if (matrixform == MATRIX_UPPER_ROW ||
                                matrixform == MATRIX_UPPER_DIAG_ROW ||
@@ -1066,14 +1091,33 @@ CCdatagroup *dat;
                             }
                             if (matrixform == MATRIX_UPPER_ROW) {
                                 tempadj[i][i] = 0;
-                                for (j = i + 1; j < *ncount; j++)
-                                    fscanf (in, "%d", &(tempadj[i][j]));
+                                for (j = i + 1; j < *ncount; j++) {
+                                    if (fscanf (in, "%d", &(tempadj[i][j])) != 1) {
+                                        CCutil_freedatagroup (*ncount, dat);
+                                        for (j = 0; j < i; j++)
+                                            CC_FREE (tempadj[j], int);
+                                        CC_FREE (tempadj, int *);
+                                        return 1;
+                                    }
+                                }
                             } else if (matrixform == MATRIX_UPPER_DIAG_ROW) {
-                                for (j = i; j < *ncount; j++)
-                                    fscanf (in, "%d", &(tempadj[i][j]));
+                                for (j = i; j < *ncount; j++) {
+                                    if (fscanf (in, "%d", &(tempadj[i][j])) != 1) {
+                                        CCutil_freedatagroup (*ncount, dat);
+                                        for (j = 0; j < i; j++)
+                                            CC_FREE (tempadj[j], int);
+                                        CC_FREE (tempadj, int *);
+                                    }
+                                }
                             } else {
-                                for (j = 0; j < *ncount; j++)
-                                    fscanf (in, "%d", &(tempadj[i][j]));
+                                for (j = 0; j < *ncount; j++) {
+                                    if (fscanf (in, "%d", &(tempadj[i][j])) != 1) {
+                                        CCutil_freedatagroup (*ncount, dat);
+                                        for (j = 0; j < i; j++)
+                                            CC_FREE (tempadj[j], int);
+                                        CC_FREE (tempadj, int *);
+                                    }
+                                }
                             }
                         }
                         for (i = 0; i < *ncount; i++) {
@@ -1291,7 +1335,11 @@ int *outcycle;
         return 1;
     }
 
-    fscanf (cycin, "%d %d", &i, &k);
+    if (fscanf (cycin, "%d %d", &i, &k) != 2) {
+        CC_FREE (elist, int);
+        fclose (cycin);
+        return 1;
+    }
     if (i != ncount || k != ncount) {
         fprintf (stderr, "file is not a cycle-edge file for this problem\n");
         CC_FREE (elist, int);
@@ -1299,8 +1347,13 @@ int *outcycle;
         return 1;
     }
 
-    for (i = 0; i < ncount; i++)
-        fscanf (cycin, "%d %d %*d", &(elist[2 * i]), &(elist[(2 * i) + 1]));
+    for (i = 0; i < ncount; i++) {
+        if (fscanf (cycin, "%d %d %*d", &(elist[2 * i]), &(elist[(2 * i) + 1])) != 2) {
+            CC_FREE (elist, int);
+            fclose (cycin);
+            return 1;
+        }
+    }
 
     if (CCutil_edge_to_cycle (ncount, elist, outcycle)) {
         fprintf (stderr, "CCutil_edge_to_cycle failed\n");
@@ -1720,13 +1773,22 @@ CCdatagroup *dat;
 
         for (i = 0; i < 3; i++) {
             for (j = 0; j < 3; j++) {
-                fscanf (datin, "%lf", &orient[i][j]);
+                if (fscanf (datin, "%lf", &orient[i][j]) != 1) {
+                    fclose (datin);
+                    return 1;
+                }
             }
         }
-        fscanf (datin, "%lf", &lambda);
+        if (fscanf (datin, "%lf", &lambda) != 1) {
+            fclose (datin);
+            return 1;
+        }
         for (i = 0; i < 3; i++) {
             for (j = 0; j < 2; j++) {
-                fscanf (datin, "%d", &(bounds[i][j]));
+                if (fscanf (datin, "%d", &(bounds[i][j])) != 1) {
+                    fclose (datin);
+                    return 1;
+                }
             }
         }
         fclose (datin);
diff --git a/contrib/kbipack/gmp_matrix_io.cpp b/contrib/kbipack/gmp_matrix_io.cpp
index 14b0fe36501793ff3736e0424ffdb83230aa80b9..2300f9437bd24435f6e81c8adeabef35c44d49ba 100644
--- a/contrib/kbipack/gmp_matrix_io.cpp
+++ b/contrib/kbipack/gmp_matrix_io.cpp
@@ -44,10 +44,10 @@ gmp_matrix * gmp_matrix_read_coord(char* filename)
 
   /* Matlab and Octave typically include comments on ascii files. */
   /* They start with # */
-  fgets(buffer, 999, p_file);
+  if(!fgets(buffer, 999, p_file)) return NULL;
   while(strncmp(buffer, "#",1) == 0)
     {
-      fgets(buffer, 999, p_file);
+      if(!fgets(buffer, 999, p_file)) return NULL;
     }
 
   /* First read the size of the matrix */
diff --git a/contrib/mmg3d/build/sources/inout.c b/contrib/mmg3d/build/sources/inout.c
index 09e719fab4e0d43a4d31e8adfbe1f75a0366e922..5f73da0c925f5bdbc3d9d7680a84608455747063 100644
--- a/contrib/mmg3d/build/sources/inout.c
+++ b/contrib/mmg3d/build/sources/inout.c
@@ -3,32 +3,32 @@ Logiciel initial: MMG3D Version 4.0
 Co-auteurs : Cecile Dobrzynski et Pascal Frey.
 Propriétaires :IPB - UPMC -INRIA.
 
-Copyright © 2004-2005-2006-2007-2008-2009-2010-2011, 
+Copyright © 2004-2005-2006-2007-2008-2009-2010-2011,
 diffusé sous les termes et conditions de la licence publique générale de GNU
-Version 3 ou toute version ultérieure.  
+Version 3 ou toute version ultérieure.
 
 Ce fichier est une partie de MMG3D.
 MMG3D est un logiciel libre ; vous pouvez le redistribuer et/ou le modifier
 suivant les termes de la licence publique générale de GNU
 Version 3 ou toute version ultérieure.
-MMG3D est distribué dans l'espoir qu'il sera utile, mais SANS 
-AUCUNE GARANTIE ; sans même garantie de valeur marchande.  
+MMG3D est distribué dans l'espoir qu'il sera utile, mais SANS
+AUCUNE GARANTIE ; sans même garantie de valeur marchande.
 Voir la licence publique générale de GNU pour plus de détails.
-MMG3D est diffusé en espérant qu’il sera utile, 
-mais SANS AUCUNE GARANTIE, ni explicite ni implicite, 
-y compris les garanties de commercialisation ou 
-d’adaptation dans un but spécifique. 
+MMG3D est diffusé en espérant qu’il sera utile,
+mais SANS AUCUNE GARANTIE, ni explicite ni implicite,
+y compris les garanties de commercialisation ou
+d’adaptation dans un but spécifique.
 Reportez-vous à la licence publique générale de GNU pour plus de détails.
-Vous devez avoir reçu une copie de la licence publique générale de GNU 
-en même temps que ce document. 
+Vous devez avoir reçu une copie de la licence publique générale de GNU
+en même temps que ce document.
 Si ce n’est pas le cas, aller voir <http://www.gnu.org/licenses/>.
 /****************************************************************************
 Initial software: MMG3D Version 4.0
 Co-authors: Cecile Dobrzynski et Pascal Frey.
 Owners: IPB - UPMC -INRIA.
 
-Copyright © 2004-2005-2006-2007-2008-2009-2010-2011, 
-spread under the terms and conditions of the license GNU General Public License 
+Copyright © 2004-2005-2006-2007-2008-2009-2010-2011,
+spread under the terms and conditions of the license GNU General Public License
 as published Version 3, or (at your option) any later version.
 
 This file is part of MMG3D
@@ -41,7 +41,7 @@ but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 You should have received a copy of the GNU General Public License
-along with MMG3D. If not, see <http://www.gnu.org/licenses/>.  
+along with MMG3D. If not, see <http://www.gnu.org/licenses/>.
 ****************************************************************************/
 #include "mesh.h"
 
@@ -52,16 +52,16 @@ extern short	       MMG_imprim;
 
 int MMG_swapbin(int sbin)
 {
-	int inv; 
+	int inv;
   char *p_in = (char *) &sbin;
   char *p = (char *)&inv;
 
- 
+
   p[0] = p_in[3];
   p[1] = p_in[2];
   p[2] = p_in[1];
-  p[3] = p_in[0]; 
-  
+  p[3] = p_in[0];
+
   return(inv);
   /*unsigned char c1, c2, c3, c4;
 
@@ -69,62 +69,62 @@ int MMG_swapbin(int sbin)
   c2 = (sbin >> 8) & 255;
   c3 = (sbin >> 16) & 255;
   c4 = (sbin >> 24) & 255;
-  
+
   return ((int)c1 << 24) + ((int)c2 << 16) + ((int)c3 << 8) + c4;   */
 
 }
 float MMG_swapf(float sbin)
-{ 
+{
   float out;
-  char *p_in = (char *) &sbin; 
-  char *p_out = (char *) &out; 
-  p_out[0] = p_in[3]; 
-  p_out[1] = p_in[2]; 
-  p_out[2] = p_in[1]; 
+  char *p_in = (char *) &sbin;
+  char *p_out = (char *) &out;
+  p_out[0] = p_in[3];
+  p_out[1] = p_in[2];
+  p_out[2] = p_in[1];
   p_out[3] = p_in[0];
-  
+
   return(out);
 }
 double MMG_swapd(double sbin)
 {
   float out;
-  char *p_in = (char *) &sbin; 
-  char *p_out = (char *) &out; 
+  char *p_in = (char *) &sbin;
+  char *p_out = (char *) &out;
   int i;
-  
+
   for(i=0;i<8;i++)
   {
     p_out[i] = p_in[7-i];
-  }          
+  }
   //printf("CONVERTION DOUBLE\n");
   return(out);
 }
 
 /* read mesh data */
-int MMG_loadMesh(pMesh mesh,char *filename) {  
+int MMG_loadMesh(pMesh mesh,char *filename) {
   FILE*            inm;
-  Hedge    				 hed,hed2;
+  Hedge            hed,hed2;
   pPoint       	   ppt;
   pTetra           pt;
   pTria            pt1;
-  int              k,dim,ref,bin,bpos,i,tmp;   
+  int              k,dim,ref,bin,bpos,i,tmp;
   long             posnp,posnt,posne,posnhex,posnpris,posncor,posned,posnq;
   char            *ptr,data[128],chaine[128];
   int              nhex,npris,netmp,ncor,ned,nq;
-  int              p0,p1,p2,p3,p4,p5,p6,p7;  
-  int              binch,bdim,iswp,nu1,nu2,nimp,ne;       
+  int              p0,p1,p2,p3,p4,p5,p6,p7;
+  int              binch,bdim,iswp,nu1,nu2,nimp,ne;
   float            fc;
-  
+
 
   posnp = posnt = posne = posnhex = posnpris = 0;
   netmp = ncor = ned = 0;
   bin = 0;
   iswp = 0;
-  mesh->np = mesh->nt = mesh->ne = mesh->ncor = 0; 
+  mesh->np = mesh->nt = mesh->ne = mesh->ncor = 0;
   npris = nhex = nq = 0;
-  
+
   strcpy(data,filename);
-  ptr = strstr(data,".mesh");  
+  ptr = strstr(data,".mesh");
   if ( !ptr ) {
     strcat(data,".meshb");
     if( !(inm = fopen(data,"rb")) ) {
@@ -145,88 +145,88 @@ int MMG_loadMesh(pMesh mesh,char *filename) {
       if( !(inm = fopen(data,"r")) ) {
         fprintf(stderr,"  ** %s  NOT FOUND.\n",data);
         return(0);
-      }      
+      }
     } else {
       bin = 1;
       if( !(inm = fopen(data,"rb")) ) {
         fprintf(stderr,"  ** %s  NOT FOUND.\n",data);
         return(0);
       }
-      
-    }  
+
+    }
   }
 
   fprintf(stdout,"  %%%% %s OPENED\n",data);
   if (!bin) {
-    strcpy(chaine,"D");     
-    while(fscanf(inm,"%s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) { 
+    strcpy(chaine,"D");
+    while(fscanf(inm,"%s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
       if(!strncmp(chaine,"MeshVersionFormatted",strlen("MeshVersionFormatted"))) {
-          fscanf(inm,"%d",&mesh->ver);
-          continue;
+        if(fscanf(inm,"%d",&mesh->ver)!=1){ fclose(inm); return(0); }
+        continue;
       } else if(!strncmp(chaine,"Dimension",strlen("Dimension"))) {
-          fscanf(inm,"%d",&dim);
-          if(dim!=3) {
-            fprintf(stdout,"BAD DIMENSION : %d\n",dim);
-            return(0);
-          }
-          continue;
+        if(fscanf(inm,"%d",&dim)!=1){ fclose(inm); return(0); }
+        if(dim!=3) {
+          fprintf(stdout,"BAD DIMENSION : %d\n",dim);
+          return(0);
+        }
+        continue;
       } else if(!strncmp(chaine,"Vertices",strlen("Vertices"))) {
-        fscanf(inm,"%d",&mesh->np); 
+        if(fscanf(inm,"%d",&mesh->np)!=1){ fclose(inm); return(0); }
         posnp = ftell(inm);
         continue;
       } else if(!strncmp(chaine,"Triangles",strlen("Triangles"))) {
-        fscanf(inm,"%d",&mesh->nt); 
+        if(fscanf(inm,"%d",&mesh->nt)!=1){ fclose(inm); return(0); }
         posnt = ftell(inm);
         continue;
       } else if(!strncmp(chaine,"Tetrahedra",strlen("Tetrahedra"))) {
-        fscanf(inm,"%d",&mesh->ne); 
+        if(fscanf(inm,"%d",&mesh->ne)!=1){ fclose(inm); return(0); }
         netmp = mesh->ne;
         posne = ftell(inm);
         continue;
-      } else if(!strncmp(chaine,"Hexahedra",strlen("Hexahedra"))) { 
-        assert(abs(mesh->info.option)==10);  
-        fscanf(inm,"%d",&nhex); 
+      } else if(!strncmp(chaine,"Hexahedra",strlen("Hexahedra"))) {
+        assert(abs(mesh->info.option)==10);
+        if(fscanf(inm,"%d",&nhex)!=1){ fclose(inm); return(0); }
 				//nhex=0;
         posnhex = ftell(inm);
         continue;
-      } else if(!strncmp(chaine,"Pentahedra",strlen("Pentahedra"))) { 
-        assert(abs(mesh->info.option)==10); 
-        fscanf(inm,"%d",&npris);
+      } else if(!strncmp(chaine,"Pentahedra",strlen("Pentahedra"))) {
+        assert(abs(mesh->info.option)==10);
+        if(fscanf(inm,"%d",&npris)!=1){ fclose(inm); return(0); }
 				//npris=0;
         posnpris = ftell(inm);
         continue;
-      } else if(!strncmp(chaine,"Corners",strlen("Corners"))) { 
-        fscanf(inm,"%d",&ncor); 
+      } else if(!strncmp(chaine,"Corners",strlen("Corners"))) {
+        if(fscanf(inm,"%d",&ncor)!=1){ fclose(inm); return(0); }
         posncor = ftell(inm);
         continue;
-      } else if(!strncmp(chaine,"Edges",strlen("Edges"))) { 
-	      fscanf(inm,"%d",&ned); 
-	      posned = ftell(inm);
-	      continue;
-	    } else if(abs(mesh->info.option)==10 && !strncmp(chaine,"Quadrilaterals",strlen("Quadrilaterals"))) {
-		    fscanf(inm,"%d",&nq); 
-		    posnq = ftell(inm);
-		    continue;
-		  }
-    }  
+      } else if(!strncmp(chaine,"Edges",strlen("Edges"))) {
+        if(fscanf(inm,"%d",&ned)!=1){ fclose(inm); return(0); }
+        posned = ftell(inm);
+        continue;
+      } else if(abs(mesh->info.option)==10 && !strncmp(chaine,"Quadrilaterals",strlen("Quadrilaterals"))) {
+        if(fscanf(inm,"%d",&nq)!=1){ fclose(inm); return(0); }
+        posnq = ftell(inm);
+        continue;
+      }
+    }
   } else {
     bdim = 0;
-    fread(&mesh->ver,sw,1,inm);
-    iswp=0;   
-    if(mesh->ver==16777216) 
-      iswp=1;    
+    if(fread(&mesh->ver,sw,1,inm)!=1){ fclose(inm); return(0); }
+    iswp=0;
+    if(mesh->ver==16777216)
+      iswp=1;
     else if(mesh->ver!=1) {
       fprintf(stdout,"BAD FILE ENCODING\n");
-    } 
-    fread(&mesh->ver,sw,1,inm); 
-    if(iswp) mesh->ver = MMG_swapbin(mesh->ver); 
-    while(fread(&binch,sw,1,inm)!=0 && binch!=54 ) {  
-      if(iswp) binch=MMG_swapbin(binch);      
-      if(binch==54) break;  
+    }
+    if(fread(&mesh->ver,sw,1,inm)!=1){ fclose(inm); return(0); }
+    if(iswp) mesh->ver = MMG_swapbin(mesh->ver);
+    while(fread(&binch,sw,1,inm)!=0 && binch!=54 ) {
+      if(iswp) binch=MMG_swapbin(binch);
+      if(binch==54) break;
       if(!bdim && binch==3) {  //Dimension
-        fread(&bdim,sw,1,inm);  //NulPos=>20
+        if(fread(&bdim,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos=>20
         if(iswp) bdim=MMG_swapbin(bdim);
-        fread(&bdim,sw,1,inm);  
+        if(fread(&bdim,sw,1,inm)!=1){ fclose(inm); return(0); }
         if(iswp) bdim=MMG_swapbin(bdim);
         if(bdim!=3) {
           fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim);
@@ -235,87 +235,87 @@ int MMG_loadMesh(pMesh mesh,char *filename) {
         }
         continue;
       } else if(!mesh->np && binch==4) {  //Vertices
-        fread(&bpos,sw,1,inm); //NulPos 
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
         if(iswp) bpos=MMG_swapbin(bpos);
-        fread(&mesh->np,sw,1,inm); 
+        if(fread(&mesh->np,sw,1,inm)!=1){ fclose(inm); return(0); }
         if(iswp) mesh->np=MMG_swapbin(mesh->np);
-        posnp = ftell(inm);     
+        posnp = ftell(inm);
         rewind(inm);
-        fseek(inm,bpos,SEEK_SET);        
+        fseek(inm,bpos,SEEK_SET);
         continue;
       }  else if(!mesh->nt && binch==6) {//Triangles
-        fread(&bpos,sw,1,inm); //NulPos 
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
         if(iswp) bpos=MMG_swapbin(bpos);
-        fread(&mesh->nt,sw,1,inm); 
+        if(fread(&mesh->nt,sw,1,inm)!=1){ fclose(inm); return(0); }
         if(iswp) mesh->nt=MMG_swapbin(mesh->nt);
-        posnt = ftell(inm); 
+        posnt = ftell(inm);
+        rewind(inm);
+        fseek(inm,bpos,SEEK_SET);
+        continue;
+      } else if(!mesh->ne && binch==8) {
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
+        if(iswp) bpos=MMG_swapbin(bpos);
+        if(fread(&mesh->ne,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) mesh->ne=MMG_swapbin(mesh->ne);
+        netmp = mesh->ne;
+        posne = ftell(inm);
+        rewind(inm);
+        fseek(inm,bpos,SEEK_SET);
+        continue;
+      } else if(!nhex && binch==10) {
+        assert(abs(mesh->info.option)==10);
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
+        if(iswp) bpos=MMG_swapbin(bpos);
+        if(fread(&nhex,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) nhex=MMG_swapbin(nhex);
+        posnhex = ftell(inm);
+        rewind(inm);
+        fseek(inm,bpos,SEEK_SET);
+        continue;
+      } else if(!npris && binch==9) {
+        assert(abs(mesh->info.option)==10);
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
+        if(iswp) bpos=MMG_swapbin(bpos);
+        if(fread(&npris,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) npris=MMG_swapbin(npris);
+        posnpris = ftell(inm);
+        rewind(inm);
+        fseek(inm,bpos,SEEK_SET);
+        continue;
+      } else if(!ncor && binch==13) {
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
+        if(iswp) bpos=MMG_swapbin(bpos);
+        if(fread(&ncor,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ncor=MMG_swapbin(ncor);
+        posncor = ftell(inm);
+        rewind(inm);
+        fseek(inm,bpos,SEEK_SET);
+        continue;
+      } else if(!ned && binch==5) { //Edges
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
+        if(iswp) bpos=MMG_swapbin(bpos);
+        if(fread(&ned,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ned=MMG_swapbin(ned);
+        posned = ftell(inm);
         rewind(inm);
-        fseek(inm,bpos,SEEK_SET);        
+        fseek(inm,bpos,SEEK_SET);
         continue;
-       } else if(!mesh->ne && binch==8) {  
-         fread(&bpos,sw,1,inm); //NulPos 
-         if(iswp) bpos=MMG_swapbin(bpos);
-         fread(&mesh->ne,sw,1,inm); 
-         if(iswp) mesh->ne=MMG_swapbin(mesh->ne);
-         netmp = mesh->ne;
-         posne = ftell(inm);
-         rewind(inm);
-         fseek(inm,bpos,SEEK_SET);        
-         continue;
-       } else if(!nhex && binch==10) { 
-         assert(abs(mesh->info.option)==10);
-         fread(&bpos,sw,1,inm); //NulPos 
-         if(iswp) bpos=MMG_swapbin(bpos);
-         fread(&nhex,sw,1,inm); 
-         if(iswp) nhex=MMG_swapbin(nhex);
-         posnhex = ftell(inm);
-         rewind(inm);
-         fseek(inm,bpos,SEEK_SET);        
-         continue;
-       } else if(!npris && binch==9) { 
-         assert(abs(mesh->info.option)==10);
-         fread(&bpos,sw,1,inm); //NulPos 
-         if(iswp) bpos=MMG_swapbin(bpos);
-         fread(&npris,sw,1,inm); 
-         if(iswp) npris=MMG_swapbin(npris);
-         posnpris = ftell(inm);
-         rewind(inm);
-         fseek(inm,bpos,SEEK_SET);        
-         continue;
-       } else if(!ncor && binch==13) { 
-         fread(&bpos,sw,1,inm); //NulPos 
-         if(iswp) bpos=MMG_swapbin(bpos);
-         fread(&ncor,sw,1,inm);          
-         if(iswp) ncor=MMG_swapbin(ncor);
-         posncor = ftell(inm);
-         rewind(inm);
-         fseek(inm,bpos,SEEK_SET);        
-         continue;
-        } else if(!ned && binch==5) { //Edges
-	       fread(&bpos,sw,1,inm); //NulPos 
-	       if(iswp) bpos=MMG_swapbin(bpos);
-	       fread(&ned,sw,1,inm); 
-	       if(iswp) ned=MMG_swapbin(ned);
-	       posned = ftell(inm);
-	       rewind(inm);
-	       fseek(inm,bpos,SEEK_SET);        
-	       continue;
-	      } else {
-         //printf("on traite ? %d\n",binch);
-         fread(&bpos,sw,1,inm); //NulPos 
-         if(iswp) bpos=MMG_swapbin(bpos);
-         //printf("on avance... Nulpos %d\n",bpos);         
-         rewind(inm);
-         fseek(inm,bpos,SEEK_SET);        
-       }     
-    }            
-    
+      } else {
+        //printf("on traite ? %d\n",binch);
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
+        if(iswp) bpos=MMG_swapbin(bpos);
+        //printf("on avance... Nulpos %d\n",bpos);
+        rewind(inm);
+        fseek(inm,bpos,SEEK_SET);
+      }
+    }
+
   }
 
   if ( abs(mesh->info.option)==10 ) {
-    fprintf(stdout,"  -- READING %8d HEXA %8d PRISMS\n",nhex,npris);  
-    if(!mesh->ne) netmp = 0;  
-    mesh->ne += 6*nhex + 3*npris; 
+    fprintf(stdout,"  -- READING %8d HEXA %8d PRISMS\n",nhex,npris);
+    if(!mesh->ne) netmp = 0;
+    mesh->ne += 6*nhex + 3*npris;
   }
 
   if ( abs(mesh->info.imprim) > 5 )
@@ -333,35 +333,35 @@ int MMG_loadMesh(pMesh mesh,char *filename) {
   fseek(inm,posnp,SEEK_SET);
   for (k=1; k<=mesh->np; k++) {
     ppt = &mesh->point[k];
-    if (mesh->ver < 2) { /*float*/ 
+    if (mesh->ver < 2) { /*float*/
       if (!bin) {
         for (i=0 ; i<3 ; i++) {
-          fscanf(inm,"%f",&fc);
+          if(fscanf(inm,"%f",&fc)!=1){ fclose(inm); return(0); }
           ppt->c[i] = (double) fc;
-        } 
-        fscanf(inm,"%d",&ppt->ref);
+        }
+        if(fscanf(inm,"%d",&ppt->ref)!=1){ fclose(inm); return(0); }
       } else {
         for (i=0 ; i<3 ; i++) {
-          fread(&fc,sw,1,inm); 
-          if(iswp) fc=MMG_swapf(fc);    
+          if(fread(&fc,sw,1,inm)!=1){ fclose(inm); return(0); }
+          if(iswp) fc=MMG_swapf(fc);
           ppt->c[i] = (double) fc;
-        }     
-        fread(&ppt->ref,sw,1,inm);         
-        if(iswp) ppt->ref=MMG_swapbin(ppt->ref);    
+        }
+        if(fread(&ppt->ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ppt->ref=MMG_swapbin(ppt->ref);
       }
     } else {
-      if (!bin) 
-        fscanf(inm,"%lf %lf %lf %d",&ppt->c[0],&ppt->c[1],&ppt->c[2],&ppt->ref); 
+      if (!bin)
+        if(fscanf(inm,"%lf %lf %lf %d",&ppt->c[0],&ppt->c[1],&ppt->c[2],&ppt->ref)!=4){ fclose(inm); return(0); }
       else {
-        for (i=0 ; i<3 ; i++) { 
-          fread(&ppt->c[i],sd,1,inm);
-          if(iswp) ppt->c[i]=MMG_swapd(ppt->c[i]); 
-        }   
-        fread(&ppt->ref,sw,1,inm);         
-        if(iswp) ppt->ref=MMG_swapbin(ppt->ref);    
-      }  
-    }             
-    ppt->tag  = M_UNUSED;    
+        for (i=0 ; i<3 ; i++) {
+          if(fread(&ppt->c[i],sd,1,inm)!=1){ fclose(inm); return(0); }
+          if(iswp) ppt->c[i]=MMG_swapd(ppt->c[i]);
+        }
+        if(fread(&ppt->ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ppt->ref=MMG_swapbin(ppt->ref);
+      }
+    }
+    ppt->tag  = M_UNUSED;
   }
 
   /* read mesh triangles */
@@ -369,83 +369,83 @@ int MMG_loadMesh(pMesh mesh,char *filename) {
   rewind(inm);
   fseek(inm,posnt,SEEK_SET);
   for (k=1; k<=mesh->nt; k++) {
-    pt1 = &mesh->tria[k]; 
+    pt1 = &mesh->tria[k];
     if (!bin)
-      fscanf(inm,"%d %d %d %d",&pt1->v[0],&pt1->v[1],&pt1->v[2],&pt1->ref);
+      if(fscanf(inm,"%d %d %d %d",&pt1->v[0],&pt1->v[1],&pt1->v[2],&pt1->ref)!=4){ fclose(inm); return(0); }
     else {
-      for (i=0 ; i<3 ; i++) { 
-        fread(&pt1->v[i],sw,1,inm); 
-        if(iswp) pt1->v[i]=MMG_swapbin(pt1->v[i]);    
+      for (i=0 ; i<3 ; i++) {
+        if(fread(&pt1->v[i],sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) pt1->v[i]=MMG_swapbin(pt1->v[i]);
       }
-      fread(&pt1->ref,sw,1,inm); 
-      if(iswp) pt1->ref=MMG_swapbin(pt1->ref);           
-    }  
+      if(fread(&pt1->ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+      if(iswp) pt1->ref=MMG_swapbin(pt1->ref);
+    }
   }
-  /* read mesh quads (option 10)*/ 
-	if(abs(mesh->info.option)==10) { 
+  /* read mesh quads (option 10)*/
+	if(abs(mesh->info.option)==10) {
 		fprintf(stdout,"     QUADS READING %d\n",nq);
     mesh->ntfixe += 4*nq;
     rewind(inm);
     fseek(inm,posnq,SEEK_SET);
     for (k=1; k<=nq; k++) {
       if (!bin)
-        fscanf(inm,"%d %d %d %d %d",&p0,&p1,&p2,&p3,&ref);
+        if(fscanf(inm,"%d %d %d %d %d",&p0,&p1,&p2,&p3,&ref)!=5){ fclose(inm); return(0); }
       else {
-        fread(&p0,sw,1,inm); 
-        if(iswp) p0=MMG_swapbin(p0);    
-        fread(&p1,sw,1,inm); 
-        if(iswp) p1=MMG_swapbin(p1);    
-        fread(&p2,sw,1,inm); 
-        if(iswp) p2=MMG_swapbin(p2);    
-        fread(&p3,sw,1,inm); 
-        if(iswp) p3=MMG_swapbin(p3);    
-	      fread(&pt1->ref,sw,1,inm); 
-	      if(iswp) ref=MMG_swapbin(ref);           
-      } 
-      pt1 = &mesh->tria[++mesh->nt]; 
+        if(fread(&p0,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) p0=MMG_swapbin(p0);
+        if(fread(&p1,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) p1=MMG_swapbin(p1);
+        if(fread(&p2,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) p2=MMG_swapbin(p2);
+        if(fread(&p3,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) p3=MMG_swapbin(p3);
+        if(fread(&pt1->ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ref=MMG_swapbin(ref);
+      }
+      pt1 = &mesh->tria[++mesh->nt];
 			pt1->v[0] = p0;
 			pt1->v[1] = p1;
 			pt1->v[2] = p2;
 			pt1->ref  = ref;
-      pt1 = &mesh->tria[++mesh->nt]; 
+      pt1 = &mesh->tria[++mesh->nt];
 			pt1->v[0] = p0;
 			pt1->v[1] = p2;
 			pt1->v[2] = p3;
 			pt1->ref  = ref;
-      pt1 = &mesh->tria[++mesh->nt]; 
+      pt1 = &mesh->tria[++mesh->nt];
 			pt1->v[0] = p0;
 			pt1->v[1] = p1;
 			pt1->v[2] = p3;
 			pt1->ref  = ref;
-      pt1 = &mesh->tria[++mesh->nt]; 
+      pt1 = &mesh->tria[++mesh->nt];
 			pt1->v[0] = p1;
 			pt1->v[1] = p2;
 			pt1->v[2] = p3;
 			pt1->ref  = ref;
- 
+
     }
   }
 
 	/*read and store edges*/
-  if (ned) {         
+  if (ned) {
 	  if ( !MMG_zaldy4(&hed,ned) ) {
-      if ( mesh->info.ddebug )  fprintf(stdout,"  ## MEMORY ALLOCATION PROBLEM : EDGES IGNORED\n"); 
+      if ( mesh->info.ddebug )  fprintf(stdout,"  ## MEMORY ALLOCATION PROBLEM : EDGES IGNORED\n");
 			ned = 0;
-    }   
+    }
 		mesh->ned = ned;
     rewind(inm);
-    fseek(inm,posned,SEEK_SET); 
-    for (k=1; k<=ned; k++) { 
+    fseek(inm,posned,SEEK_SET);
+    for (k=1; k<=ned; k++) {
       if (!bin)
-        fscanf(inm,"%d %d %d",&nu1,&nu2,&ref);
+        if(fscanf(inm,"%d %d %d",&nu1,&nu2,&ref)!=3){ fclose(inm); return(0); }
       else {
-        fread(&nu1,sw,1,inm); 
-        if(iswp) nu1=MMG_swapbin(nu1);    
-        fread(&nu2,sw,1,inm); 
-        if(iswp) nu2=MMG_swapbin(nu2);    
-        fread(&ref,sw,1,inm); 
-        if(iswp) ref=MMG_swapbin(ref);    
-      }  
+        if(fread(&nu1,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) nu1=MMG_swapbin(nu1);
+        if(fread(&nu2,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) nu2=MMG_swapbin(nu2);
+        if(fread(&ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ref=MMG_swapbin(ref);
+      }
 			if(MMG_edgePut(&hed,nu1,nu2,2)>1) {
 				fprintf(stdout,"  ## WARNING DOUBLE EDGE : %d %d\n",nu1,nu2);
 			}
@@ -456,97 +456,96 @@ int MMG_loadMesh(pMesh mesh,char *filename) {
   mesh->nefixe = mesh->ne;
   rewind(inm);
   fseek(inm,posne,SEEK_SET);
-  for (k=1; k<=netmp; k++) { 
+  for (k=1; k<=netmp; k++) {
     pt = &mesh->tetra[k];
-    if (!bin) 
-      fscanf(inm,"%d %d %d %d %d",&pt->v[0],&pt->v[1],&pt->v[2],&pt->v[3],&ref); 
-    else {                                                                        
-	
-      for (i=0 ; i<4 ; i++) { 
-        fread(&pt->v[i],sw,1,inm); 
-        if(iswp) pt->v[i]=MMG_swapbin(pt->v[i]);    
+    if (!bin)
+      if(fscanf(inm,"%d %d %d %d %d",&pt->v[0],&pt->v[1],&pt->v[2],&pt->v[3],&ref)!=5){ fclose(inm); return(0); }
+    else {
+
+      for (i=0 ; i<4 ; i++) {
+        if(fread(&pt->v[i],sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) pt->v[i]=MMG_swapbin(pt->v[i]);
       }
-      fread(&ref,sw,1,inm);         
-      if(iswp) ref=MMG_swapbin(ref);           
-    }  
-    pt->ref  = ref;//0;//ref ;  
+      if(fread(&ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+      if(iswp) ref=MMG_swapbin(ref);
+    }
+    pt->ref  = ref;//0;//ref ;
     for(i=0 ; i<4 ; i++)
-      pt->bdryref[i] = -1;  
-    
+      pt->bdryref[i] = -1;
+
 		if (ned) {
-	    for(i=0 ; i<6 ; i++) {                         
+	    for(i=0 ; i<6 ; i++) {
 				nu1 = pt->v[MMG_iare[i][0]];
 				nu2 = pt->v[MMG_iare[i][1]];
 	      pt->bdryinfo[i] = MMG_edgePoint(&hed,nu1,nu2);
-			}  			
-			
+			}
+
 		} else {
 	    for(i=0 ; i<6 ; i++)
-	      pt->bdryinfo[i] = 0;  			
+	      pt->bdryinfo[i] = 0;
 		}
   }
-  if (ned) M_free(hed.item); 
+  if (ned) M_free(hed.item);
 
-  /*read corners*/ 
+  /*read corners*/
   if (ncor) {
     rewind(inm);
-    fseek(inm,posncor,SEEK_SET); 
+    fseek(inm,posncor,SEEK_SET);
     mesh->ncor = ncor;
-    for (k=1; k<=ncor; k++) { 
+    for (k=1; k<=ncor; k++) {
       if (!bin)
-        fscanf(inm,"%d",&ref);
+        if(fscanf(inm,"%d",&ref)!=1){ fclose(inm); return(0); }
       else {
-        fread(&ref,sw,1,inm); 
-        if(iswp) ref=MMG_swapbin(ref);    
-      }  
+        if(fread(&ref,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) ref=MMG_swapbin(ref);
+      }
       ppt = &mesh->point[ref];
-      ppt->geom = M_CORNER; 
-    } 
+      ppt->geom = M_CORNER;
+    }
   }
-   
-	
-  if ( abs(mesh->info.option)==10 ) { 
+
+
+  if ( abs(mesh->info.option)==10 ) {
     if(bin) {
       printf("NOT SUPPORTED\n");
       exit(0);
-    } 
-	  if ( !MMG_zaldy4(&hed2,3*npris+6*nhex) ) {
-      if ( mesh->info.ddebug )  fprintf(stdout,"  ## MEMORY ALLOCATION PROBLEM : PRISM IGNORED\n"); 
-			npris = 0;
-			nhex  = 0;
-    }   
+    }
+    if ( !MMG_zaldy4(&hed2,3*npris+6*nhex) ) {
+      if ( mesh->info.ddebug )  fprintf(stdout,"  ## MEMORY ALLOCATION PROBLEM : PRISM IGNORED\n");
+      npris = 0;
+      nhex  = 0;
+    }
 
     /*read hexa and transform to tetra*/
     rewind(inm);
     fseek(inm,posnhex,SEEK_SET);
-    for (k=1; k<=nhex; k++) {   
-      fscanf(inm,"%d %d %d %d %d %d %d %d %d",&p0,&p1,&p2,&p3,&p4,&p5,&p6,&p7,&ref); 
-      //fscanf(inm,"%d %d %d %d %d %d %d %d %d",&p0,&p4,&p2,&p1,&p3,&p5,&p6,&p7,&ref); 
-      //printf("hex %d : %d %d %d %d %d %d %d %d\n",k,p0,p1,p2,p3,p4,p5,p6,p7);   
-			MMG_cuthex(mesh,&hed2,netmp+(k-1)*6,p0,p1,p2,p3,p4,p5,p6,p7,ref);
-    }  
-    
+    for (k=1; k<=nhex; k++) {
+      if(fscanf(inm,"%d %d %d %d %d %d %d %d %d",&p0,&p1,&p2,&p3,&p4,&p5,&p6,&p7,&ref)!=9){ fclose(inm); return(0); }
+      //if(fscanf(inm,"%d %d %d %d %d %d %d %d %d",&p0,&p4,&p2,&p1,&p3,&p5,&p6,&p7,&ref)!=9){ fclose(inm); return(0); }
+      //printf("hex %d : %d %d %d %d %d %d %d %d\n",k,p0,p1,p2,p3,p4,p5,p6,p7);
+      MMG_cuthex(mesh,&hed2,netmp+(k-1)*6,p0,p1,p2,p3,p4,p5,p6,p7,ref);
+    }
+
     /*read prism and transform to tetra
-		---> compatibility pbs ==> hash edge and switch case*/  
+		---> compatibility pbs ==> hash edge and switch case*/
     rewind(inm);
-    fseek(inm,posnpris,SEEK_SET); 
-		nimp = 0; 
-		ne = netmp+6*nhex;
+    fseek(inm,posnpris,SEEK_SET);
+    nimp = 0;
+    ne = netmp+6*nhex;
     for (k=1; k<=npris; k++) {
-      fscanf(inm,"%d %d %d %d %d %d %d",&p0,&p1,&p2,&p3,&p4,&p5,&ref); 
-			if(!MMG_cutprism(mesh,&hed2,ne,p0,p1,p2,p3,p4,p5,ref))
-			{
-				if(mesh->info.imprim < 0 ) fprintf(stdout,"DECOMPOSITION PRISM INVALID \n\n"); 
-				mesh->ne += 5;
-				ne += 8;
-				nimp++; 
-				continue;
-			}
-			ne += 3;
+      if(fscanf(inm,"%d %d %d %d %d %d %d",&p0,&p1,&p2,&p3,&p4,&p5,&ref)!=7){ fclose(inm); return(0); }
+      if(!MMG_cutprism(mesh,&hed2,ne,p0,p1,p2,p3,p4,p5,ref)) {
+        if(mesh->info.imprim < 0 ) fprintf(stdout,"DECOMPOSITION PRISM INVALID \n\n");
+        mesh->ne += 5;
+        ne += 8;
+        nimp++;
+        continue;
+      }
+      ne += 3;
     }
-		if(abs(mesh->info.imprim) > 3 )fprintf(stdout,"     %d INVALID DECOMPOSITION\n\n",nimp);
+    if(abs(mesh->info.imprim) > 3 )fprintf(stdout,"     %d INVALID DECOMPOSITION\n\n",nimp);
   }
-  
+
   if ( abs(mesh->info.imprim) > 3 ) {
     fprintf(stdout,"     NUMBER OF GIVEN VERTICES   %8d\n",mesh->npfixe);
     if ( mesh->ntfixe )
@@ -563,18 +562,18 @@ int MMG_loadMesh(pMesh mesh,char *filename) {
 
 
 /* load solution (metric) */
-int MMG_loadSol(pSol sol,char *filename,int npmax) { 
-  FILE       *inm;   
+int MMG_loadSol(pSol sol,char *filename,int npmax) {
+  FILE       *inm;
   float       fsol;
-  double      tmp;       
+  double      tmp;
   int         binch,bdim,iswp;
   int         k,i,isol,type,bin,dim,btyp,bpos;
   long        posnp;
   char        *ptr,data[128],chaine[128];
-  
-  posnp = 0; 
+
+  posnp = 0;
   bin   = 0;
-  iswp  = 0; 
+  iswp  = 0;
 
   strcpy(data,filename);
   ptr = strstr(data,".mesh");
@@ -593,45 +592,45 @@ int MMG_loadSol(pSol sol,char *filename,int npmax) {
   }
   fprintf(stdout,"  %%%% %s OPENED\n",data);
 
-   
-  if(!bin) {   
+
+  if(!bin) {
     strcpy(chaine,"DDD");
-    while(fscanf(inm,"%s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) { 
+    while(fscanf(inm,"%s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
       if(!strncmp(chaine,"Dimension",strlen("Dimension"))) {
-          fscanf(inm,"%d",&dim);
+          if(fscanf(inm,"%d",&dim)!=1){ fclose(inm); return(0); }
           if(dim!=3) {
-            fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim); 
+            fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim);
             return(1);
           }
           continue;
       } else if(!strncmp(chaine,"SolAtVertices",strlen("SolAtVertices"))) {
-        fscanf(inm,"%d",&sol->np); 
-        fscanf(inm,"%d",&type); 
+        if(fscanf(inm,"%d",&sol->np)!=1){ fclose(inm); return(0); }
+        if(fscanf(inm,"%d",&type)!=1){ fclose(inm); return(0); }
         if(type!=1) {
           fprintf(stdout,"SEVERAL SOLUTION => IGNORED : %d\n",type);
           return(1);
         }
-        fscanf(inm,"%d",&btyp);
+        if(fscanf(inm,"%d",&btyp)!=1){ fclose(inm); return(0); }
         posnp = ftell(inm);
         break;
-      } 
-    }            
-  } else {     
-    fread(&binch,sw,1,inm); 
-    iswp=0;   
-    if(binch==16777216) iswp=1;    
+      }
+    }
+  } else {
+    if(fread(&binch,sw,1,inm)!=1){ fclose(inm); return(0); }
+    iswp=0;
+    if(binch==16777216) iswp=1;
     else if(binch!=1) {
       fprintf(stdout,"BAD FILE ENCODING\n");
-    } 
-    fread(&sol->ver,sw,1,inm); 
-    if(iswp) sol->ver = MMG_swapbin(sol->ver); 
+    }
+    if(fread(&sol->ver,sw,1,inm)!=1){ fclose(inm); return(0); }
+    if(iswp) sol->ver = MMG_swapbin(sol->ver);
     while(fread(&binch,sw,1,inm)!=EOF && binch!=54 ) {
-      if(iswp) binch=MMG_swapbin(binch);      
-      if(binch==54) break;  
+      if(iswp) binch=MMG_swapbin(binch);
+      if(binch==54) break;
       if(binch==3) {  //Dimension
-        fread(&bdim,sw,1,inm);  //NulPos=>20
+        if(fread(&bdim,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos=>20
         if(iswp) bdim=MMG_swapbin(bdim);
-        fread(&bdim,sw,1,inm);
+        if(fread(&bdim,sw,1,inm)!=1){ fclose(inm); return(0); }
         if(iswp) bdim=MMG_swapbin(bdim);
         if(bdim!=3) {
           fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim);
@@ -640,29 +639,29 @@ int MMG_loadSol(pSol sol,char *filename,int npmax) {
         }
         continue;
       } else if(binch==62) {  //SolAtVertices
-        fread(&binch,sw,1,inm); //NulPos
+        if(fread(&binch,sw,1,inm)!=1){ fclose(inm); return(0); } //NulPos
         if(iswp) binch=MMG_swapbin(binch);
-        fread(&sol->np,sw,1,inm); 
+        if(fread(&sol->np,sw,1,inm)!=1){ fclose(inm); return(0); }
         if(iswp) sol->np=MMG_swapbin(sol->np);
-        fread(&binch,sw,1,inm); //nb sol
+        if(fread(&binch,sw,1,inm)!=1){ fclose(inm); return(0); } //nb sol
         if(iswp) binch=MMG_swapbin(binch);
         if(binch!=1) {
           fprintf(stdout,"SEVERAL SOLUTION => IGNORED : %d\n",type);
           return(1);
         }
-        fread(&btyp,sw,1,inm); //typsol
+        if(fread(&btyp,sw,1,inm)!=1){ fclose(inm); return(0); } //typsol
         if(iswp) btyp=MMG_swapbin(btyp);
         posnp = ftell(inm);
         break;
       } else {
-        fread(&bpos,sw,1,inm); //Pos 
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //Pos
         if(iswp) bpos=MMG_swapbin(bpos);
         rewind(inm);
-        fseek(inm,bpos,SEEK_SET);        
-      } 
-    }            
-    
-  }       
+        fseek(inm,bpos,SEEK_SET);
+      }
+    }
+
+  }
   if ( !sol->np ) {
     fprintf(stdout,"  ** MISSING DATA\n");
     return(1);
@@ -673,7 +672,7 @@ int MMG_loadSol(pSol sol,char *filename,int npmax) {
     sol->np = 0;
     return(1);
   }
-  
+
   sol->offset = (btyp==1) ? 1 : 6;
 
   if ( abs(MMG_imprim) > 5 )
@@ -690,31 +689,31 @@ int MMG_loadSol(pSol sol,char *filename,int npmax) {
   /* read mesh solutions */
   sol->npfixe = sol->np;
   rewind(inm);
-  fseek(inm,posnp,SEEK_SET);  
+  fseek(inm,posnp,SEEK_SET);
   for (k=1; k<=sol->np; k++) {
     isol = (k-1) * sol->offset + 1;
-    if (sol->ver == 1) { 
+    if (sol->ver == 1) {
       for (i=0; i<sol->offset; i++) {
         if(!bin){
-          fscanf(inm,"%f",&fsol);    
+          if(fscanf(inm,"%f",&fsol)!=1){ fclose(inm); return(0); }
           sol->met[isol + i] = (double) fsol;
         } else {
-          fread(&fsol,sw,1,inm);             
+          if(fread(&fsol,sw,1,inm)!=1){ fclose(inm); return(0); }
           if(iswp) fsol=MMG_swapf(fsol);
           sol->met[isol + i] = (double) fsol;
         }
-      } 
+      }
     } else {
       for (i=0; i<sol->offset; i++) {
         if(!bin){
-          fscanf(inm,"%lf",&sol->met[isol + i]); 
+          if(fscanf(inm,"%lf",&sol->met[isol + i])!=1){ fclose(inm); return(0); }
 
         } else {
-          fread(&sol->met[isol + i],sd,1,inm);       
+          if(fread(&sol->met[isol + i],sd,1,inm)!=1){ fclose(inm); return(0); }
           if(iswp) sol->met[isol + i]=MMG_swapd(sol->met[isol + i]);
-        } 
-      } 
-    }             
+        }
+      }
+    }
     /* MMG_swap data */
     if ( sol->offset == 6 ) {
       tmp                = sol->met[isol + 2];
@@ -726,24 +725,24 @@ int MMG_loadSol(pSol sol,char *filename,int npmax) {
   if ( abs(MMG_imprim) > 3 )
     fprintf(stdout,"     NUMBER OF GIVEN DATA       %8d\n",sol->npfixe);
 
-  fclose(inm);   
-  return(1);  
+  fclose(inm);
+  return(1);
 }
 
 
 int MMG_loadVect(pMesh mesh,char *filename,int npmax) {
-  FILE       *inm;   
+  FILE       *inm;
   pDispl       pd;
   float       fsol;
   int         binch,bdim,iswp;
   int         k,i,type,bin,dim,btyp,bpos,iadr;
   long        posnp;
   char        *ptr,data[128],chaine[128];
-  
+
   pd = mesh->disp;
-  
-  posnp = 0; 
-  bin   = 0; 
+
+  posnp = 0;
+  bin   = 0;
   iswp  = 0;
 
   strcpy(data,filename);
@@ -763,46 +762,46 @@ int MMG_loadVect(pMesh mesh,char *filename,int npmax) {
   }
   fprintf(stdout,"  %%%% %s OPENED\n",data);
 
-   
-  if(!bin) {   
+
+  if(!bin) {
     strcpy(chaine,"DDD");
-    while(fscanf(inm,"%s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) { 
+    while(fscanf(inm,"%s",&chaine[0])!=EOF && strncmp(chaine,"End",strlen("End")) ) {
       if(!strncmp(chaine,"Dimension",strlen("Dimension"))) {
-          fscanf(inm,"%d",&dim);
+          if(fscanf(inm,"%d",&dim)!=1){ fclose(inm); return(0); }
           if(dim!=3) {
-            fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim); 
+            fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim);
             return(1);
           }
           continue;
       } else if(!strncmp(chaine,"SolAtVertices",strlen("SolAtVertices"))) {
-        fscanf(inm,"%d",&pd->np); 
-        fscanf(inm,"%d",&type); 
+        if(fscanf(inm,"%d",&pd->np)!=1){ fclose(inm); return(0); }
+        if(fscanf(inm,"%d",&type)!=1){ fclose(inm); return(0); }
         if(type!=1) {
           fprintf(stdout,"SEVERAL SOLUTION => IGNORED : %d\n",type);
           return(1);
         }
-        fscanf(inm,"%d",&btyp);
+        if(fscanf(inm,"%d",&btyp)!=1){ fclose(inm); return(0); }
         posnp = ftell(inm);
         break;
-      } 
-    }            
-  } else {     
-    fread(&pd->ver,sw,1,inm); 
-    iswp=0;   
-    if(pd->ver==16777216) iswp=1;    
+      }
+    }
+  } else {
+    if(fread(&pd->ver,sw,1,inm)!=1){ fclose(inm); return(0); }
+    iswp=0;
+    if(pd->ver==16777216) iswp=1;
     else if(pd->ver!=1) {
       fprintf(stdout,"BAD FILE ENCODING\n");
-    } 
-    fread(&pd->ver,sw,1,inm); 
-    if(iswp) pd->ver = MMG_swapbin(pd->ver); 
+    }
+    if(fread(&pd->ver,sw,1,inm)!=1){ fclose(inm); return(0); }
+    if(iswp) pd->ver = MMG_swapbin(pd->ver);
     while(fread(&binch,sw,1,inm)!=EOF && binch!=54 ) {
-      if(iswp) binch=MMG_swapbin(binch);      
-      if(binch==54) break;  
+      if(iswp) binch=MMG_swapbin(binch);
+      if(binch==54) break;
       if(binch==3) {  //Dimension
-        fread(&bdim,sw,1,inm);  //Pos=>20
-        if(iswp) bdim=MMG_swapbin(bdim);      
-        fread(&bdim,sw,1,inm);
-        if(iswp) bdim=MMG_swapbin(bdim);      
+        if(fread(&bdim,sw,1,inm)!=1){ fclose(inm); return(0); } //Pos=>20
+        if(iswp) bdim=MMG_swapbin(bdim);
+        if(fread(&bdim,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) bdim=MMG_swapbin(bdim);
         if(bdim!=3) {
           fprintf(stdout,"BAD SOL DIMENSION : %d\n",dim);
           exit(0);
@@ -810,29 +809,29 @@ int MMG_loadVect(pMesh mesh,char *filename,int npmax) {
         }
         continue;
       } else if(binch==62) {  //SolAtVertices
-        fread(&binch,sw,1,inm); //Pos
-        if(iswp) binch=MMG_swapbin(binch);      
-        fread(&pd->np,sw,1,inm); 
-        if(iswp) pd->np=MMG_swapbin(pd->np);      
-        fread(&binch,sw,1,inm); //nb sol
-        if(iswp) binch=MMG_swapbin(binch);      
+        if(fread(&binch,sw,1,inm)!=1){ fclose(inm); return(0); } //Pos
+        if(iswp) binch=MMG_swapbin(binch);
+        if(fread(&pd->np,sw,1,inm)!=1){ fclose(inm); return(0); }
+        if(iswp) pd->np=MMG_swapbin(pd->np);
+        if(fread(&binch,sw,1,inm)!=1){ fclose(inm); return(0); } //nb sol
+        if(iswp) binch=MMG_swapbin(binch);
         if(binch!=1) {
           fprintf(stdout,"SEVERAL SOLUTION => IGNORED : %d\n",type);
           return(1);
         }
-        fread(&btyp,sw,1,inm); //typsol
-        if(iswp) btyp=MMG_swapbin(btyp);      
+        if(fread(&btyp,sw,1,inm)!=1){ fclose(inm); return(0); } //typsol
+        if(iswp) btyp=MMG_swapbin(btyp);
         posnp = ftell(inm);
         break;
       } else {
-        fread(&bpos,sw,1,inm); //Pos 
-        if(iswp) bpos=MMG_swapbin(bpos);      
+        if(fread(&bpos,sw,1,inm)!=1){ fclose(inm); return(0); } //Pos
+        if(iswp) bpos=MMG_swapbin(bpos);
         rewind(inm);
-        fseek(inm,bpos,SEEK_SET);        
-      } 
-    }            
-    
-  }       
+        fseek(inm,bpos,SEEK_SET);
+      }
+    }
+
+  }
   if ( !pd->np ) {
     fprintf(stdout,"  ** MISSING DATA\n");
     return(0);
@@ -855,46 +854,46 @@ int MMG_loadVect(pMesh mesh,char *filename,int npmax) {
   fseek(inm,posnp,SEEK_SET);
   for (k=1; k<=pd->np; k++) {
     iadr = (k - 1) * 3 + 1;
-    if (pd->ver < 2) { 
+    if (pd->ver < 2) {
       for (i=0; i<3; i++) {
         if(!bin){
-          fscanf(inm,"%f",&fsol); 
+          if(fscanf(inm,"%f",&fsol)!=1){ fclose(inm); return(0); }
           pd->mv[iadr + i] = (double) fsol;
         } else {
-          fread(&fsol,sw,1,inm);             
-          if(iswp) fsol=MMG_swapf(fsol);      
+          if(fread(&fsol,sw,1,inm)!=1){ fclose(inm); return(0); }
+          if(iswp) fsol=MMG_swapf(fsol);
           pd->mv[iadr + i] = (double) fsol;
         }
-      } 
+      }
     } else {
       for (i=0; i<3; i++) {
         if(!bin){
-          fscanf(inm,"%lf",&pd->mv[iadr + i]); 
+          if(fscanf(inm,"%lf",&pd->mv[iadr + i])!=1){ fclose(inm); return(0); }
         } else {
-          fread(&pd->mv[iadr + i],sd,1,inm);
-          if(iswp) pd->mv[iadr + i]=MMG_swapd(pd->mv[iadr + i]);      
-        } 
-      } 
-    }             
+          if(fread(&pd->mv[iadr + i],sd,1,inm)!=1){ fclose(inm); return(0); }
+          if(iswp) pd->mv[iadr + i]=MMG_swapd(pd->mv[iadr + i]);
+        }
+      }
+    }
   }
 
   if ( abs(mesh->info.imprim) > 3 )
     fprintf(stdout,"     NUMBER OF GIVEN DATA       %8d\n",pd->np);
 
-  fclose(inm); 
+  fclose(inm);
   return(1);
 }
 
 
 /* save mesh to disk */
-int MMG_saveMesh(pMesh mesh,char *filename) {  
-  FILE*        inm; 
+int MMG_saveMesh(pMesh mesh,char *filename) {
+  FILE*        inm;
 	Hedge				 hed;
   pPoint       ppt;
   pTria        pt1;
   pTetra       pt;
   int          i,k,np,ne,nc,ned,*cor,*ed,ref,bin,bpos;
-  char        *ptr,data[128],chaine[128]; 
+  char        *ptr,data[128],chaine[128];
   int          binch,nu1,nu2;
   mesh->ver = 2; //double precision
   bin = 0;
@@ -911,24 +910,24 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
         return(0);
       }
     } else {
-      bin = 1;   
+      bin = 1;
     }
   }
-  else { 
+  else {
     ptr = strstr(data,".meshb");
     if( ptr ) bin = 1;
     if( !(inm = fopen(data,"w")) ) {
       fprintf(stderr,"  ** UNABLE TO OPEN %s.\n",data);
       return(0);
-    } 
+    }
   }
   fprintf(stdout,"  %%%% %s OPENED\n",data);
 
   /*entete fichier*/
   if(!bin) {
-    strcpy(&chaine[0],"MeshVersionFormatted 2\n"); 
+    strcpy(&chaine[0],"MeshVersionFormatted 2\n");
     fprintf(inm,"%s",chaine);
-    strcpy(&chaine[0],"\n\nDimension 3\n"); 
+    strcpy(&chaine[0],"\n\nDimension 3\n");
     fprintf(inm,"%s ",chaine);
   } else {
     binch = 1; //MeshVersionFormatted
@@ -941,33 +940,33 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
     fwrite(&bpos,sw,1,inm);
     binch = 3;
     fwrite(&binch,sw,1,inm);
-    
+
   }
 
   /* compact vertices */
-  if(mesh->ncor) {   
+  if(mesh->ncor) {
     cor = (int*) M_calloc(mesh->ncor,sizeof(int),"MMG_savemesh");
-    assert(cor);   
+    assert(cor);
   }
-  if(mesh->ned) {   
+  if(mesh->ned) {
 	  if ( !MMG_zaldy4(&hed,mesh->ned) ) {
-      if ( mesh->info.ddebug )  fprintf(stdout,"  ## MEMORY ALLOCATION PROBLEM : EXPORT EDGES IGNORED\n"); 
+      if ( mesh->info.ddebug )  fprintf(stdout,"  ## MEMORY ALLOCATION PROBLEM : EXPORT EDGES IGNORED\n");
 			mesh->ned = 0;
-    }   
+    }
     ed = (int*)M_calloc(2*mesh->ned,sizeof(int),"MMG_savemesh");
-    assert(ed);   
+    assert(ed);
   }
-  np = 0; 
+  np = 0;
   nc = 0;
   for (k=1; k<=mesh->np; k++) {
     ppt = &mesh->point[k];
-    if ( ppt->tag & M_UNUSED )  continue;  
-		ppt->tmp = ++np;  
+    if ( ppt->tag & M_UNUSED )  continue;
+		ppt->tmp = ++np;
     if ( ppt->geom & M_CORNER )  cor[nc++] = ppt->tmp;
-  } 
+  }
 	assert(mesh->ncor==nc);
   if(!bin) {
-    strcpy(&chaine[0],"\n\nVertices\n"); 
+    strcpy(&chaine[0],"\n\nVertices\n");
     fprintf(inm,"%s",chaine);
     fprintf(inm,"%d\n",np);
   } else {
@@ -975,27 +974,27 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
     fwrite(&binch,sw,1,inm);
     bpos += 12+(1+3*mesh->ver)*4*np; //NullPos
     fwrite(&bpos,sw,1,inm);
-    fwrite(&np,sw,1,inm);    
+    fwrite(&np,sw,1,inm);
   }
   for(k=1; k<=mesh->np; k++) {
     ppt = &mesh->point[k];
-    if ( ppt->tag & M_UNUSED )  continue;  
+    if ( ppt->tag & M_UNUSED )  continue;
 		//if(ppt->tmp==52453) printf("point %d --> %d\n",ppt->tmp,k);
     if(!bin) {
       fprintf(inm,"%.15lg %.15lg %.15lg %d\n",ppt->c[0],ppt->c[1],ppt->c[2],ppt->ref);
     } else {
-      fwrite((unsigned char*)&ppt->c[0],sd,1,inm);    
-      fwrite((unsigned char*)&ppt->c[1],sd,1,inm);    
-      fwrite((unsigned char*)&ppt->c[2],sd,1,inm);    
-      fwrite((unsigned char*)&ppt->ref,sw,1,inm);    
+      fwrite((unsigned char*)&ppt->c[0],sd,1,inm);
+      fwrite((unsigned char*)&ppt->c[1],sd,1,inm);
+      fwrite((unsigned char*)&ppt->c[2],sd,1,inm);
+      fwrite((unsigned char*)&ppt->ref,sw,1,inm);
     }
   }
 
-  /* rebuild triangles tabular and write triangles */ 
+  /* rebuild triangles tabular and write triangles */
   mesh->nt = 0;
   if(MMG_markBdry(mesh)) {
     if(!bin) {
-      strcpy(&chaine[0],"\n\nTriangles\n"); 
+      strcpy(&chaine[0],"\n\nTriangles\n");
       fprintf(inm,"%s",chaine);
       fprintf(inm,"%d \n",mesh->nt);
     } else {
@@ -1003,26 +1002,26 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
       fwrite(&binch,sw,1,inm);
       bpos += 12+16*mesh->nt; //Pos
       fwrite(&bpos,sw,1,inm);
-      fwrite(&mesh->nt,sw,1,inm);    
+      fwrite(&mesh->nt,sw,1,inm);
     }
     for (k=1; k<=mesh->nt; k++) {
       pt1  = &mesh->tria[k];
-  	    ref  = pt1->ref;    
+  	    ref  = pt1->ref;
       if(!bin) {
         fprintf(inm,"%d %d %d %d\n",mesh->point[pt1->v[0]].tmp,mesh->point[pt1->v[1]].tmp
     							  ,mesh->point[pt1->v[2]].tmp,ref);
       } else {
-        fwrite(&mesh->point[pt1->v[0]].tmp,sw,1,inm);    
-        fwrite(&mesh->point[pt1->v[1]].tmp,sw,1,inm);    
-        fwrite(&mesh->point[pt1->v[2]].tmp,sw,1,inm);    
-        fwrite(&ref,sw,1,inm);    
+        fwrite(&mesh->point[pt1->v[0]].tmp,sw,1,inm);
+        fwrite(&mesh->point[pt1->v[1]].tmp,sw,1,inm);
+        fwrite(&mesh->point[pt1->v[2]].tmp,sw,1,inm);
+        fwrite(&ref,sw,1,inm);
       }
     }
-  }   
- 
+  }
+
   /* write tetrahedra */
-  ne = 0; 
-	ned = 0;  
+  ne = 0;
+	ned = 0;
 	//printf("avt %d\n",mesh->ned);
   for (k=1; k<=mesh->ne; k++) {
     pt = &mesh->tetra[k];
@@ -1036,15 +1035,15 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
 		  			ed[2*ned] = (mesh->point[nu1]).tmp;
 		  			ed[2*ned + 1] = (mesh->point[nu2]).tmp;
 		  			ned++;
-		  		} 
+		  		}
 		  	}
-		  } 
+		  }
 		}
-	  ne++;  
+	  ne++;
   }
 	//printf("ned %d\n",ned);
   if(!bin) {
-    strcpy(&chaine[0],"\n\nTetrahedra\n"); 
+    strcpy(&chaine[0],"\n\nTetrahedra\n");
     fprintf(inm,"%s",chaine);
     fprintf(inm,"%d\n",ne);
   } else {
@@ -1052,29 +1051,29 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
     fwrite(&binch,sw,1,inm);
     bpos += 12 + 20*ne;//Pos
     fwrite(&bpos,sw,1,inm);
-    fwrite((unsigned char*)&ne,sw,1,inm);    
-  } 
+    fwrite((unsigned char*)&ne,sw,1,inm);
+  }
 	ne=0;
   for (k=1; k<=mesh->ne; k++) {
     pt = &mesh->tetra[k];
-    if ( !pt->v[0] )  continue;  
-		ne++; 
-    ref = pt->ref;    
+    if ( !pt->v[0] )  continue;
+		ne++;
+    ref = pt->ref;
     if(!bin) {
       fprintf(inm,"%d %d %d %d %d\n",mesh->point[pt->v[0]].tmp,mesh->point[pt->v[1]].tmp
   							   ,mesh->point[pt->v[2]].tmp,mesh->point[pt->v[3]].tmp,ref);
     } else {
-      fwrite(&mesh->point[pt->v[0]].tmp,sw,1,inm);    
-      fwrite(&mesh->point[pt->v[1]].tmp,sw,1,inm);    
-      fwrite(&mesh->point[pt->v[2]].tmp,sw,1,inm);    
-      fwrite(&mesh->point[pt->v[3]].tmp,sw,1,inm);    
-      fwrite(&ref,sw,1,inm);    
+      fwrite(&mesh->point[pt->v[0]].tmp,sw,1,inm);
+      fwrite(&mesh->point[pt->v[1]].tmp,sw,1,inm);
+      fwrite(&mesh->point[pt->v[2]].tmp,sw,1,inm);
+      fwrite(&mesh->point[pt->v[3]].tmp,sw,1,inm);
+      fwrite(&ref,sw,1,inm);
     }
-  }  
-     
+  }
+
   if(mesh->ned) {
     if(!bin) {
-      strcpy(&chaine[0],"\n\nEdges\n"); 
+      strcpy(&chaine[0],"\n\nEdges\n");
       fprintf(inm,"%s",chaine);
       fprintf(inm,"%d\n",ned);
     } else {
@@ -1082,49 +1081,49 @@ int MMG_saveMesh(pMesh mesh,char *filename) {
       fwrite(&binch,sw,1,inm);
       bpos += 12 + 3*4*ned;//Pos
       fwrite(&bpos,sw,1,inm);
-      fwrite((unsigned char*)&ned,sw,1,inm);    
-    } 
+      fwrite((unsigned char*)&ned,sw,1,inm);
+    }
   	  for (k=0; k<ned; k++) {
-   	    ref = 0;    
+   	    ref = 0;
   	    if(!bin) {
   	      fprintf(inm,"%d %d %d \n",ed[2*k],ed[2*k+1],ref);
   	    } else {
-  	      fwrite(&ed[2*k],sw,1,inm);    
-  	      fwrite(&ed[2*k+1],sw,1,inm);    
-  	      fwrite(&ref,sw,1,inm);    
+  	      fwrite(&ed[2*k],sw,1,inm);
+  	      fwrite(&ed[2*k+1],sw,1,inm);
+  	      fwrite(&ref,sw,1,inm);
   	    }
   	  }
   	  M_free(hed.item);
   }
-  
+
   /* write corners */
   if(!bin) {
-    strcpy(&chaine[0],"\n\nCorners\n"); 
+    strcpy(&chaine[0],"\n\nCorners\n");
     fprintf(inm,"%s",chaine);
     fprintf(inm,"%d\n",mesh->ncor);
   } else {
     binch = 13; //Corners
     fwrite(&binch,sw,1,inm);
-    bpos += 12 + 4*mesh->ncor;//Pos  
+    bpos += 12 + 4*mesh->ncor;//Pos
     fwrite(&bpos,sw,1,inm);
-    fwrite((unsigned char*)&mesh->ncor,sw,1,inm);    
+    fwrite((unsigned char*)&mesh->ncor,sw,1,inm);
   }
   for (k=0; k<mesh->ncor; k++) {
     if(!bin) {
       fprintf(inm,"%d \n",cor[k]);
     } else {
-      fwrite(&cor[k],sw,1,inm);    
+      fwrite(&cor[k],sw,1,inm);
     }
-  }  
+  }
   /*fin fichier*/
   if(!bin) {
-    strcpy(&chaine[0],"\n\nEnd\n"); 
+    strcpy(&chaine[0],"\n\nEnd\n");
     fprintf(inm,"%s",chaine);
   } else {
     binch = 54; //End
     fwrite(&binch,sw,1,inm);
   }
-  fclose(inm); 
+  fclose(inm);
   if(mesh->ncor) M_free(cor);
   if ( mesh->info.imprim ) {
     fprintf(stdout,"     NUMBER OF GIVEN VERTICES   %8d\n",mesh->npfixe);
@@ -1151,7 +1150,7 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
   float        fsol;
   double       tmp;
   int          i,k,nbl,isol,bin,bpos,typ;
-  char        *ptr,data[128],chaine[128]; 
+  char        *ptr,data[128],chaine[128];
   int          binch;
 
   if ( !sol->np )  return(1);
@@ -1168,21 +1167,21 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
 	    ptr = strstr(data,".solb");
 	    if ( ptr ) {
 	      *ptr = '\0';
-	      bin  = 1;	
+	      bin  = 1;
       } else {
 			  ptr = strstr(data,".sol");
 			  if ( ptr ) {
 			    *ptr = '\0';
-			    bin  = 0;	
+			    bin  = 0;
 			  }
 			}
-    } 
+    }
   }
-  if ( bin ) 
+  if ( bin )
     strcat(data,".solb");
   else
     strcat(data,".sol");
-  
+
   sol->ver = 2;
   if( bin && !(inm = fopen(data,"wb")) ) {
     fprintf(stderr,"  ** UNABLE TO OPEN %s.\n",data);
@@ -1197,9 +1196,9 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
 
   /*entete fichier*/
   if(!bin) {
-    strcpy(&chaine[0],"MeshVersionFormatted 2\n"); 
+    strcpy(&chaine[0],"MeshVersionFormatted 2\n");
     fprintf(inm,"%s",chaine);
-    strcpy(&chaine[0],"\n\nDimension 3\n"); 
+    strcpy(&chaine[0],"\n\nDimension 3\n");
     fprintf(inm,"%s ",chaine);
   } else {
     binch = 1; //MeshVersionFormatted
@@ -1212,7 +1211,7 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
     fwrite(&bpos,sw,1,inm);
     binch = 3;
     fwrite(&binch,sw,1,inm);
-    
+
   }
 
 
@@ -1235,9 +1234,9 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
     if ( ppt->tag & M_UNUSED )  continue;
 	nbl++;
   }
-  
+
   if(!bin) {
-    strcpy(&chaine[0],"\n\nSolAtVertices\n"); 
+    strcpy(&chaine[0],"\n\nSolAtVertices\n");
     fprintf(inm,"%s",chaine);
     fprintf(inm,"%d\n",nbl);
     fprintf(inm,"%d %d\n",1,typ);
@@ -1246,7 +1245,7 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
     fwrite(&binch,sw,1,inm);
     bpos += 20+(sol->offset*sol->ver)*4*nbl; //Pos
     fwrite(&bpos,sw,1,inm);
-    fwrite(&nbl,sw,1,inm);    
+    fwrite(&nbl,sw,1,inm);
     binch = 1; //nb sol
     fwrite(&binch,sw,1,inm);
     binch = typ; //typ sol
@@ -1263,34 +1262,34 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
       sol->met[isol + 3] = tmp;
     }
     if (sol->ver < 2) {
-      if(!bin) { 
+      if(!bin) {
         for (i=0; i<sol->offset; i++) {
           fsol = (float) sol->met[isol + i];
           fprintf(inm,"%f ",fsol);
-        } 
-        fprintf(inm,"\n");  
+        }
+        fprintf(inm,"\n");
       } else {
-        for (i=0; i<sol->offset; i++) { 
+        for (i=0; i<sol->offset; i++) {
           fsol = (float) sol->met[isol + i];
           fwrite(&fsol,sw,1,inm);
-        }    
+        }
       }
     } else {
-      if(!bin) { 
+      if(!bin) {
         for (i=0; i<sol->offset; i++)
-          fprintf(inm,"%.15lg ",sol->met[isol + i]); 
-        fprintf(inm,"\n");  
+          fprintf(inm,"%.15lg ",sol->met[isol + i]);
+        fprintf(inm,"\n");
       } else {
         for (i=0; i<sol->offset; i++)
-          fwrite(&sol->met[isol + i],sd,1,inm);    
+          fwrite(&sol->met[isol + i],sd,1,inm);
       }
-      
+
     }
   }
-  
+
   /*fin fichier*/
   if(!bin) {
-    strcpy(&chaine[0],"\n\nEnd\n"); 
+    strcpy(&chaine[0],"\n\nEnd\n");
     fprintf(inm,"%s",chaine);
   } else {
     binch = 54; //End
@@ -1302,12 +1301,12 @@ int MMG_saveSol(pMesh mesh,pSol sol,char *filename) {
 
 /*save the node speed : coornew-coorold/dt*/
 int MMG_saveVect(pMesh mesh,char *filename) {
-  FILE*        inm;  
+  FILE*        inm;
   pDispl        pd;
   pPoint       ppt;
   double       dsol,dd;
   int          i,k,nbl,bin,bpos,typ;
-  char        *ptr,data[128],chaine[128]; 
+  char        *ptr,data[128],chaine[128];
   unsigned char binch;
 
   pd      = mesh->disp;
@@ -1324,7 +1323,7 @@ int MMG_saveVect(pMesh mesh,char *filename) {
       bin  = 0;
     }
   }
-  if ( bin ) 
+  if ( bin )
     strcat(data,".o.solb");
   else
     strcat(data,".o.sol");
@@ -1341,9 +1340,9 @@ int MMG_saveVect(pMesh mesh,char *filename) {
 
   /*entete fichier*/
   if(!bin) {
-    strcpy(&chaine[0],"MeshVersionFormatted 2\n"); 
+    strcpy(&chaine[0],"MeshVersionFormatted 2\n");
     fprintf(inm,"%s",chaine);
-    strcpy(&chaine[0],"\n\nDimension 3\n"); 
+    strcpy(&chaine[0],"\n\nDimension 3\n");
     fprintf(inm,"%s ",chaine);
   } else {
     binch = 1; //MeshVersionFormatted
@@ -1356,7 +1355,7 @@ int MMG_saveVect(pMesh mesh,char *filename) {
     fwrite(&bpos,sw,1,inm);
     binch = 3;
     fwrite(&binch,sw,1,inm);
-    
+
   }
 	typ = 2;
 
@@ -1367,9 +1366,9 @@ int MMG_saveVect(pMesh mesh,char *filename) {
     if ( ppt->tag & M_UNUSED )  continue;
 	nbl++;
   }
-  
+
   if(!bin) {
-    strcpy(&chaine[0],"\n\nSolAtVertices\n"); 
+    strcpy(&chaine[0],"\n\nSolAtVertices\n");
     fprintf(inm,"%s",chaine);
     fprintf(inm,"%d\n",nbl);
     fprintf(inm,"%d %d\n",1,typ);
@@ -1378,34 +1377,34 @@ int MMG_saveVect(pMesh mesh,char *filename) {
     fwrite(&binch,sw,1,inm);
     bpos += 20+(3*pd->ver)*4*nbl; //Pos
     fwrite(&bpos,sw,1,inm);
-    fwrite(&nbl,sw,1,inm);    
+    fwrite(&nbl,sw,1,inm);
     binch = 1; //nb sol
     fwrite(&binch,sw,1,inm);
     binch = typ; //typ sol
     fwrite(&binch,sw,1,inm);
-  } 
-  
-  
-  dd = mesh->info.delta / (double)PRECI;  
+  }
+
+
+  dd = mesh->info.delta / (double)PRECI;
   fprintf(stdout,"        DT %e\n",mesh->info.dt);
   for (k=1; k<=mesh->np; k++) {
     ppt = &mesh->point[k];
-    if ( ppt->tag & M_UNUSED )  continue; 
+    if ( ppt->tag & M_UNUSED )  continue;
     for (i=0 ; i<3 ; i++) {
-      dsol = (ppt->c[i] - mesh->disp->cold[3*(k-1) + 1 + i]* dd - mesh->info.min[i])/mesh->info.dt; 
-      if(!bin) { 
-        fprintf(inm,"%.15lg ",dsol); 
+      dsol = (ppt->c[i] - mesh->disp->cold[3*(k-1) + 1 + i]* dd - mesh->info.min[i])/mesh->info.dt;
+      if(!bin) {
+        fprintf(inm,"%.15lg ",dsol);
       } else {
-        fwrite(&dsol,sd,1,inm);    
+        fwrite(&dsol,sd,1,inm);
       }
     }
-    if (!bin) fprintf(inm,"\n");  
+    if (!bin) fprintf(inm,"\n");
   }
-  
-  
+
+
   /*fin fichier*/
   if(!bin) {
-    strcpy(&chaine[0],"\n\nEnd\n"); 
+    strcpy(&chaine[0],"\n\nEnd\n");
     fprintf(inm,"%s",chaine);
   } else {
     binch = 54; //End
diff --git a/contrib/mpeg_encode/specifics.cpp b/contrib/mpeg_encode/specifics.cpp
index 67dddc0fa03f5f3c5723bb676d764b26cdc8b9e6..48a4e46402574b29c1321a3c75d5ca5b5c66c269 100644
--- a/contrib/mpeg_encode/specifics.cpp
+++ b/contrib/mpeg_encode/specifics.cpp
@@ -159,10 +159,10 @@ void Specifics_Init()
   FILE *specificsFP;
   
   sprintf(command, "/bin/rm -f %s.cpp", specificsFile);
-  system(command);
+  if(system(command));
   sprintf(command, "%s -P %s %s %s.cpp",
 	  CPP_LOC, specificsDefines, specificsFile, specificsFile);
-  system(command);
+  if(system(command));
   strcat(specificsFile, ".cpp");
   if ((specificsFP = fopen(specificsFile, "r")) == NULL) {
     throw "Cannot open specifics file";
@@ -170,7 +170,7 @@ void Specifics_Init()
   printf("Specifics file: %s\n", specificsFile);
   Parse_Specifics_File(specificsFP);
   sprintf(command, "/bin/rm -f %s.cpp", specificsFile);
-  system(command);
+  if(system(command));
 
 }
 
diff --git a/contrib/onelab/OnelabClients.cpp b/contrib/onelab/OnelabClients.cpp
index e79c8c62238dc2158ecf5811978dbf0e7695bff0..ad0e7bd41c9a8d973872c3a0059d991ae9a6d415 100644
--- a/contrib/onelab/OnelabClients.cpp
+++ b/contrib/onelab/OnelabClients.cpp
@@ -10,6 +10,18 @@
 #define PCLOSE pclose
 #endif
 
+#if defined(WIN32)
+static char dirSep='\\';
+static std::string cmdSep(" & ");
+static std::string removeCmd("del ");
+static std::string lsCmd("dir ");
+#else
+static char dirSep='/';
+static std::string cmdSep(" ; ");
+static std::string removeCmd("rm -rf ");
+static std::string lsCmd("ls ");
+#endif
+
 class onelabMetaModelServer : public GmshServer{
  private:
   localNetworkSolverClient *_client;
diff --git a/contrib/onelab/OnelabClients.h b/contrib/onelab/OnelabClients.h
index 4e4a20986c4a76f051783fbef3cda1d726bb5d01..0ef3c769e4dd58fb1016e460fe15c45ebe47692f 100644
--- a/contrib/onelab/OnelabClients.h
+++ b/contrib/onelab/OnelabClients.h
@@ -22,18 +22,6 @@ static std::string localFileTag("_");
 // Possible actions for clients
 enum parseMode {REGISTER, ANALYZE, COMPUTE, EXIT};
 
-#if defined(WIN32)
-static char dirSep='\\';
-static std::string cmdSep(" & ");
-static std::string removeCmd("del ");
-static std::string lsCmd("dir ");
-#else
-static char dirSep='/';
-static std::string cmdSep(" ; ");
-static std::string removeCmd("rm -rf ");
-static std::string lsCmd("ls ");
-#endif
-
 // TOOLS 
 
 std::string itoa(const int i);
diff --git a/contrib/onelab/metamodel.cpp b/contrib/onelab/metamodel.cpp
index 20535cff22aa7a2157f58640d792b93eb440ec44..d0de66ea9c78b17535689a1c568f9104180146d6 100644
--- a/contrib/onelab/metamodel.cpp
+++ b/contrib/onelab/metamodel.cpp
@@ -68,8 +68,8 @@ int metamodel(const std::string &action){
     std::string mystderr = FixWindowsQuotes(workingDir + "stderr.txt");
     OLMsg::Info("Redirecting stdout into <%s>",mystdout.c_str());
     OLMsg::Info("Redirecting stderr into <%s>",mystderr.c_str());
-    freopen(mystdout.c_str(),"w",stdout);
-    freopen(mystderr.c_str(),"w",stderr);
+    if(!freopen(mystdout.c_str(),"w",stdout)) return 0;
+    if(!freopen(mystderr.c_str(),"w",stderr)) return 0;
   }
 
   if( myModel->isTodo(ANALYZE)){