diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 1cb76601e4bb296cc083ae4200f291ea8e9521ef..67f79ae38c945af0af3a403e27c5e3a82f16b1b4 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -45,10 +45,6 @@
 #include "PView.h"
 #endif
 
-#if defined(HAVE_MESH)
-#include "periodical.h"
-#endif
-
 #if defined(HAVE_PARSER)
 #include "Parser.h"
 #endif
@@ -541,16 +537,6 @@ void GetOptions(int argc, char *argv[])
         else
           Msg::Fatal("Missing number of lloyd iterations");
       }
-#if defined(HAVE_MESH)
-      else if(!strcmp(argv[i] + 1, "microstructure")) {
-        i++;
-        if(argv[i]) microstructure(argv[i++]);
-      }
-      else if(!strcmp(argv[i] + 1, "computeBestSeeds")) {
-        i++;
-        if(argv[i]) computeBestSeeds(argv[i++]);
-      }
-#endif
       else if(!strcmp(argv[i] + 1, "nopopup")) {
         CTX::instance()->noPopup = 1;
         i++;
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 62a90298b27a74fe0a51bf1f5aac20a378fbe98b..ffbe77be989cf9b362d399bd02ad5694f50507ca 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -43,7 +43,6 @@ set(SRC
   DivideAndConquer.cpp
     Voronoi3D.cpp
     Levy3D.cpp
-    periodical.cpp
     directions3D.cpp
     filterElements.cpp
     yamakawa.cpp
diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index 83343629bac83735ce48bb3965f4eedcdf267959..5f37603afc7d6190cd15c071f95f51cc69265ba3 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -38,6 +38,7 @@ set(SRC
   CVTRemesh.cpp
   ShowNeighborElements.cpp
   GaussPoints.cpp
+  VoroMetal.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index 25599d209443687f15210a43603932b956bbc866..bfdd58b0def96e5d3853114bcad62f28765a8053 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -67,6 +67,7 @@
 #include "CVTRemesh.h"
 #include "ShowNeighborElements.h"
 #include "GaussPoints.h"
+#include "VoroMetal.h"
 
 // for testing purposes only :-)
 #undef HAVE_DLOPEN
@@ -272,6 +273,8 @@ void PluginManager::registerDefaultPlugins()
 #if defined(HAVE_MESH)
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
                       ("AnalyseCurvedMesh", GMSH_RegisterAnalyseCurvedMeshPlugin()));
+    allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+                      ("VoroMetal", GMSH_RegisterVoroMetalPlugin()));
 #endif
 #if defined(HAVE_REVOROPT)
     allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
diff --git a/Plugin/Tetrahedralize.cpp b/Plugin/Tetrahedralize.cpp
index 4145df62973b345a002117337bb6a29892073fa9..7bc315cff0a7eecf9d3f470b89fe6592aad246e8 100644
--- a/Plugin/Tetrahedralize.cpp
+++ b/Plugin/Tetrahedralize.cpp
@@ -152,7 +152,7 @@ PView *GMSH_TetrahedralizePlugin::execute(PView *v)
 
 PView *GMSH_TetrahedralizePlugin::execute(PView *v)
 {
-  Msg::Error("Plugin(Tetrahedralize requires mesh module");
+  Msg::Error("Plugin(Tetrahedralize) requires mesh module");
   return v;
 }
 
diff --git a/Mesh/periodical.cpp b/Plugin/VoroMetal.cpp
similarity index 76%
rename from Mesh/periodical.cpp
rename to Plugin/VoroMetal.cpp
index 21833c682cbcf132e462f2bcd2c0c1133aa6a6a7..130708bea81509321c3cc73f11884af25c95d420 100644
--- a/Mesh/periodical.cpp
+++ b/Plugin/VoroMetal.cpp
@@ -2,82 +2,83 @@
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@onelab.info>.
-//
 // Contributor(s):
-//   Tristan Carrier Maxime Melchior
+//   Tristan Carrier
+//   Maxime Melchior
 
-#include "periodical.h"
-#include "GModel.h"
-#include "meshGRegion.h"
+#include <vector>
 #include <fstream>
 #include <algorithm>
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "GmshMessage.h"
 #include "MElement.h"
-#if defined(HAVE_VORO3D)
-#include "voro++.hh"
-#endif
-
-#if defined(HAVE_VORO3D)
-using namespace voro;
-#endif
-
-/*********definitions*********/
-
-class geo_cell{
- public:
-  std::vector<std::pair<int,int> > lines;
-  std::vector<std::vector<int> > line_loops;
-  std::vector<std::vector<int> > orientations;
-
-  std::vector<int> points2;
-  std::vector<int> lines2;
-  std::vector<int> line_loops2;
-  std::vector<int> faces2;
-  int face_loops2;
+#include "VoroMetal.h"
 
-  geo_cell();
-  ~geo_cell();
+StringXNumber VoroMetalOptions_Number[] = {
+  {GMSH_FULLRC, "ComputeBestSeeds", NULL, 0.},
+  {GMSH_FULLRC, "ComputeMicrostructure", NULL, 1.}
+};
 
-  int search_line(std::pair<int,int>);
+StringXString VoroMetalOptions_String[] = {
+  {GMSH_FULLRC, "SeedsFile", NULL, "seeds.txt"},
 };
 
-/*********class geo_cell*********/
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterVoroMetalPlugin()
+  {
+    return new GMSH_VoroMetalPlugin();
+  }
+}
 
-geo_cell::geo_cell(){}
+std::string GMSH_VoroMetalPlugin::getHelp() const
+{
+  return "Plugin(VoroMetal) creates microstructures using Voronoi diagrams.\n\n";
+}
 
-geo_cell::~geo_cell(){}
+int GMSH_VoroMetalPlugin::getNbOptions() const
+{
+  return sizeof(VoroMetalOptions_Number) / sizeof(StringXNumber);
+}
 
-int geo_cell::search_line(std::pair<int,int> line){
-  unsigned int i;
+int GMSH_VoroMetalPlugin::getNbOptionsStr() const
+{
+  return sizeof(VoroMetalOptions_String) / sizeof(StringXString);
+}
 
-  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;
-  }
+StringXNumber *GMSH_VoroMetalPlugin::getOption(int iopt)
+{
+  return &VoroMetalOptions_Number[iopt];
+}
 
-  return -1;
+StringXString *GMSH_VoroMetalPlugin::getOptionStr(int iopt)
+{
+  return &VoroMetalOptions_String[iopt];
 }
 
-/*********class voroMetal3D*********/
+#if defined(HAVE_MESH) && defined(HAVE_VORO3D)
 
-voroMetal3D::voroMetal3D(){}
+#include "meshGRegion.h"
+#include "voro++.hh"
 
-voroMetal3D::~voroMetal3D(){}
+using namespace voro;
 
-void voroMetal3D::execute(double h){
+void voroMetal3D::execute(double h)
+{
   GRegion* gr;
   GModel* model = GModel::current();
   GModel::riter it;
-
-  for(it=model->firstRegion();it!=model->lastRegion();it++)
-  {
+  for(it = model->firstRegion(); it != model->lastRegion(); it++){
     gr = *it;
-    if(gr->getNumMeshElements()>0){
-      execute(gr,h);
+    if(gr->getNumMeshElements() > 0){
+      execute(gr, h);
     }
   }
 }
 
-void voroMetal3D::execute(GRegion* gr,double h){
+void voroMetal3D::execute(GRegion* gr,double h)
+{
   unsigned int i;
   int j;
   MElement* element;
@@ -91,16 +92,15 @@ void voroMetal3D::execute(GRegion* gr,double h){
   radii.clear();
   vertices.clear();
 
-  for(i=0;i<gr->getNumMeshElements();i++){
+  for(i = 0; i < gr->getNumMeshElements(); i++){
     element = gr->getMeshElement(i);
-
-    for(j=0;j<element->getNumVertices();j++){
+    for(j = 0; j < element->getNumVertices(); j++){
       vertex = element->getVertex(j);
       vertices.insert(vertex);
     }
   }
 
-  for(it=vertices.begin();it!=vertices.end();it++){
+  for(it = vertices.begin(); it != vertices.end(); it++){
     vertices2.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
     radii.push_back(1.0);
   }
@@ -111,24 +111,22 @@ void voroMetal3D::execute(GRegion* gr,double h){
   execute(vertices2,radii,0,h,xMax,yMax,zMax);
 }
 
-void voroMetal3D::execute(std::vector<double>& properties,int radical,double h, double xMax, double yMax, double zMax){
+void voroMetal3D::execute(std::vector<double>& properties, int radical, double h,
+                          double xMax, double yMax, double zMax)
+{
   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]));
+  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,xMax,yMax,zMax);
+  execute(vertices, radii, radical, h, xMax, yMax, zMax);
 }
 
-void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& radii,int radical,double h, double xMax, double yMax, double zMax){
-#if defined(HAVE_VORO3D)
+void voroMetal3D::execute(std::vector<SPoint3>& vertices, std::vector<double>& radii,
+                          int radical, double h, double xMax, double yMax, double zMax)
+{
   unsigned int i;
   unsigned int j;
   unsigned int k;
@@ -178,21 +176,20 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
 
   delta = 0.2*(max_x - min_x);
   min_x=min_y=min_z = 0;
-//  max_x=max_y=max_z = 1;
+  // max_x=max_y=max_z = 1;
   max_x = xMax;
   max_y = yMax;
   max_z = zMax;
   delta = 0;
 
-  container contA(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,true,true,true,vertices.size());
-  //container contA(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,false,false,false,vertices.size());
-  container_poly contB(min_x-delta,max_x+delta,min_y-delta,max_y+delta,min_z-delta,max_z+delta,6,6,6,true,true,true,vertices.size());
-
-  //  wall_cylinder cyl(.5,.5,-3,0,0,6,.4);
-  //  contA.add_wall(cyl);
-
+  container contA(min_x-delta, max_x+delta, min_y-delta, max_y+delta,
+                  min_z-delta, max_z+delta, 6, 6, 6, true, true, true,
+                  vertices.size());
+  container_poly contB(min_x-delta, max_x+delta, min_y-delta, max_y+delta,
+                       min_z-delta, max_z+delta, 6, 6, 6, true, true, true,
+                       vertices.size());
 
-  for(i=0;i<vertices.size();i++){
+  for(i = 0; i < vertices.size(); i++){
     if(radical==0){
       contA.put(i,vertices[i].x(),vertices[i].y(),vertices[i].z());
     }
@@ -206,7 +203,7 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
   c_loop_all loopA(contA);
   c_loop_all loopB(contB);
 
-  if(radical==0){
+  if(radical == 0){
     loopA.start();
     do{
       contA.compute_cell(cell,loopA);
@@ -215,10 +212,9 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
       *pointer = cell;
       pointers.push_back(pointer);
       generators.push_back(SPoint3(x,y,z));
-      //      printf("%d %d (%f,%f,%f)\n",loopA.pid()+1,number+1,x,y,z);
       table.insert(std::pair<int,int>(loopA.pid(),number));
       number++;
-    }while(loopA.inc());
+    } while(loopA.inc());
   }
   else{
     loopB.start();
@@ -229,15 +225,18 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
       *pointer = cell;
       pointers.push_back(pointer);
       generators.push_back(SPoint3(x,y,z));
-      //	  printf("%d %d (%f,%f,%f)\n",loopB.pid()+1,number+1,x,y,z);
       table.insert(std::pair<int,int>(loopB.pid(),number));
       number++;
-    }while(loopB.inc());
+    } while(loopB.inc());
   }
 
   std::ofstream file6("table.txt");
+  if(!file6.is_open()){
+    Msg::Error("Could not open file 'table.txt'");
+    return;
+  }
 
-  for(i=0;i<vertices.size();i++){
+  for(i = 0; i < vertices.size(); i++){
     file6 << i+1 << " " << table[i]+1 << "\n";
   }
 
@@ -245,42 +244,48 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
 
   min_area = 1000000000.0;
 
-  for(i=0;i<pointers.size();i++){
+  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){
+    for(j = 0; j < areas.size(); 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");
+  if(!file.is_open()){
+    Msg::Error("Could not open file 'MicrostructurePolycrystal3D.pos'");
+    return;
+  }
   file << "View \"test\" {\n";
 
   std::ofstream file2("MicrostructurePolycrystal3D.geo");
+  if(!file2.is_open()){
+    Msg::Error("Could not open file 'MicrostructurePolycrystal3D.geo'");
+    return;
+  }
   std::ofstream file5("SET.map");
+  if(!file5.is_open()){
+    Msg::Error("Could not open file 'SET.map'");
+    return;
+  }
   file2 << "c=" << h << ";\n";
   int countPeriodSurf=0;
   int countVolume=0;
-  for(i=0;i<pointers.size();i++){
+  for(i = 0; i < pointers.size(); i++){
     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);
-
+    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());
-
     face_number = 0;
     end = 0;
-    while(end<faces.size()){
+    while(end < faces.size()){
       start = end + 1;
       end = start + faces[end];
       for(j=start;j<end;j++){
@@ -314,15 +319,18 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
         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){
+        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){
+        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){
+        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);
         }
@@ -336,30 +344,32 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
     }
 
     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);
+      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);
+      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++){
+    for(j = 0; j < obj.line_loops.size(); j++){
       temp.clear();
       temp2.clear();
-      for(k=0;k<obj.line_loops[j].size();k++){
+      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);
+      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++){
+    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";
@@ -367,46 +377,51 @@ void voroMetal3D::execute(std::vector<SPoint3>& vertices,std::vector<double>& ra
       increase_counter();
     }
 
-    print_geo_face_loop(get_counter(),obj.faces2,file2);
+    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);
+    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);
+    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];
-#endif
+  for(i = 0; i < pointers.size(); i++) delete pointers[i];
 }
 
-void voroMetal3D::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
+void voroMetal3D::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file)
+{
   file << "SL ("
   << p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
   << p2.x() << ", " << p2.y() << ", " << p2.z()
   << "){10, 20};\n";
 }
 
-void voroMetal3D::initialize_counter(){
+void voroMetal3D::initialize_counter()
+{
   counter = 12;
 }
 
-void voroMetal3D::increase_counter(){
+void voroMetal3D::increase_counter()
+{
   counter = counter+1;
 }
 
-int voroMetal3D::get_counter(){
+int voroMetal3D::get_counter()
+{
   return counter;
 }
 
-void voroMetal3D::print_geo_point(int index,double x,double y,double z,std::ofstream& file){
+void voroMetal3D::print_geo_point(int index,double x,double y,double z,
+                                  std::ofstream& file)
+{
   file.precision(17);
 
   file << "Point(" << index << ")={"
@@ -414,69 +429,74 @@ void voroMetal3D::print_geo_point(int index,double x,double y,double z,std::ofst
   << ",c};\n";
 }
 
-void voroMetal3D::print_geo_line(int index1,int index2,int index3,std::ofstream& file){
+void voroMetal3D::print_geo_line(int index1,int index2,int index3,
+                                 std::ofstream& file)
+{
   file << "Line(" << index1 << ")={"
   << index2 << "," << index3
   << "};\n";
 }
 
-void voroMetal3D::print_geo_face(int index1,int index2,std::ofstream& file){
+void voroMetal3D::print_geo_face(int index1,int index2,std::ofstream& file)
+{
   file << "Plane Surface(" << index1 << ")={"
   << index2
   << "};\n";
 }
 
-void voroMetal3D::print_geo_physical_face(int index1,int index2,std::ofstream& file){
+void voroMetal3D::print_geo_physical_face(int index1,int index2,std::ofstream& file)
+{
   file << "Physical Surface(" << index1 << ")={"
   << index2
   << "};\n";
 }
 
-void voroMetal3D::print_geo_volume(int index1,int index2,std::ofstream& file){
+void voroMetal3D::print_geo_volume(int index1,int index2,std::ofstream& file)
+{
   file << "Volume(" << index1 << ")={"
   << index2
   << "};\n";
 }
 
-void voroMetal3D::print_geo_physical_volume(int index1,int index2,std::ofstream& file){
+void voroMetal3D::print_geo_physical_volume(int index1,int index2,std::ofstream& file)
+{
   file << "Physical Volume(" << index1 << ")={"
   << index2
   << "};\n";
 }
 
-void voroMetal3D::print_geo_line_loop(int index,std::vector<int>& indices,std::vector<int>& orientations,std::ofstream& file){
+void voroMetal3D::print_geo_line_loop(int index,std::vector<int>& indices,
+                                      std::vector<int>& orientations,
+                                      std::ofstream& file)
+{
   unsigned int i;
-
   file << "Line Loop(" << index << ")={";
-
-  for(i=0;i<indices.size();i++){
+  for(i = 0; i < indices.size(); i++){
     if(orientations[i]==1) file << "-";
     file << indices[i];
     if(i<indices.size()-1) file << ",";
   }
-
   file << "};\n";
 }
 
-void voroMetal3D::print_geo_face_loop(int index,std::vector<int>& indices,std::ofstream& file){
+void voroMetal3D::print_geo_face_loop(int index,std::vector<int>& indices,
+                                      std::ofstream& file)
+{
   unsigned int i;
-
   file << "Surface Loop(" << index << ")={";
-
   for(i=0;i<indices.size();i++){
     file << indices[i];
     if(i<indices.size()-1) file << ",";
   }
-
   file << "};\n";
 }
 
-void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax){
+void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax)
+{
   unsigned int i;
   unsigned int j;
   int count;
   int val;
-  //int normal;
   int phase;
   bool flag;
   bool flag1;
@@ -524,8 +544,7 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
 
   faces.clear();
 
-  for(it=model->firstFace();it!=model->lastFace();it++)
-  {
+  for(it = model->firstFace(); it != model->lastFace(); it++){
     gf = *it;
     if(gf->numRegions()==1){
       faces.push_back(gf);
@@ -537,165 +556,182 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
   pairs.clear();
   categories.clear();
 
-  for(i=0;i<faces.size();i++){
+  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 = 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],false));
+  for(i = 0; i < faces.size(); i++){
+    markings.insert(std::pair<GFace*,bool>(faces[i], false));
   }
 
   count = 0;
-
   std::ofstream file("MicrostructurePolycrystal3D.pos");
+  if(!file.is_open()){
+    Msg::Error("Could not open file 'MicrostructurePolycrystal3D.pos'");
+    return;
+  }
   file << "View \"test\" {\n";
 
   std::ofstream file2("PERIODIC.map");
+  if(!file2.is_open()){
+    Msg::Error("Could not open file 'PERIODIC.map'");
+    return;
+  }
 
-  for(i=0;i<faces.size();i++){
-    for(j=0;j<faces.size();j++){
+  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,xMax,yMax,zMax);
-
       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 (fabs((p2.x()-p1.x()-1.0))<0.0001){
-            if (fabs((p2.y()-p1.y()))<0.0001){
-              if (fabs((p2.z()-p1.z()))<0.0001){
+          if (fabs((p2.x()-p1.x()-1.0)) < 0.0001){
+            if (fabs((p2.y()-p1.y())) < 0.0001){
+              if (fabs((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(fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if(fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p2.y()-p1.y()-1.0))<0.0001){
+            }
+            else if (fabs((p2.y()-p1.y()-1.0))<0.0001){
               if (fabs((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 (fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((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 (fabs((p1.y()-p2.y()-1.0))<0.0001){
+            }
+            else if (fabs((p1.y()-p2.y()-1.0))<0.0001){
               if (fabs((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 (fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.x()-p2.x()-1.0))<0.0001){
+          }
+          else if (fabs((p1.x()-p2.x()-1.0))<0.0001){
             if (fabs((p2.y()-p1.y()))<0.0001){
               if (fabs((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(fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if(fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p2.y()-p1.y()-1.0))<0.0001){
+            }
+            else if (fabs((p2.y()-p1.y()-1.0))<0.0001){
               if (fabs((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 (fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.y()-p2.y()-1.0))<0.0001){
+            }
+            else if (fabs((p1.y()-p2.y()-1.0))<0.0001){
               if (fabs((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 (fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.x()-p2.x()))<0.0001){
+          }
+          else if (fabs((p1.x()-p2.x()))<0.0001){
             if (fabs((p2.y()-p1.y()-1.0))<0.0001){
               if (fabs((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 (fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.y()-p2.y()-1.0))<0.0001){
+            }
+            else if (fabs((p1.y()-p2.y()-1.0))<0.0001){
               if (fabs((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 (fabs((p2.z()-p1.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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 (fabs((p1.y()-p2.y()))<0.0001){
+            }
+            else if (fabs((p1.y()-p2.y()))<0.0001){
               if (fabs((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 (fabs((p1.z()-p2.z()-1.0))<0.0001){
+              }
+              else if (fabs((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";
               }
@@ -709,21 +745,25 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
 
   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);
+  if(!file3.is_open()){
+    Msg::Error("Could not open file 'MicrostructurePolycrystal3D.geo'");
+    return;
+  }
   file3.precision(17);
 
   std::ofstream file4("MicrostructurePolycrystal3D2.pos");
+  if(!file4.is_open()){
+    Msg::Error("Could not open file 'MicrostructurePolycrystal3D2.pos'");
+    return;
+  }
   file4 << "View \"test\" {\n";
 
   file3 << "Physical Surface(11)={";
 
   count = 0;
-  for(it=model->firstFace();it!=model->lastFace();it++)
-  {
+  for(it=model->firstFace();it!=model->lastFace();it++){
     gf = *it;
     if(count>0) file3 << ",";
     file3 << gf->tag();
@@ -732,39 +772,28 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
 
   file3 << "};\n";
 
-  for(i=0;i<pairs.size();i++){
+  for(i = 0; i < pairs.size(); i++){
     gf1 = pairs[i].first;
     gf2 = pairs[i].second;
-
-    std::list<GVertex*> gv1 = gf1->vertices(); 
-    std::list<GVertex*> gv2 = gf2->vertices(); 
-
+    std::list<GVertex*> gv1 = gf1->vertices();
+    std::list<GVertex*> gv2 = gf2->vertices();
     std::list<GVertex*>::iterator it1 = gv1.begin();
     std::list<GVertex*>::iterator it2 = gv2.begin();
-
     SPoint3 cg1 (0,0,0);
     SPoint3 cg2 (0,0,0);
-    
     for (; it1 != gv1.end(); it1++,it2++){
       cg1 += SPoint3((*it1)->x(),(*it1)->y(),(*it1)->z());
       cg2 += SPoint3((*it2)->x(),(*it2)->y(),(*it2)->z());
     }
-
     SVector3 dx = (cg2-cg1) * (1./gv1.size());
-    
-    //    printf("%g %g %g\n",dx.x(),dx.y(),dx.z());
-        
     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();
@@ -775,31 +804,30 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
         indices3.push_back(-(*it8)->tag());
     }
     int countReverseEdge=0;
-    for(it7=edges1.begin();it7!=edges1.end();it7++,it9++){
+    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,xMax,yMax,zMax);
-        correspondance(fabs(v4->x()-v2->x()),fabs(v4->y()-v2->y()),fabs(v4->z()-v2->z()),e,categories[i],flag2,xMax,yMax,zMax);
-
-        correspondance(fabs(v4->x()-v1->x()),fabs(v4->y()-v1->y()),fabs(v4->z()-v1->z()),e,categories[i],flag3,xMax,yMax,zMax);
-        correspondance(fabs(v3->x()-v2->x()),fabs(v3->y()-v2->y()),fabs(v3->z()-v2->z()),e,categories[i],flag4,xMax,yMax,zMax);
-
+        correspondance(fabs(v3->x()-v1->x()),fabs(v3->y()-v1->y()),fabs(v3->z()-v1->z()),
+                       e,categories[i],flag1,xMax,yMax,zMax);
+        correspondance(fabs(v4->x()-v2->x()),fabs(v4->y()-v2->y()),fabs(v4->z()-v2->z()),
+                       e,categories[i],flag2,xMax,yMax,zMax);
+
+        correspondance(fabs(v4->x()-v1->x()),fabs(v4->y()-v1->y()),fabs(v4->z()-v1->z()),
+                       e,categories[i],flag3,xMax,yMax,zMax);
+        correspondance(fabs(v3->x()-v2->x()),fabs(v3->y()-v2->y()),fabs(v3->z()-v2->z()),
+                       e,categories[i],flag4,xMax,yMax,zMax);
         if(flag1 && flag2){
           if(phase==1){
             mem = it8;
@@ -807,14 +835,12 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
           }
           else if(phase==2){
             mem++;
-
             /*if(it8==mem){
               normal = 1;
             }
             else{
               normal = -1;
             }*/
-
             phase = 3;
           }
           if (*it9==1)
@@ -833,14 +859,12 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
           }
           else if(phase==2){
             mem++;
-
             /*if(it8==mem){
               normal = 1;
             }
             else{
               normal = -1;
             }*/
-
             phase = 3;
           }
           if (*it9==1)
@@ -858,33 +882,15 @@ void voroMetal3D::correspondance(double e, double xMax, double yMax, double zMax
     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++){
-      if(j>0) file3 << ",";
-      file3 << indices2[j];
-    }
-
-    file3 << "};\n";
-    */
-
-    
-    file3 << "Periodic Surface {" << gf1->tag() << " }={ " << gf2->tag() << " } Translate { " << -dx.x() << "," << -dx.y() << "," << -dx.z() << "};\n"; 
-
-    
+    file3 << "Periodic Surface {" << gf1->tag() << " }={ " << gf2->tag()
+          << " } Translate { " << -dx.x() << "," << -dx.y() << "," << -dx.z() << "};\n";
   }
-
   file4 << "};\n";
 }
 
-bool voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,double e,int& val,double xMax,double yMax,double zMax){
+bool voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,double e,
+                                 int& val,double xMax,double yMax,double zMax)
+{
   bool flag;
 
   flag = 0;
@@ -920,11 +926,12 @@ bool voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,do
     flag = 1;
     val = 7;
   }
-
   return flag;
 }
 
-void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,double e,int val,bool& flag,double xMax,double yMax,double zMax){
+void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,
+                                 double e,int val,bool& flag,double xMax,double yMax,double zMax)
+{
   flag = 0;
 
   if(val==1 && equal(delta_x,xMax,e) && equal(delta_y,0.0,e) && equal(delta_z,0.0,e)){
@@ -952,21 +959,19 @@ void voroMetal3D::correspondance(double delta_x,double delta_y,double delta_z,do
   }
 }
 
-bool voroMetal3D::equal(double x,double y,double e){
+bool voroMetal3D::equal(double x,double y,double e)
+{
   bool flag;
-
   if(x>y-e && x<y+e){
     flag = 1;
   }
   else{
     flag = 0;
   }
-
   return flag;
 }
 
-
-void microstructure(const char *filename)
+static void microstructure(const char *filename)
 {
   int j;
   int radical;
@@ -977,6 +982,10 @@ void microstructure(const char *filename)
   std::vector<double> properties;
   if(filename){
     std::ifstream file(filename);
+    if(!file.is_open()){
+      Msg::Error("Could not open file '%s'", filename);
+      return;
+    }
     file >> max;
     file >> radical;
     file >> xMax;
@@ -998,7 +1007,7 @@ void microstructure(const char *filename)
   }
 }
 
-void computeBestSeeds(const char *filename)
+static void computeBestSeeds(const char *filename)
 {
   int j;
   int radical;
@@ -1010,6 +1019,10 @@ void computeBestSeeds(const char *filename)
   std::cout<<"entree dans computeBestSeeds"<<std::endl;
   if(filename){
     std::ifstream file(filename);
+    if(!file.is_open()){
+      Msg::Error("Could not open file '%s'", filename);
+      return;
+    }
     file >> max;
     file >> radical;
     file >> xMax;
@@ -1059,7 +1072,9 @@ void computeBestSeeds(const char *filename)
           GEdge* eTmp = (*ite);
           GVertex* vTmp1 = eTmp->getBeginVertex();
           GVertex* vTmp2 = eTmp->getEndVertex();
-          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) +
+                                (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) +
+                                (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
           if (distTmp < distMinTmp){
             distMinTmp = distTmp;
           }
@@ -1072,7 +1087,7 @@ void computeBestSeeds(const char *filename)
         }
         delete m;
       }
-      for(j=0;j<max;j++){
+      for(j = 0; j < max; j++){
         std::cout<<"j "<<j<<std::endl;
         std::vector<double> propertiesModified;
         propertiesModified.clear();
@@ -1097,12 +1112,14 @@ void computeBestSeeds(const char *filename)
           GEdge* eTmp = (*ite);
           GVertex* vTmp1 = eTmp->getBeginVertex();
           GVertex* vTmp2 = eTmp->getEndVertex();
-          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) +
+                                (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) +
+                                (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
           if (distTmp < distMinTmp){
             distMinTmp = distTmp;
           }
         }
-        if (distMinTmp > distMinGlobal){
+        if(distMinTmp > distMinGlobal){
           distMinGlobal = distMinTmp;
           jMinGlobal = j;
           xORyORz = 2;
@@ -1110,7 +1127,7 @@ void computeBestSeeds(const char *filename)
         }
         delete m;
       }
-      for(j=0;j<max;j++){
+      for(j = 0; j < max; j++){
         std::cout<<"j "<<j<<std::endl;
         std::vector<double> propertiesModified;
         propertiesModified.clear();
@@ -1135,7 +1152,9 @@ void computeBestSeeds(const char *filename)
           GEdge* eTmp = (*ite);
           GVertex* vTmp1 = eTmp->getBeginVertex();
           GVertex* vTmp2 = eTmp->getEndVertex();
-          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) +
+                                (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) +
+                                (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
           if (distTmp < distMinTmp){
             distMinTmp = distTmp;
           }
@@ -1148,7 +1167,7 @@ void computeBestSeeds(const char *filename)
         }
         delete m;
       }
-      for(j=0;j<max;j++){
+      for(j = 0; j < max; j++){
         std::cout<<"j "<<j<<std::endl;
         std::vector<double> propertiesModified;
         propertiesModified.clear();
@@ -1173,7 +1192,9 @@ void computeBestSeeds(const char *filename)
           GEdge* eTmp = (*ite);
           GVertex* vTmp1 = eTmp->getBeginVertex();
           GVertex* vTmp2 = eTmp->getEndVertex();
-          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) +
+                                (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) +
+                                (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
           if (distTmp < distMinTmp){
             distMinTmp = distTmp;
           }
@@ -1211,7 +1232,9 @@ void computeBestSeeds(const char *filename)
           GEdge* eTmp = (*ite);
           GVertex* vTmp1 = eTmp->getBeginVertex();
           GVertex* vTmp2 = eTmp->getEndVertex();
-          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) +
+                                (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) +
+                                (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
           if (distTmp < distMinTmp){
             distMinTmp = distTmp;
           }
@@ -1249,7 +1272,9 @@ void computeBestSeeds(const char *filename)
           GEdge* eTmp = (*ite);
           GVertex* vTmp1 = eTmp->getBeginVertex();
           GVertex* vTmp2 = eTmp->getEndVertex();
-          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) + (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) + (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
+          double distTmp = sqrt((vTmp1->x() -  vTmp2->x()) * (vTmp1->x() -  vTmp2->x()) +
+                                (vTmp1->y() -  vTmp2->y()) * (vTmp1->y() -  vTmp2->y()) +
+                                (vTmp1->z() -  vTmp2->z()) * (vTmp1->z() -  vTmp2->z()));
           if (distTmp < distMinTmp){
             distMinTmp = distTmp;
           }
@@ -1294,8 +1319,28 @@ void computeBestSeeds(const char *filename)
     }
     std::cout<<"liste des nouveaux seeds :"<<std::endl;
     for(unsigned int iTmp = 0; iTmp < max;iTmp++){
-      std::cout<<properties[4*iTmp]<<" "<<properties[4*iTmp + 1]<<" "<<properties[4*iTmp + 2]<<" "<<properties[4*iTmp + 3]<<std::endl;
-
+      std::cout<<properties[4*iTmp]<<" "<<properties[4*iTmp + 1]<<" "
+               <<properties[4*iTmp + 2]<<" "<<properties[4*iTmp + 3]<<std::endl;
     }
   }
 }
+
+PView *GMSH_VoroMetalPlugin::execute(PView *v)
+{
+  int runBestSeeds = (int)VoroMetalOptions_Number[0].def;
+  int runMicrostructure = (int)VoroMetalOptions_Number[1].def;
+  std::string seedsFile = VoroMetalOptions_String[0].def;
+  if(runBestSeeds) computeBestSeeds(seedsFile.c_str());
+  if(runMicrostructure) microstructure(seedsFile.c_str());
+  return v;
+}
+
+#else
+
+PView *GMSH_VoroMetalPlugin::execute(PView *v)
+{
+  Msg::Error("Plugin(VoroMetal) requires mesh module and voro++");
+  return v;
+}
+
+#endif
diff --git a/Mesh/periodical.h b/Plugin/VoroMetal.h
similarity index 50%
rename from Mesh/periodical.h
rename to Plugin/VoroMetal.h
index 7ff7377a91588783b5f3ed4de16fd575735667b2..71104f099f215ed4ac3d5ceb2da28ff2aef35fc3 100644
--- a/Mesh/periodical.h
+++ b/Plugin/VoroMetal.h
@@ -2,18 +2,42 @@
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to the public mailing list <gmsh@onelab.info>.
-//
-// Contributor(s):
-//   Tristan Carrier Maxime Melchior
 
-#include "GRegion.h"
+#ifndef _VOROMETAL_H_
+#define _VOROMETAL_H_
+
+#include <vector>
+#include "Plugin.h"
+
+class geo_cell{
+ public:
+  std::vector<std::pair<int,int> > lines;
+  std::vector<std::vector<int> > line_loops;
+  std::vector<std::vector<int> > orientations;
+  std::vector<int> points2;
+  std::vector<int> lines2;
+  std::vector<int> line_loops2;
+  std::vector<int> faces2;
+  int face_loops2;
+  geo_cell(){}
+  ~geo_cell(){}
+  int search_line(std::pair<int,int> line)
+  {
+    unsigned int i;
+    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;
+    }
+    return -1;
+  }
+};
 
 class voroMetal3D{
  private:
   int counter;
  public:
-  voroMetal3D();
-  ~voroMetal3D();
+  voroMetal3D(){}
+  ~voroMetal3D(){}
   void execute(double);
   void execute(GRegion*,double);
   void execute(std::vector<SPoint3>&,std::vector<double>&,int,double,double,double,double);
@@ -36,8 +60,27 @@ class voroMetal3D{
   bool equal(double,double,double);
 };
 
+extern "C"
+{
+  GMSH_Plugin *GMSH_RegisterVoroMetalPlugin();
+}
 
-void microstructure(const char *filename);
-void computeBestSeeds(const char *filename);
-
+class GMSH_VoroMetalPlugin : public GMSH_PostPlugin
+{
+ public:
+  GMSH_VoroMetalPlugin(){}
+  std::string getName() const { return "VoroMetal"; }
+  std::string getShortHelp() const
+  {
+    return "Voronoi microstructures";
+  }
+  std::string getHelp() const;
+  std::string getAuthor() const { return "Tristan Carrier & Maxime Melchior"; }
+  int getNbOptions() const;
+  StringXNumber* getOption(int iopt);
+  int getNbOptionsStr() const;
+  StringXString *getOptionStr(int iopt);
+  PView *execute(PView *);
+};
 
+#endif
diff --git a/benchmarks/misc/seeds64doubletw.txt b/benchmarks/misc/seeds64doubletw.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cefa5fc6bede7f46cdd15eb94e1b2dd2ed37d956
--- /dev/null
+++ b/benchmarks/misc/seeds64doubletw.txt
@@ -0,0 +1,69 @@
+   68    0    1.    1.    1.

+0.5602  0.5768  0.4357         1

+0.5713  0.5688  0.4418         1

+0.5491  0.5848  0.4296         1

+0.5379  0.6113  0.6941         1

+0.5261  0.6204  0.6959         1

+0.5497  0.6022  0.6923         1

+0.6873  0.7394  0.5293         1

+0.6133  0.3627  0.5324         1

+0.1458  0.1328  0.3812         1

+0.7003  0.5497  0.8256         1

+0.5373  0.6623  0.0532         1

+0.1715  0.6187  0.7587         1

+0.9101  0.7140  0.7400         1

+0.3567  0.8219  0.5521         1

+0.6508  0.1876  0.8183         1

+0.0224  0.6572  0.1196         1

+0.3854  0.2426  0.9653         1

+0.4853  0.9923  0.2813         1

+0.9784  0.9529  0.2431         1

+0.1387  0.4549  0.9844         1

+0.4249  0.7048  0.3654         1

+0.2004  0.4096  0.2392         1

+0.7050  0.0791  0.2038         1

+0.2546  0.8227  0.2816         1

+0.6248  0.7458  0.8739         1

+0.0452  0.1465  0.8069         1

+0.0928  0.6517  0.3952         1

+0.6199  0.8252  0.3391         1

+0.2802  0.9931  0.7612         1

+0.5263  0.9937  0.7263         1

+0.8715  0.5572  0.3972         1

+0.5489  0.3050  0.1544         1

+0.0010  0.4132  0.8038         1

+0.8779  0.1408  0.0431         1

+0.6007  0.0963  0.5492         1

+0.3143  0.0031  0.9996         1

+0.8804  0.2438  0.5439         1

+0.9771  0.0028  0.4942         1

+0.9610  0.2860  0.2800         1

+0.1864  0.0621  0.1669         1

+0.3355  0.0698  0.4926         1

+0.8340  0.9321  0.7665         1

+0.4538  0.2884  0.7151         1

+0.1844  0.4269  0.6276         1

+0.1173  0.8548  0.5805         1

+0.7669  0.8185  0.0563         1

+0.7686  0.3834  0.0918         1

+0.7193  0.9982  0.9369         1

+0.4796  0.1995  0.3542         1

+0.9097  0.7816  0.4809         1

+0.6418  0.6288  0.2492         1

+0.3331  0.5907  0.0614         1

+0.5100  0.8816  0.0027         1

+0.1555  0.7684  0.9837         1

+0.7948  0.5999  0.0487         1

+0.3215  0.3967  0.4377         1

+0.4659  0.4677  0.9396         1

+0.2469  0.2219  0.7829         1

+0.0729  0.3849  0.4411         1

+0.9547  0.5321  0.6121         1

+0.3699  0.7445  0.8400         1

+0.9925  0.9447  0.9161         1

+0.7929  0.4079  0.6849         1

+0.1046  0.2266  0.0321         1

+0.7516  0.2727  0.3346         1

+0.1198  0.2221  0.6062         1

+0.3474  0.2320  0.1859         1

+0.5263  0.1019  0.0603         1

diff --git a/benchmarks/misc/vorometal.geo b/benchmarks/misc/vorometal.geo
new file mode 100644
index 0000000000000000000000000000000000000000..a48e5916479151126f74f9789dff6246f0da447d
--- /dev/null
+++ b/benchmarks/misc/vorometal.geo
@@ -0,0 +1,4 @@
+Plugin(VoroMetal).SeedsFile = "seeds64doubletw.txt";
+Plugin(VoroMetal).ComputeBestSeeds = 0;
+Plugin(VoroMetal).ComputeMicrostructure = 1;
+Plugin(VoroMetal).Run;