diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..2b8dc3f
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+build
+sources/libmeshb7.mod
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bc087c1..94b6e28 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,12 +9,13 @@ set (libMeshb_VERSION_MAJOR 7)
 set (libMeshb_VERSION_MINOR 54)
 project(libMeshb VERSION ${libMeshb_VERSION_MAJOR}.${libMeshb_VERSION_MINOR} LANGUAGES C)
 
-option(WITH_AIO   "Use Unix low-level and asynchronous I/O for higher speed"       OFF)
-option(WITH_CPACK "Enable cpack target to generate a zip file containing binaries" OFF)
+option(WITH_GMF_AIO     "Use Unix low-level and asynchronous I/O for higher speed"       OFF)
+option(WITH_GMF_CPACK   "Enable cpack target to generate a zip file containing binaries" OFF)
+option(WITH_GMF_FORTRAN "Build the Fortran API"                                          ON )
 
 include (CheckLanguage)
 check_language (Fortran)
-if(CMAKE_Fortran_COMPILER)
+if(CMAKE_Fortran_COMPILER AND WITH_GMF_FORTRAN)
   enable_language(Fortran)
   set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules)
   if (CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
@@ -34,14 +35,18 @@ if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
    if (WIN32)
       set(CMAKE_INSTALL_PREFIX "$ENV{HOMEPATH}/cmakebuilds" CACHE PATH "..." FORCE)
    else ()
-      set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/cmakebuilds/$ENV{OSTYPE}-$ENV{MACHTYPE}" CACHE PATH "..." FORCE)
+      if (DEFINED ENV{OSTYPE} AND DEFINED ENV{MACHTYPE})
+         set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/cmakebuilds/$ENV{OSTYPE}-$ENV{MACHTYPE}" CACHE PATH "..." FORCE)
+      else ()
+         set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/cmakebuilds" CACHE PATH "..." FORCE)
+      endif ()
    endif()
 endif ()
 
 list (APPEND CMAKE_PREFIX_PATH "${CMAKE_INSTALL_PREFIX}")
 
-if(WITH_AIO)
-   set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DWITH_AIO")
+if(WITH_GMF_AIO)
+   set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DWITH_GMF_AIO")
    IF(${CMAKE_SYSTEM_NAME} MATCHES Linux)
       set(AIO_LIBRARIES rt)
    endif()
@@ -66,7 +71,7 @@ install (DIRECTORY sample_meshes DESTINATION share/libMeshb)
 # SET PACKAGE AND DEPLOYMENT VARIABLES
 ######################################
 
-if (WITH_CPACK)
+if (WITH_GMF_CPACK)
    include (InstallRequiredSystemLibraries)
    set (CPACK_GENERATOR TXZ)
    set (CPACK_RESOURCE_FILE_LICENSE ${CMAKE_CURRENT_SOURCE_DIR}/copyright.txt)
@@ -106,6 +111,7 @@ install(FILES
 ##################################
 
 message("-- Build mode            : " ${CMAKE_BUILD_TYPE})
-message("-- cpack target enabled  : " ${WITH_CPACK})
-message("-- Asynchronous IO       : " ${WITH_AIO})
+message("-- cpack target enabled  : " ${WITH_GMF_CPACK})
+message("-- Asynchronous IO       : " ${WITH_GMF_AIO})
+message("-- Fortran API           : " ${WITH_GMF_FORTRAN})
 message("-- Install directory     : " ${CMAKE_INSTALL_PREFIX})
diff --git a/Documentation/libMeshb7.tex b/Documentation/libMeshb7.tex
index 3f44a53..99bd6eb 100644
--- a/Documentation/libMeshb7.tex
+++ b/Documentation/libMeshb7.tex
@@ -467,7 +467,7 @@ \subsubsection{Writing mode}
 \normalfont
 
 \paragraph{Asynchronous Input Output:}
-to get best performance out of flash storage SSD (more than a GB/s), the library needs to access the filesystem through low-level functions and perform the reading and processing tasks in parallel using asynchronous I/O. To do so, you need to compile the library with the {\tt -DWITH\_AIO} option and link the final executable with the {\tt -lrt} library under Linux.
+to get best performance out of flash storage SSD (more than a GB/s), the library needs to access the filesystem through low-level functions and perform the reading and processing tasks in parallel using asynchronous I/O. To do so, you need to compile the library with the {\tt -DWITH\_GMF\_AIO} option and link the final executable with the {\tt -lrt} library under Linux.
 
 
 \subsection{GmfCloseMesh}
diff --git a/LICENSE.txt b/LICENSE.txt
index 1e5e8c0..6df08a8 100644
--- a/LICENSE.txt
+++ b/LICENSE.txt
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2001-2018 Loïc Maréchal / INRIA
+Copyright (c) 2001-2022 Loïc Maréchal / INRIA
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/README.md b/README.md
index 6b68bdc..6f63254 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
 
-## libMeshb version 7.60
+## libMeshb version 7.62
 A library to handle the *.meshb file format.
 
 ## Overview
diff --git a/copyright.txt b/copyright.txt
index 8fdc0e3..560ca26 100644
--- a/copyright.txt
+++ b/copyright.txt
@@ -1,4 +1,4 @@
-All libMeshb code is Copyright 2001-2018 - by Loïc Maréchal / INRIA.
+All libMeshb code is Copyright 2001-2022 - by Loïc Maréchal / INRIA.
 This program is a free software.
 You can redistribute it and/or modify it under the terms of the MIT License
 as published by the Open Source Initiative.
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 3ab8a55..7b38774 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -37,7 +37,7 @@ install (TARGETS mesh2poly DESTINATION share/libMeshb/examples COMPONENT example
 ########################
 
 
-if (CMAKE_Fortran_COMPILER)
+if (CMAKE_Fortran_COMPILER AND WITH_GMF_FORTRAN)
    add_executable(test_libmeshb_f77 test_libmeshb.f)
    target_link_libraries(test_libmeshb_f77 Meshbf.7 ${AIO_LIBRARIES})
    install (TARGETS test_libmeshb_f77 DESTINATION share/libMeshb/examples COMPONENT examples)
diff --git a/examples/mesh2poly.c b/examples/mesh2poly.c
index 5992f74..fe1a5c3 100644
--- a/examples/mesh2poly.c
+++ b/examples/mesh2poly.c
@@ -1,11 +1,38 @@
+
+
+/*----------------------------------------------------------------------------*/
+/* Includes                                                                   */
+/*----------------------------------------------------------------------------*/
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <libmeshb7.h>
 
+
+/*----------------------------------------------------------------------------*/
+/* Topologic tables                                                           */
+/*----------------------------------------------------------------------------*/
+
+int tettvpe[6][2] = { {0,1}, {1,2}, {2,0}, {3,0}, {3,1}, {3,2} };
+int tettvpf[4][3] = { {1,2,3}, {2,0,3}, {3,0,1}, {0,2,1} };
+int pyrtvpe[8][2] = { {0,1}, {3,2}, {0,3}, {1,2}, {0,4}, {1,4}, {2,4}, {3,4} };
+int pyrtvpf[5][4] = { {0,1,4,-1}, {1,2,4,-1}, {2,3,4,-1}, {3,0,4,-1}, {3,2,1,0} };
+int pritvpe[9][2] = { {0,1}, {1,2}, {2,0}, {3,4}, {4,5}, {5,3}, {0,3}, {1,4}, {2,5} };
+int pritvpf[5][4] = { {3,5,2,0}, {1,2,5,4}, {0,1,4,3}, {2,1,0,-1}, {3,4,5,-1} };
+int hextvpe[12][2] = {  {3,2}, {0,1}, {4,5}, {7,6}, {3,7}, {2,6},
+                        {1,5}, {0,4}, {3,0}, {7,4}, {6,5}, {2,1} };
+int hextvpf[6][4] = {   {3,0,4,7}, {6,5,1,2}, {3,2,1,0},
+                        {4,5,6,7},{3,7,6,2}, {1,5,4,0} };
+
+
+/*----------------------------------------------------------------------------*/
+/* Transform a surface mesh into polygons and a tet mesh into polyhedra       */ 
+/*----------------------------------------------------------------------------*/
+
 int main(int ArgCnt, char **ArgVec)
 {
-   int      i, j, ver, dim, ref, nod[8], ArgIdx = 1;
+   int      i, j, k, ver, dim, ref, nod[8], ArgIdx = 1;
    int      NmbVer, NmbTri, NmbQad, NmbTet, NmbPyr, NmbPri, NmbHex;
    int64_t  InpMsh, OutMsh;
    float    xf, yf, zf;
@@ -14,7 +41,7 @@ int main(int ArgCnt, char **ArgVec)
 
    if(ArgCnt != 3)
    {
-      puts("\nMESH2POLY 1.0, may 12 2021, Loic MARECHAL / INRIA\n");
+      puts("\nMESH2POLY 1.1, september 15 2021, Loic MARECHAL / INRIA\n");
       puts(" Usage    : mesh2poly source_name destination_name\n");
       exit(0);
    }
@@ -62,6 +89,7 @@ int main(int ArgCnt, char **ArgVec)
 
    printf(" InpMsh: idx = %lld, version: %d, dimension: %d\n", InpMsh, ver, dim);
 
+   // Copy all the vertices
    if(NmbVer)
    {
       printf(" Copying %d vertices\n", NmbVer);
@@ -84,6 +112,7 @@ int main(int ArgCnt, char **ArgVec)
       }
    }
 
+   // Transform triangles and quads into boundary polygons
    if(NmbTri || NmbQad)
    {
       GmfSetKwd(OutMsh, GmfBoundaryPolygonHeaders, NmbTri + NmbQad);
@@ -91,25 +120,17 @@ int main(int ArgCnt, char **ArgVec)
       if(NmbTri)
       {
          printf(" Adding  %d triangles to the polygons headers\n", NmbTri);
-         GmfGotoKwd(InpMsh, GmfTriangles);
 
          for(i=1;i<=NmbTri;i++)
-         {
-            GmfGetLin(InpMsh, GmfTriangles, &nod[0], &nod[1], &nod[2], &ref);
             GmfSetLin(OutMsh, GmfBoundaryPolygonHeaders, (i-1) * 3 + 1, ref);
-         }
       }
 
       if(NmbQad)
       {
          printf(" Adding  %d quads to the polygons headers\n", NmbQad);
-         GmfGotoKwd(InpMsh, GmfQuadrilaterals);
 
          for(i=1;i<=NmbQad;i++)
-         {
-            GmfGetLin(InpMsh, GmfQuadrilaterals, &nod[0], &nod[1], &nod[2], &nod[3], &ref);
             GmfSetLin(OutMsh, GmfBoundaryPolygonHeaders, (i-1) * 4 + NmbTri * 3 + 1, ref);
-         }
       }
 
       GmfSetKwd(OutMsh, GmfBoundaryPolygonVertices, NmbTri * 3 + NmbQad * 4);
@@ -143,6 +164,63 @@ int main(int ArgCnt, char **ArgVec)
       }
    }
 
+   // Transform tets into inner polygons and polyhedra
+   if(NmbTet)
+   {
+      GmfSetKwd(OutMsh, GmfInnerPolygonHeaders, 4 * NmbTet);
+
+      if(NmbTet)
+      {
+         printf(" Adding  %d tet's faces to the polygons headers\n", 4 * NmbTet);
+
+         for(i=0;i<NmbTet;i++)
+            for(j=0;j<4;j++)
+               GmfSetLin(OutMsh, GmfInnerPolygonHeaders, i * 12 + j * 3 + 1, ref);
+      }
+
+      GmfSetKwd(OutMsh, GmfInnerPolygonVertices, 12 * NmbTet);
+
+      if(NmbTet)
+      {
+         printf(" Adding  %d tet's faces to the polygons faces\n", NmbTet);
+         GmfGotoKwd(InpMsh, GmfTetrahedra);
+
+         for(i=1;i<=NmbTet;i++)
+         {
+            GmfGetLin(InpMsh, GmfTetrahedra, &nod[0], &nod[1], &nod[2], &nod[3], &ref);
+
+            for(j=0;j<4;j++)
+               for(k=0;k<3;k++)
+                  GmfSetLin(OutMsh, GmfInnerPolygonVertices, nod[ tettvpf[j][k] ]);
+         }
+      }
+
+      GmfSetKwd(OutMsh, GmfPolyhedraHeaders, NmbTet);
+
+      if(NmbTet)
+      {
+         printf(" Adding  %d tets to the polyhedra headers\n", NmbTet);
+         GmfGotoKwd(InpMsh, GmfTetrahedra);
+
+         for(i=0;i<NmbTet;i++)
+         {
+            GmfGetLin(InpMsh, GmfTetrahedra, &nod[0], &nod[1], &nod[2], &nod[3], &ref);
+            GmfSetLin(OutMsh, GmfPolyhedraHeaders, i * 4 + 1, ref);
+         }
+      }
+
+      GmfSetKwd(OutMsh, GmfPolyhedraFaces, NmbTet * 4);
+
+      if(NmbTet)
+      {
+         printf(" Adding  %d tet's face to the polyhedra faces\n", NmbTet * 4);
+
+         for(i=0;i<NmbTet;i++)
+            for(j=0;j<4;j++)
+               GmfSetLin(OutMsh, GmfPolyhedraFaces, i * 4 + j + 1);
+      }
+   }
+
    GmfCloseMesh(InpMsh);
    GmfCloseMesh(OutMsh);
 
diff --git a/examples/test_libmeshb_block_pipeline.f b/examples/test_libmeshb_block_pipeline.f
index 940d634..e54a85c 100644
--- a/examples/test_libmeshb_block_pipeline.f
+++ b/examples/test_libmeshb_block_pipeline.f
@@ -1,5 +1,5 @@
 
-c     libmeshb example : transform a quadrilateral mesh into a triangular one
+c     libmeshb example: transform a quadrilateral mesh into a triangular one
 c     using fast block transfer and pipelined post processing
 
       include 'libmeshb7.ins'
@@ -39,14 +39,15 @@
 
 c     Read the vertices
       res = gmfgetblock(InpMsh, GmfVertices, 1_8, NmbVer,
-     +                  movver, 1, VerTab,
+     +                  0, %val(0), movver, 1, VerTab,
      +                  GmfDouble, VerTab(1,1), VerTab(1,2),
      +                  GmfDouble, VerTab(2,1), VerTab(2,2),
      +                  GmfDouble, VerTab(3,1), VerTab(3,2),
      +                  GmfInt,    RefTab(1),   RefTab(2))
 
 c     Read the quads
-      res = gmfgetblock(InpMsh, GmfQuadrilaterals, 1_8, NmbQad,0,
+      res = gmfgetblock(InpMsh, GmfQuadrilaterals, 1_8, NmbQad,
+     +                  0, %val(0), %val(0),
      +                  GmfInt, QadTab(1,1), QadTab(1,2),
      +                  GmfInt, QadTab(2,1), QadTab(2,2),
      +                  GmfInt, QadTab(3,1), QadTab(3,2),
@@ -79,7 +80,7 @@
 c     Write the triangles
       res = gmfsetkwd(OutMsh, GmfTriangles, 2*NmbQad, 0, 0)
       res = gmfsetblock(OutMsh, GmfTriangles, 1_8, 2*NmbQad,
-     +                  0,%val(0),%val(0),
+     +                  0,%val(0),
      +                  qad2tri, 2, QadTab, TriTab,
      +                  GmfInt, TriTab(1,1), TriTab(1,2),
      +                  GmfInt, TriTab(2,1), TriTab(2,2),
diff --git a/utilities/bb2sol.c b/legacy_sources/utilities/bb2sol.c
similarity index 100%
rename from utilities/bb2sol.c
rename to legacy_sources/utilities/bb2sol.c
diff --git a/utilities/mesh2patran.c b/legacy_sources/utilities/mesh2patran.c
similarity index 100%
rename from utilities/mesh2patran.c
rename to legacy_sources/utilities/mesh2patran.c
diff --git a/utilities/noboiteb2meshb.c b/legacy_sources/utilities/noboiteb2meshb.c
similarity index 100%
rename from utilities/noboiteb2meshb.c
rename to legacy_sources/utilities/noboiteb2meshb.c
diff --git a/utilities/sol2bb.c b/legacy_sources/utilities/sol2bb.c
similarity index 100%
rename from utilities/sol2bb.c
rename to legacy_sources/utilities/sol2bb.c
diff --git a/sources/CMakeLists.txt b/sources/CMakeLists.txt
index 8180b98..ed2e4c8 100644
--- a/sources/CMakeLists.txt
+++ b/sources/CMakeLists.txt
@@ -3,9 +3,11 @@
 # BUILD THE STATIC LIBRARY
 ##########################
 
-add_library(Meshb.7 libmeshb7.c)
-install (FILES libmeshb7.h DESTINATION include COMPONENT headers)
-target_include_directories(Meshb.7 PUBLIC $<INSTALL_INTERFACE:include>)
+add_library(Meshb.7 libmeshb7.c ../utilities/libmeshb7_helpers.c)
+install (FILES libmeshb7.h ../utilities/libmeshb7_helpers.h DESTINATION include COMPONENT headers)
+target_include_directories(
+  Meshb.7 PUBLIC $<INSTALL_INTERFACE:include>
+                 $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/sources>)
 install (TARGETS Meshb.7 EXPORT meshb-target DESTINATION lib COMPONENT libraries)
 install (EXPORT meshb-target NAMESPACE ${PROJECT_NAME}::
          DESTINATION lib/cmake/${PROJECT_NAME})
@@ -18,7 +20,7 @@ export  (PACKAGE libMeshb)
 # BUILD THE LIBRARY WITH FORTRAN API
 ####################################
 
-if (CMAKE_Fortran_COMPILER)
+if (CMAKE_Fortran_COMPILER AND WITH_GMF_FORTRAN)
    file(GLOB_RECURSE SOURCES *.[chfF] *.[fF]90)
    add_library(Meshbf.7 STATIC ${SOURCES})
    set_target_properties(Meshbf.7 PROPERTIES COMPILE_FLAGS "-DF77API")
diff --git a/sources/libmeshb7.c b/sources/libmeshb7.c
index b7d64a8..d2fa30d 100644
--- a/sources/libmeshb7.c
+++ b/sources/libmeshb7.c
@@ -2,14 +2,14 @@
 
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
-/*                               LIBMESHB V7.60                               */
+/*                               LIBMESHB V7.62                               */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
 /*   Description:        handles .meshb file format I/O                       */
 /*   Author:             Loic MARECHAL                                        */
 /*   Creation date:      dec 09 1999                                          */
-/*   Last modification:  aug 04 2021                                          */
+/*   Last modification:  jan 07 2022                                          */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 
@@ -128,7 +128,7 @@
 
 // AIO: hardware or software mockup are both encapsulated into my_aio functions
 
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
 
 #include <aio.h>
 
@@ -461,6 +461,15 @@ const char *GmfKwdFmt[ GmfMaxKwd + 1 ][3] =
    {"InnerPolygonVertices",                     "i", "i"},
    {"PolyhedraHeaders",                         "i", "ii"},
    {"PolyhedraFaces",                           "i", "i"},
+   {"Domains",                                  "",  "ii"},
+   {"VerticesGID",                              "i", "iii"},
+   {"EdgesGID",                                 "i", "iii"},
+   {"TrianglesGID",                             "i", "iii"},
+   {"QuadrilateralsGID",                        "i", "iii"},
+   {"TetrahedraGID",                            "i", "iii"},
+   {"PyramidsGID",                              "i", "iii"},
+   {"PrismsGID",                                "i", "iii"},
+   {"HexahedraGID",                             "i", "iii"},
 };
 
 #ifdef TRANSMESH
@@ -599,7 +608,7 @@ int64_t GmfOpenMesh(const char *FilNam, int mod, ...)
       if(msh->typ & Bin)
       {
          // Create the name string and open the file
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
          // [Bruno] added binary flag (necessary under Windows)
          msh->FilDes = open(msh->FilNam, OPEN_READ_FLAGS, OPEN_READ_MODE);
 
@@ -646,7 +655,7 @@ int64_t GmfOpenMesh(const char *FilNam, int mod, ...)
 
          do
          {
-            res = fscanf(msh->hdl, "%s", str);
+            res = fscanf(msh->hdl, "%100s", str);
          }while( (res != EOF) && strcmp(str, "MeshVersionFormatted") );
 
          if(res == EOF)
@@ -659,7 +668,7 @@ int64_t GmfOpenMesh(const char *FilNam, int mod, ...)
 
          do
          {
-            res = fscanf(msh->hdl, "%s", str);
+            res = fscanf(msh->hdl, "%100s", str);
          }while( (res != EOF) && strcmp(str, "Dimension") );
 
          if(res == EOF)
@@ -728,7 +737,7 @@ int64_t GmfOpenMesh(const char *FilNam, int mod, ...)
           * with a call to open(), because Windows needs the
           * binary flag to be specified.
           */
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
          msh->FilDes = open(msh->FilNam, OPEN_WRITE_FLAGS, OPEN_WRITE_MODE);
 
          if(msh->FilDes <= 0)
@@ -794,7 +803,7 @@ int GmfCloseMesh(int64_t MshIdx)
 
    // Close the file and free the mesh structure
    if(msh->typ & Bin)
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
       close(msh->FilDes);
 #else
       fclose(msh->hdl);
@@ -1162,6 +1171,16 @@ int NAMF77(GmfSetLin, gmfsetlin)(TYPF77(int64_t) MshIdx, TYPF77(int) KwdCod, ...
       return(0);
    }
 
+   // Save the current stack environment for longjmp
+   // This is needed in RecBlk()
+   if( (err = setjmp(msh->err)) != 0)
+   {
+#ifdef GMFDEBUG
+      printf("libMeshb : mesh %p : error %d\n", msh, err);
+#endif
+      return(0);
+   }
+
    // Start decoding the arguments
    va_start(VarArg, KwdCod);
 
@@ -1386,7 +1405,7 @@ int GmfCpyLin(int64_t InpIdx, int64_t OutIdx, int KwdCod)
          if(InpMsh->typ & Asc)
             safe_fgets(s, WrdSiz * FilStrSiz, InpMsh->hdl, InpMsh->err);
          else
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
             read(InpMsh->FilDes, s, WrdSiz * FilStrSiz);
 #else
             safe_fread(s, WrdSiz, FilStrSiz, InpMsh->hdl, InpMsh->err);
@@ -1394,7 +1413,7 @@ int GmfCpyLin(int64_t InpIdx, int64_t OutIdx, int KwdCod)
          if(OutMsh->typ & Asc)
             fprintf(OutMsh->hdl, "%s ", s);
          else
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
             write(OutMsh->FilDes, s, WrdSiz * FilStrSiz);
 #else
             fwrite(s, WrdSiz, FilStrSiz, OutMsh->hdl);
@@ -1664,7 +1683,7 @@ int NAMF77(GmfGetBlock, gmfgetblock)(  TYPF77(int64_t) MshIdx,
       memset(&aio, 0, sizeof(struct aiocb));
       FilBuf = BckBuf;
       aio.aio_buf = BckBuf;
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
       aio.aio_fildes = msh->FilDes;
 #else
       aio.aio_fildes = msh->hdl;
@@ -1726,7 +1745,7 @@ int NAMF77(GmfGetBlock, gmfgetblock)(  TYPF77(int64_t) MshIdx,
             {
                printf("block      = %zd / %zd\n", b+1, NmbBlk+1);
                printf("size       = "INT64_T_FMT" lines\n", BlkNmbLin);
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
                printf("aio_fildes = %d\n",aio.aio_fildes);
 #else
                printf("aio_fildes = %p\n",aio.aio_fildes);
@@ -2114,7 +2133,7 @@ int NAMF77(GmfSetBlock, gmfsetblock)(  TYPF77(int64_t) MshIdx,
       // Setup the asynchronous parameters
       memset(&aio, 0, sizeof(struct aiocb));
       FilBuf = BckBuf;
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
       aio.aio_fildes = msh->FilDes;
 #else
       aio.aio_fildes = msh->hdl;
@@ -2134,7 +2153,7 @@ int NAMF77(GmfSetBlock, gmfsetblock)(  TYPF77(int64_t) MshIdx,
             
             if(my_aio_write(&aio) == -1)
             {
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
                printf("aio_fildes = %d\n",aio.aio_fildes);
 #else
                printf("aio_fildes = %p\n",aio.aio_fildes);
@@ -2369,6 +2388,27 @@ int GmfSetHONodesOrdering(int64_t MshIdx, int KwdCod, int *BasTab, int *OrdTab)
       }
    }
 
+   // Check the ordering consistency
+   for(i=0;i<NmbNod;i++)
+   {
+      flg = 0;
+
+      for(j=0;j<NmbNod;j++)
+         if(kwd->OrdTab[j] == i)
+         {
+            flg = 1;
+            break;
+         }
+
+      if(!flg)
+      {
+         for(j=0;j<NmbNod;j++)
+            kwd->OrdTab[j] = j;
+
+         return(0);
+      }
+   }
+
    return(1);
 }
 
@@ -2509,7 +2549,7 @@ static int ScaKwdTab(GmfMshSct *msh)
    if(msh->typ & Asc)
    {
       // Scan each string in the file until the end
-      while(fscanf(msh->hdl, "%s", str) != EOF)
+      while(fscanf(msh->hdl, "%100s", str) != EOF)
       {
          // Fast test in order to reject quickly the numeric values
          if(isalpha(str[0]))
@@ -2729,7 +2769,7 @@ static void ExpFmt(GmfMshSct *msh, int KwdCod)
 
 static void ScaWrd(GmfMshSct *msh, void *ptr)
 {
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
    if(read(msh->FilDes, ptr, WrdSiz) != WrdSiz)
 #else
    if(fread(ptr, WrdSiz, 1, msh->hdl) != 1)
@@ -2747,7 +2787,7 @@ static void ScaWrd(GmfMshSct *msh, void *ptr)
 
 static void ScaDblWrd(GmfMshSct *msh, void *ptr)
 {
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
    if(read(msh->FilDes, ptr, WrdSiz * 2) != WrdSiz * 2)
 #else
    if( fread(ptr, WrdSiz, 2, msh->hdl) != 2 )
@@ -2787,7 +2827,7 @@ static int64_t GetPos(GmfMshSct *msh)
 static void RecWrd(GmfMshSct *msh, const void *wrd)
 {
    // [Bruno] added error control
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
    if(write(msh->FilDes, wrd, WrdSiz) != WrdSiz)
 #else
    if(fwrite(wrd, WrdSiz, 1, msh->hdl) != 1)
@@ -2803,7 +2843,7 @@ static void RecWrd(GmfMshSct *msh, const void *wrd)
 static void RecDblWrd(GmfMshSct *msh, const void *wrd)
 {
    // [Bruno] added error control
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
    if(write(msh->FilDes, wrd, WrdSiz * 2) != WrdSiz*2)
 #else
    if(fwrite(wrd, WrdSiz, 2, msh->hdl) != 2)
@@ -2841,14 +2881,14 @@ static void RecBlk(GmfMshSct *msh, const void *blk, int siz)
        * the cache size is 10000 words, this is much much smaller than 4Gb
        * so there is probably no problem.
        */
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
       if(write(msh->FilDes, msh->blk, (int)msh->pos) != (ssize_t)msh->pos)
 #else      
       if(fwrite(msh->blk, 1, (size_t)msh->pos, msh->hdl) != msh->pos)
 #endif      
          longjmp(msh->err, -30);
 #else      
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
       if(write(msh->FilDes, msh->blk, msh->pos) != (ssize_t)msh->pos)
 #else      
       if(fwrite(msh->blk, 1, msh->pos, msh->hdl) != msh->pos)
@@ -2902,7 +2942,7 @@ static void SwpWrd(char *wrd, int siz)
 
 static int SetFilPos(GmfMshSct *msh, int64_t pos)
 {
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
    if(msh->typ & Bin)
       return((lseek(msh->FilDes, (off_t)pos, 0) != -1));
    else
@@ -2919,7 +2959,7 @@ static int SetFilPos(GmfMshSct *msh, int64_t pos)
 
 static int64_t GetFilPos(GmfMshSct *msh)
 {
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
    if(msh->typ & Bin)
       return(lseek(msh->FilDes, 0, 1));
    else
@@ -2940,7 +2980,7 @@ static int64_t GetFilSiz(GmfMshSct *msh)
 
    if(msh->typ & Bin)
    {
-#ifdef WITH_AIO
+#ifdef WITH_GMF_AIO
       CurPos = lseek(msh->FilDes, 0, 1);
       EndPos = lseek(msh->FilDes, 0, 2);
       lseek(msh->FilDes, (off_t)CurPos, 0);
diff --git a/sources/libmeshb7.h b/sources/libmeshb7.h
index 7e81ba0..bd67928 100644
--- a/sources/libmeshb7.h
+++ b/sources/libmeshb7.h
@@ -2,14 +2,14 @@
 
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
-/*                               LIBMESHB V7.60                               */
+/*                               LIBMESHB V7.61                               */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
 /*    Description:        handle .meshb file format I/O                       */
 /*    Author:             Loic MARECHAL                                       */
 /*    Creation date:      dec 09 1999                                         */
-/*    Last modification:  mar 24 2021                                         */
+/*    Last modification:  sep 27 2021                                         */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 
@@ -259,6 +259,15 @@ enum GmfKwdCod
    GmfInnerPolygonVertices,
    GmfPolyhedraHeaders,
    GmfPolyhedraFaces,
+   GmfDomains,
+   GmfVerticesGID,
+   GmfEdgesGID,
+   GmfTrianglesGID,
+   GmfQuadrilateralsGID,
+   GmfTetrahedraGID,
+   GmfPyramidsGID,
+   GmfPrismsGID,
+   GmfHexahedraGID,
    GmfLastKeyword
 };
 
diff --git a/todolist.md b/todolist.md
index c991da8..5a5305f 100644
--- a/todolist.md
+++ b/todolist.md
@@ -1,8 +1,7 @@
 ## HIGH PRIORITY
 
-### Handle arbitrary degree high-order elements
-- Setup a keyword to store Pk elements.
-- Give along a routine to convert to and from well-known high-order numberings.
+### Develop a helper to build the list of Ghost elements
+- Read and allocate the required ElementGID, then build the list of elements that do not belong to the local domain and return it to the caller.
 
 ### Distributed parallel write
 - Open a mesh file in writing mode but only create the skeleton of the mesh structure to enable further concurrent write access.
@@ -10,6 +9,10 @@
 
 ## STANDARD PRIORITY
 
+### Handle arbitrary degree high-order elements
+- Setup a keyword to store Pk elements.
+- Give along a routine to convert to and from well-known high-order numberings.
+
 ### Handle arbitrary degree polygons and polyhedra
 - Add a helper that cuts a polyhedron through a plane and generates the intersection's triangulated mesh in an STL-like format.
 
@@ -50,12 +53,13 @@ for example
 - Setup a polyhedron keyword that provides an arbitrary degree and number of oriented polygons.
 - Setup a boundary polygon keyword that lists the boundary faces among the ensemble of volume polygons.
 - Create a separate helper file to be optionally compiled along the libMeshb.
-- Add a helper to allocate and read the boundary polygons.
-- Add a helper to allocate and read the polyhedra and inner polygons.
-- Add a helper that returns a polygon's degree and nodes list.
-- Add a helper that returns a polyhedron's degree and face indices list.
+- Added a helper to allocate and read the boundary polygons.
+- Added a helper to allocate and read the polyhedra and inner polygons.
+- Added a helper that returns a polygon's degree and nodes list.
+- Added a helper that returns a polyhedron's degree and face indices list.
 - Write the helper's documentation.
 - Added an example that converts and fuses all surfaces and volume elements to polygons and polyhedra.
-- Add a helper to tesselate a polygon into a set of triangles made of nodes and edge visibility flags.
-- Add a set of keywords for each kind of element that stores the number of vertices and their barycentric coordinates.
-- Add a section about the AIO mode (Asynchronous Input Output).
+- Added a helper to tesselate a polygon into a set of triangles made of nodes and edge visibility flags.
+- Added a set of keywords for each kind of element that stores the number of vertices and their barycentric coordinates.
+- Added a section about the AIO mode (Asynchronous Input Output).
+- Added a helper that evaluates the quality of a mesh numbering in terms of efficient cache reuse and inner concurrency.
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index bf28305..632c0ea 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -7,3 +7,15 @@ add_executable(transmesh transmesh.c ../sources/libmeshb7.c)
 target_compile_options (transmesh PRIVATE -DTRANSMESH)
 target_link_libraries (transmesh Meshb.7 ${AIO_LIBRARIES})
 install (TARGETS transmesh DESTINATION bin COMPONENT applications)
+
+add_executable(mshinfo mshinfo.c)
+target_link_libraries (mshinfo Meshb.7 ${AIO_LIBRARIES})
+install (TARGETS mshinfo DESTINATION bin COMPONENT applications)
+
+add_executable(speed_test speed_test.c)
+target_link_libraries (speed_test Meshb.7 ${AIO_LIBRARIES})
+install (TARGETS speed_test DESTINATION bin COMPONENT applications)
+
+add_executable(ugrid2mesh ugrid2mesh.c)
+target_link_libraries (ugrid2mesh Meshb.7 ${AIO_LIBRARIES})
+install (TARGETS ugrid2mesh DESTINATION bin COMPONENT applications)
diff --git a/utilities/libmeshb7_helpers.c b/utilities/libmeshb7_helpers.c
index 7260445..e4a369b 100644
--- a/utilities/libmeshb7_helpers.c
+++ b/utilities/libmeshb7_helpers.c
@@ -2,14 +2,14 @@
 
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
-/*                            LIBMESHB-HELPERS V0.9                           */
+/*                            LIBMESHB-HELPERS V0.91                          */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
-/* Description:         set of helpers functions useful for the libMeshb      */
+/* Description:         Set of helpers functions useful for the libMeshb      */
 /* Author:              Loic MARECHAL                                         */
 /* Creation date:       mar 24 2021                                           */
-/* Last modification:   aug 05 2021                                           */
+/* Last modification:   nov 05 2021                                           */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 
@@ -274,7 +274,13 @@ itg GmfGetBoundaryPolygon(PolMshSct *pol, itg EleIdx, itg *UsrTab)
 
 itg GmfGetInnerPolygon(PolMshSct *pol, itg EleIdx, itg *UsrTab)
 {
-   itg i, BegIdx, EndIdx;
+   itg i, BegIdx, EndIdx, ord = 1;
+
+   if(EleIdx < 0)
+   {
+      EleIdx = -EleIdx;
+      ord = -1;
+   }
 
    if (!pol || (EleIdx <= 0) || (EleIdx > pol->NmbInrHdr))
       return (0);
@@ -288,8 +294,12 @@ itg GmfGetInnerPolygon(PolMshSct *pol, itg EleIdx, itg *UsrTab)
 
    EndIdx = MIN(EndIdx, BegIdx + 256);
 
-   for (i = BegIdx; i <= EndIdx; i++)
-      UsrTab[i - BegIdx] = pol->InrVerTab[i];
+   if(ord == 1)
+      for (i = BegIdx; i <= EndIdx; i++)
+         UsrTab[i - BegIdx] = pol->InrVerTab[i];
+   else
+      for (i = EndIdx; i >= BegIdx; i--)
+         UsrTab[EndIdx - i] = pol->InrVerTab[i];
 
    return (EndIdx - BegIdx + 1);
 }
@@ -367,3 +377,37 @@ itg GmfTesselatePolygon(PolMshSct *pol, itg EleIdx, itg (*Tri)[3], itg (*VisEdg)
 
    return NmbTri;
 }
+
+
+/*----------------------------------------------------------------------------*/
+/* Return the percentage of cache hit while doing element to node accesses    */
+/*----------------------------------------------------------------------------*/
+
+float GmfEvaluateNumbering(int NmbEle, int NmbNod, int *PtrEle, int *PtrEnd)
+{
+   int i, j, EleStp, HshTab[256] = {0}, NmbHit = 0;
+
+   if(!NmbEle || !NmbNod || !PtrEle || !PtrEnd || (PtrEnd < PtrEle))
+      return(0.);
+
+   // Compute the user data structure stride between two lines of elements
+   EleStp = (NmbEle > 1) ? (PtrEnd - PtrEle) / (NmbEle - 1) : 0;
+
+   // Simulate the search and insert in a 256 entries hash table
+   for(i=0;i<NmbEle;i++)
+   {
+      // Get the next element's pointer on its nodes
+      PtrEle += EleStp;
+
+      // Search for each node in the hash table:
+      // if its found, increase the hit counter and insert it otherwise
+      for(j=0;j<NmbNod;j++)
+         if(HshTab[ PtrEle[j] & 0xff ] == PtrEle[j])
+            NmbHit++;
+         else
+            HshTab[ PtrEle[j] & 0xff ] = PtrEle[j];
+   }
+
+   // Return the percentage of cache hits
+   return( (100. * NmbHit) / (NmbEle * NmbNod) );
+}
diff --git a/utilities/libmeshb7_helpers.h b/utilities/libmeshb7_helpers.h
index ca78410..d4c2d67 100644
--- a/utilities/libmeshb7_helpers.h
+++ b/utilities/libmeshb7_helpers.h
@@ -2,14 +2,14 @@
 
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
-/*                            LIBMESHB-HELPERS V0.9                           */
+/*                            LIBMESHB-HELPERS V0.91                          */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 /*                                                                            */
 /* Description:         libmeshb's helper functions' headers                  */
 /* Author:              Loic MARECHAL                                         */
 /* Creation date:       mar 24 2021                                           */
-/* Last modification:   aug 13 2021                                           */
+/* Last modification:   nov 05 2021                                           */
 /*                                                                            */
 /*----------------------------------------------------------------------------*/
 
@@ -64,6 +64,7 @@ extern itg        GmfGetBoundaryPolygon         (PolMshSct *, itg, itg *);
 extern itg        GmfGetInnerPolygon            (PolMshSct *, itg, itg *);
 extern itg        GmfGetPolyhedron              (PolMshSct *, itg, itg *);
 extern itg        GmfTesselatePolygon           (PolMshSct *, itg, itg (*)[3], itg (*)[3]);
+extern float      GmfEvaluateNumbering          (int, int, int *, int *);
 
 #ifdef __cplusplus
 } // end extern "C"
diff --git a/utilities/libmeshb7_helpers.md b/utilities/libmeshb7_helpers.md
index 70edc6c..ecaff06 100644
--- a/utilities/libmeshb7_helpers.md
+++ b/utilities/libmeshb7_helpers.md
@@ -87,3 +87,12 @@ Synopsis: extract a polyhedron's face list from the main global table and copy t
 
 - input: a pointer to an allocated PolMshSct structure, the index of a polyhedral element, a pointer to a table of **itg** that you have previously allocated with at least 256 entries.
 - output: the polyhedron's number of faces is the function's returned value (0 in case of failure) and the provided table is filled with the polyhedron's face indices.
+
+### GmlEvaluateNumbering
+Synopsis: Evaluate the quality of a mesh elements and nodes numbering. Usually, Delaunay meshes are quite bad (quality < 50) and need renumbering, octree or advancing front meshes are tolerable (from 50 to 90), and Hilbert renumbered meshes have a quality greater than 90.
+
+-input: the number of elements.
+-input: the number of nodes that make-up this kind of element (edge = 2, tetrahedron = 4).
+-input pointer to the first node of the first element.
+-input pointer to the first node of the last element.
+-output: a quality factor that ranges from 0 (bad) to 100 (perfect).
diff --git a/utilities/mshinfo.c b/utilities/mshinfo.c
new file mode 100644
index 0000000..0a593aa
--- /dev/null
+++ b/utilities/mshinfo.c
@@ -0,0 +1,338 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "libmeshb7.h"
+
+
+/* simulate exception */
+#include <setjmp.h>
+jmp_buf ex_buf__;
+
+extern const char *GmfKwdFmt[ GmfMaxKwd + 1 ][3];
+
+int main(int argc, char *argv[])
+{
+  int        NmbVer, dim, deg, FilVer, NmbCor, NmbRidg;
+  int        NmbTet, NmbHex, NmbPri, NmbPyr, NmbTri, NmbQua, NmbEdg;
+  int        NmbTetP2, NmbHexQ2, NmbPriP2, NmbPyrP2, NmbTriP2, NmbQuaQ2, NmbEdgP2, NmbNod;
+  int        NmbTetP3, NmbHexQ3, NmbPriP3, NmbPyrP3, NmbTriP3, NmbQuaQ3, NmbEdgP3;
+  int        NmbTetP4, NmbHexQ4, NmbPriP4, NmbPyrP4, NmbTriP4, NmbQuaQ4, NmbEdgP4;
+  int        OrdTet, OrdHex, OrdPri, OrdPyr, OrdTri, OrdQua, OrdEdg;
+  int        OrdTetP2, OrdHexQ2, OrdPriP2, OrdPyrP2, OrdTriP2, OrdQuaQ2, OrdEdgP2;
+  int        OrdTetP3, OrdHexQ3, OrdPriP3, OrdPyrP3, OrdTriP3, OrdQuaQ3, OrdEdgP3;
+  int        OrdTetP4, OrdHexQ4, OrdPriP4, OrdPyrP4, OrdTriP4, OrdQuaQ4, OrdEdgP4;
+  int64_t    InpMsh = 0, InpSol = 0;
+  int        i, j, ite = 0;
+  float      flt;
+  double     dbl, time = 0.0;
+  int        NbrLin, SolSiz, NbrTyp, TypTab[ GmfMaxTyp ];
+  
+
+  deg    = 1;
+  NmbNod = 0;
+
+
+  if (argc <= 1) {
+    printf(" USAGE : mshinfo name [-check]\n");
+    exit(1);
+  }
+
+  char name[1024];
+  char sol[1024];
+
+  strcpy(name,argv[1]);
+  strcpy(sol,argv[1]);
+
+
+  char *sub = NULL;
+  sub = strstr(name,".mesh"); // check if it has the extension .mesh[b]
+
+
+  char *subsol = NULL;
+  subsol = strstr(name,".sol"); // check if it has the extension .sol[b]
+
+// Mesh part
+
+//-- test the reading of the file as a mesh
+  if (subsol == NULL) { //--- no extension sol[b], check if it is a mesh name without extnsion
+    if(sub != NULL) sol[sub-name]='\0'; // change file name to be opened as sol even if mesh extension is given
+
+    if (sub != NULL)
+      InpMsh = GmfOpenMesh(name, GmfRead, &FilVer, &dim);
+    else {
+      strcat(name,".meshb");
+      InpMsh = GmfOpenMesh(name, GmfRead, &FilVer, &dim);
+      if(InpMsh == 0){
+        int ln = strlen(name);
+        name[ln -1] = '\0';
+        InpMsh = GmfOpenMesh(name, GmfRead, &FilVer, &dim);
+      }
+    }
+  }
+
+  if ( InpMsh != 0 && subsol == NULL ) {
+    printf("Mesh informations :\n");
+   /* Get number of entities*/
+       
+    NmbVer   = GmfStatKwd(InpMsh, GmfVertices);
+    NmbCor   = GmfStatKwd(InpMsh, GmfCorners);
+    NmbRidg  = GmfStatKwd(InpMsh, GmfRidges);
+    NmbTet   = GmfStatKwd(InpMsh, GmfTetrahedra);
+    NmbHex   = GmfStatKwd(InpMsh, GmfHexahedra);
+    NmbPyr   = GmfStatKwd(InpMsh, GmfPyramids);
+    NmbPri   = GmfStatKwd(InpMsh, GmfPrisms);
+    NmbTri   = GmfStatKwd(InpMsh, GmfTriangles);
+    NmbQua   = GmfStatKwd(InpMsh, GmfQuadrilaterals);
+    NmbEdg   = GmfStatKwd(InpMsh, GmfEdges);
+    NmbTetP2 = GmfStatKwd(InpMsh, GmfTetrahedraP2);
+    NmbHexQ2 = GmfStatKwd(InpMsh, GmfHexahedraQ2);
+    NmbPyrP2 = GmfStatKwd(InpMsh, GmfPyramidsP2);
+    NmbPriP2 = GmfStatKwd(InpMsh, GmfPrismsP2);
+    NmbTriP2 = GmfStatKwd(InpMsh, GmfTrianglesP2);
+    NmbQuaQ2 = GmfStatKwd(InpMsh, GmfQuadrilateralsQ2);
+    NmbEdgP2 = GmfStatKwd(InpMsh, GmfEdgesP2);
+    NmbTetP3 = GmfStatKwd(InpMsh, GmfTetrahedraP3);
+    NmbHexQ3 = GmfStatKwd(InpMsh, GmfHexahedraQ3);
+    NmbPyrP3 = GmfStatKwd(InpMsh, GmfPyramidsP3);
+    NmbPriP3 = GmfStatKwd(InpMsh, GmfPrismsP3);
+    NmbTriP3 = GmfStatKwd(InpMsh, GmfTrianglesP3);
+    NmbQuaQ3 = GmfStatKwd(InpMsh, GmfQuadrilateralsQ3);
+    NmbEdgP3 = GmfStatKwd(InpMsh, GmfEdgesP3);
+    NmbTetP4 = GmfStatKwd(InpMsh, GmfTetrahedraP4);
+    NmbHexQ4 = GmfStatKwd(InpMsh, GmfHexahedraQ4);
+    NmbPyrP4 = GmfStatKwd(InpMsh, GmfPyramidsP4);
+    NmbPriP4 = GmfStatKwd(InpMsh, GmfPrismsP4);
+    NmbTriP4 = GmfStatKwd(InpMsh, GmfTrianglesP4);
+    NmbQuaQ4 = GmfStatKwd(InpMsh, GmfQuadrilateralsQ4);
+    NmbEdgP4 = GmfStatKwd(InpMsh, GmfEdgesP4);
+
+    //-- Nodes ordering
+    OrdTet   = GmfStatKwd(InpMsh, GmfTetrahedraP1Ordering);
+    OrdHex   = GmfStatKwd(InpMsh, GmfHexahedraQ1Ordering);
+    OrdPyr   = GmfStatKwd(InpMsh, GmfPyramidsP1Ordering);
+    OrdPri   = GmfStatKwd(InpMsh, GmfPrismsP1Ordering);
+    OrdTri   = GmfStatKwd(InpMsh, GmfTrianglesP1Ordering);
+    OrdQua   = GmfStatKwd(InpMsh, GmfQuadrilateralsQ1Ordering);
+    OrdEdg   = GmfStatKwd(InpMsh, GmfEdgesP1Ordering);
+    OrdTetP2 = GmfStatKwd(InpMsh, GmfTetrahedraP2Ordering);
+    OrdHexQ2 = GmfStatKwd(InpMsh, GmfHexahedraQ2Ordering);
+    OrdPyrP2 = GmfStatKwd(InpMsh, GmfPyramidsP2Ordering);
+    OrdPriP2 = GmfStatKwd(InpMsh, GmfPrismsP2Ordering);
+    OrdTriP2 = GmfStatKwd(InpMsh, GmfTrianglesP2Ordering);
+    OrdQuaQ2 = GmfStatKwd(InpMsh, GmfQuadrilateralsQ2Ordering);
+    OrdEdgP2 = GmfStatKwd(InpMsh, GmfEdgesP2Ordering);
+    OrdTetP3 = GmfStatKwd(InpMsh, GmfTetrahedraP3Ordering);
+    OrdHexQ3 = GmfStatKwd(InpMsh, GmfHexahedraQ3Ordering);
+    OrdPyrP3 = GmfStatKwd(InpMsh, GmfPyramidsP3Ordering);
+    OrdPriP3 = GmfStatKwd(InpMsh, GmfPrismsP3Ordering);
+    OrdTriP3 = GmfStatKwd(InpMsh, GmfTrianglesP3Ordering);
+    OrdQuaQ3 = GmfStatKwd(InpMsh, GmfQuadrilateralsQ3Ordering);
+    OrdEdgP3 = GmfStatKwd(InpMsh, GmfEdgesP3Ordering);
+    OrdTetP4 = GmfStatKwd(InpMsh, GmfTetrahedraP4Ordering);
+    OrdHexQ4 = GmfStatKwd(InpMsh, GmfHexahedraQ4Ordering);
+    OrdPyrP4 = GmfStatKwd(InpMsh, GmfPyramidsP4Ordering);
+    OrdPriP4 = GmfStatKwd(InpMsh, GmfPrismsP4Ordering);
+    OrdTriP4 = GmfStatKwd(InpMsh, GmfTrianglesP4Ordering);
+    OrdQuaQ4 = GmfStatKwd(InpMsh, GmfQuadrilateralsQ4Ordering);
+    OrdEdgP4 = GmfStatKwd(InpMsh, GmfEdgesP4Ordering);
+ 
+    printf("version = %d; dim = %d; nbv = %d; (nbc = %d)\n", FilVer, dim, NmbVer, NmbCor);
+    
+    if ( NmbTet > 0 || NmbTri > 0 || NmbEdg > 0 || NmbHex > 0|| NmbPri > 0 || NmbPyr > 0|| NmbQua > 0 ) {
+      printf("P1/Q1 elements :\n");
+      if ( NmbTet >0 || NmbTri > 0 || NmbEdg > 0 ) 
+        printf("nbtet = %d; nbtri = %d; nbedg = %d; (nbrdg = %d)\n", NmbTet, NmbTri, NmbEdg, NmbRidg);
+      if ( NmbHex > 0|| NmbPri > 0 || NmbPyr > 0|| NmbQua > 0 )
+        printf("nbhex = %d; nbpri = %d; nbpyr = %d; nbqua = %d\n", NmbHex, NmbPri, NmbPyr, NmbQua);
+    }
+    
+    if ( OrdTet >0 || OrdTri > 0 || OrdEdg > 0 || OrdHex > 0|| OrdPri > 0 || OrdPyr > 0|| OrdQua > 0 ) {
+      printf("P1/Q1 ordering :\n");
+      if ( OrdTet > 0 || OrdTri > 0 || OrdEdg > 0 ) {
+        if ( OrdTet > 0 )
+          printf("tet ");
+        if ( OrdTri > 0 )
+          printf("tri ");
+        if ( OrdEdg > 0 )
+          printf("edg ");
+        printf("\n");
+      }
+      if ( OrdHex > 0|| OrdPri > 0 || OrdPyr > 0|| OrdQua > 0 ) {
+        if ( OrdHex > 0 )
+          printf("hex ");
+        if ( OrdPri > 0 )
+          printf("pri ");
+        if ( OrdPyr > 0 )
+          printf("pyr ");
+        if ( OrdQua > 0 )
+          printf("qua ");
+        printf("\n");
+      }
+    }
+
+    if ( NmbTetP2 > 0 || NmbTriP2 > 0 || NmbEdgP2 > 0 || NmbHexQ2 > 0|| NmbPriP2 > 0 || NmbPyrP2 > 0|| NmbQuaQ2 > 0 ) {
+      printf("P2/Q2 elements :\n"); 
+      if ( NmbTetP2 > 0 || NmbTriP2 > 0 || NmbEdgP2 > 0 )
+        printf("nbtetP2 = %d; nbtriP2 = %d; nbedgP2 = %d; (nbrdg = %d)\n", NmbTetP2, NmbTriP2, NmbEdgP2, NmbRidg);
+      if ( NmbHexQ2 > 0|| NmbPriP2 > 0 || NmbPyrP2 > 0|| NmbQuaQ2 > 0 )
+        printf("nbhexQ2 = %d; nbpriP2 = %d; nbpyrP2 = %d; nbquaQ2 = %d\n", NmbHexQ2, NmbPriP2, NmbPyrP2, NmbQuaQ2);
+    }
+    
+    if ( OrdTetP2 >0 || OrdTriP2 > 0 || OrdEdgP2 > 0 || OrdHexQ2 > 0|| OrdPriP2 > 0 || OrdPyrP2 > 0|| OrdQuaQ2 > 0 ) {
+      printf("P2/Q2 ordering :\n");
+      if ( OrdTetP2 > 0 || OrdTriP2 > 0 || OrdEdgP2 > 0 ) {
+        if ( OrdTetP2 > 0 )
+          printf("tetP2 ");
+        if ( OrdTriP2 > 0 )
+          printf("triP2 ");
+        if ( OrdEdgP2 > 0 )
+          printf("edgP2 ");
+        printf("\n");
+      }
+      if ( OrdHexQ2 > 0|| OrdPriP2 > 0 || OrdPyrP2 > 0|| OrdQuaQ2 > 0 ) {
+        if ( OrdHexQ2 > 0 )
+          printf("hexQ2 ");
+        if ( OrdPriP2 > 0 )
+          printf("priP2 ");
+        if ( OrdPyrP2 > 0 )
+          printf("pyrP2 ");
+        if ( OrdQuaQ2 > 0 )
+          printf("quaQ2 ");
+        printf("\n");
+      }
+    }
+
+    if ( NmbTetP3 > 0 || NmbTriP3 > 0 || NmbEdgP3 > 0 || NmbHexQ3 > 0|| NmbPriP3 > 0 || NmbPyrP3 > 0|| NmbQuaQ3 > 0 ) {
+      printf("P3/Q3 elements :\n"); 
+      if ( NmbTetP3 > 0 || NmbTriP3 > 0 || NmbEdgP3 > 0 )
+        printf("nbtetP3 = %d; nbtriP3 = %d; nbedgP3 = %d; (nbrdg = %d)\n", NmbTetP3, NmbTriP3, NmbEdgP3, NmbRidg);
+      if ( NmbHexQ3 > 0|| NmbPriP3 > 0 || NmbPyrP3 > 0|| NmbQuaQ3 > 0 )
+        printf("nbhexQ3 = %d; nbpriP3 = %d; nbpyrP3 = %d; nbquaQ3 = %d\n", NmbHexQ3, NmbPriP3, NmbPyrP3, NmbQuaQ3);
+    }
+
+    if ( OrdTetP3 >0 || OrdTriP3 > 0 || OrdEdgP3 > 0 || OrdHexQ3 > 0|| OrdPriP3 > 0 || OrdPyrP3 > 0|| OrdQuaQ3 > 0 ) {
+      printf("P3/Q3 ordering :\n");
+      if ( OrdTetP3 > 0 || OrdTriP3 > 0 || OrdEdgP3 > 0 ) {
+        if ( OrdTetP3 > 0 )
+          printf("tetP3 ");
+        if ( OrdTriP3 > 0 )
+          printf("triP3 ");
+        if ( OrdEdgP3 > 0 )
+          printf("edgP3 ");
+        printf("\n");
+      }
+      if ( OrdHexQ3 > 0|| OrdPriP3 > 0 || OrdPyrP3 > 0|| OrdQuaQ3 > 0 ) {
+        if ( OrdHexQ3 > 0 )
+          printf("hexQ3 ");
+        if ( OrdPriP3 > 0 )
+          printf("priP3 ");
+        if ( OrdPyrP3 > 0 )
+          printf("pyrP3 ");
+        if ( OrdQuaQ3 > 0 )
+          printf("quaQ3 ");
+        printf("\n");
+      }
+    }
+
+    if ( NmbTetP4 > 0 || NmbTriP4 > 0 || NmbEdgP4 > 0 || NmbHexQ4 > 0|| NmbPriP4 > 0 || NmbPyrP4 > 0|| NmbQuaQ4 > 0 ) {
+      printf("P4/Q4 elements :\n"); 
+      if ( NmbTetP4 > 0 || NmbTriP4 > 0 || NmbEdgP4 > 0 )
+        printf("nbtetP4 = %d; nbtriP4 = %d; nbedgP4 = %d; (nbrdg = %d)\n", NmbTetP4, NmbTriP4, NmbEdgP4, NmbRidg);
+      if ( NmbHexQ4 > 0|| NmbPriP4 > 0 || NmbPyrP4 > 0|| NmbQuaQ4 > 0 )
+        printf("nbhexQ4 = %d; nbpriP4 = %d; nbpyrP4 = %d; nbquaQ4 = %d\n", NmbHexQ4, NmbPriP4, NmbPyrP4, NmbQuaQ4);
+    }
+
+    if ( OrdTetP4 >0 || OrdTriP4 > 0 || OrdEdgP4 > 0 || OrdHexQ4 > 0|| OrdPriP4 > 0 || OrdPyrP4 > 0|| OrdQuaQ4 > 0 ) {
+      printf("P4/Q4 ordering :\n");
+      if ( OrdTetP4 > 0 || OrdTriP4 > 0 || OrdEdgP4 > 0 ) {
+        if ( OrdTetP4 > 0 )
+          printf("tetP4 ");
+        if ( OrdTriP4 > 0 )
+          printf("triP4 ");
+        if ( OrdEdgP4 > 0 )
+          printf("edgP4 ");
+        printf("\n");
+      }
+      if ( OrdHexQ4 > 0|| OrdPriP4 > 0 || OrdPyrP4 > 0|| OrdQuaQ4 > 0 ) {
+        if ( OrdHexQ4 > 0 )
+          printf("hexQ4 ");
+        if ( OrdPriP4 > 0 )
+          printf("priP4 ");
+        if ( OrdPyrP4 > 0 )
+          printf("pyrP4 ");
+        if ( OrdQuaQ4 > 0 )
+          printf("quaQ4 ");
+        printf("\n");
+      }
+    }
+
+    GmfCloseMesh(InpMsh);
+     
+  }
+  
+
+  // solution part
+
+//-- test the reading of the file as a solution
+  if (subsol != NULL) { //--- a file containing char ".sol" is given in input
+    InpSol = GmfOpenMesh(sol, GmfRead, &FilVer, &dim);
+    if ( InpSol == 0 ) {
+      strcat(sol,"b");
+      InpSol = GmfOpenMesh(sol, GmfRead, &FilVer, &dim);
+    }
+  }
+  else { //-- check if a solution with the same name as the mesh one exists
+    strcat(sol,".sol");
+    InpSol = GmfOpenMesh(sol, GmfRead, &FilVer, &dim);
+    if ( InpSol == 0 ) {
+      strcat(sol,"b");
+      InpSol = GmfOpenMesh(sol, GmfRead, &FilVer, &dim);
+    }
+  }
+
+  if (InpSol!= 0 ) {
+    printf("Solution informations :\n");
+    if ( GmfStatKwd(InpSol, GmfTime) ) {
+      GmfGotoKwd(InpSol, GmfTime);
+      if ( FilVer == GmfFloat ) {	// read 32 bits float
+       GmfGetLin(InpSol, GmfTime, &flt);
+       time = (double)flt;
+      }  
+      else if ( FilVer == GmfDouble ) {	// read 64 bits float
+        GmfGetLin(InpSol, GmfTime, &dbl);
+        time = dbl;
+      }
+    }
+    if ( GmfStatKwd(InpSol, GmfIterations) ) {
+      GmfGotoKwd(InpSol, GmfIterations);
+      GmfGetLin( InpSol, GmfIterations, &ite);  
+    }  
+    printf("dim = %d; ite = %d; time = %lg\n",dim, ite, time);
+    for(i=1; i<=GmfMaxKwd; i++) {
+      if( ( (!strcmp(GmfKwdFmt[i][2], "sr")  || !strcmp(GmfKwdFmt[i][2], "hr")) ) && ( (NbrLin = GmfStatKwd(InpSol, i, &NbrTyp, &SolSiz, TypTab, &deg, &NmbNod)) ) ) {
+        printf("%s = %d\n", GmfKwdFmt[i][0], NbrLin);
+        if ( deg != 1 && NmbNod != 0 )
+          printf("deg = %d; nbnod = %d\n",deg, NmbNod);
+        printf("type = [");
+        for(j=0; j<NbrTyp; j++) {
+          if (      TypTab[j] == 1 ) printf(" scalar ");
+          else if ( TypTab[j] == 2 ) printf(" vector ");
+          else if ( TypTab[j] == 3 ) printf(" metric ");
+          else if ( TypTab[j] == 4 ) printf(" matrix ");
+        }
+        printf("]\n");
+        deg = 1;
+        NmbNod = 0;
+      }
+    }
+    GmfCloseMesh(InpSol);
+  }
+
+  if ( InpSol == 0 && InpMsh == 0 ) {
+    printf("cannot open mesh/solution file %s\n", argv[1]);
+    exit(1);
+  }
+
+  return(0);
+
+}
\ No newline at end of file
diff --git a/utilities/speed_test.c b/utilities/speed_test.c
new file mode 100644
index 0000000..cc81855
--- /dev/null
+++ b/utilities/speed_test.c
@@ -0,0 +1,110 @@
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#ifdef _WIN32
+#include <sys/timeb.h>
+#else
+#include <sys/time.h>
+#endif
+
+#include <libmeshb7.h>
+
+
+/*----------------------------------------------------------------------------*/
+/* Starts or stops the given timer                                            */
+/*----------------------------------------------------------------------------*/
+
+double GetWallClock()
+{
+#ifdef WIN32
+   struct __timeb64 tb;
+   _ftime64(&tb);
+   return((double)tb.time + (double)tb.millitm/1000.);
+#else
+   struct timeval tp;
+   gettimeofday(&tp, NULL);
+   return(tp.tv_sec + tp.tv_usec / 1000000.);
+#endif
+}
+
+
+/*----------------------------------------------------------------------------*/
+/*----------------------------------------------------------------------------*/
+
+int main()
+{
+   int   i, NmbVer, NmbTet, ver, dim, *RefTab, (*TetTab)[5];
+   int64_t InpMsh, OutMsh;
+   double (*VerTab)[3];
+   double timer;
+
+
+   /*-----------------------------------*/
+   /* Open mesh file "tets.meshb"       */
+   /*-----------------------------------*/
+
+   timer = GetWallClock();
+
+   if(!(InpMsh = GmfOpenMesh("tets.meshb", GmfRead, &ver, &dim)))
+      return(1);
+
+   printf("InpMsh : idx = %lld, version = %d, dimension = %d\n", InpMsh, ver, dim);
+
+   // Read the number of vertices and allocate memory
+   NmbVer = GmfStatKwd(InpMsh, GmfVertices);
+   printf("InpMsh : nmb vertices = %d\n", NmbVer);
+   VerTab = malloc((size_t)(NmbVer+1) * 3 * sizeof(double));
+   RefTab = malloc((size_t)(NmbVer+1)     * sizeof(int));
+
+   // Read the number of tets and allocate memory
+   NmbTet = GmfStatKwd(InpMsh, GmfTetrahedra);
+   printf("InpMsh : nmb tets = %d\n", NmbTet);
+   TetTab = malloc((size_t)(NmbTet+1) * 5 * sizeof(int));
+
+   // Read the vertices
+   GmfGetBlock(InpMsh, GmfVertices, 1, NmbVer, 0, NULL, NULL,
+               GmfDoubleVec, 3, &VerTab[1][0], &VerTab[ NmbVer ][0],
+               GmfInt,         &RefTab[1],    &RefTab[ NmbVer ] );
+
+   // Read the tets
+   GmfGetBlock(InpMsh, GmfTetrahedra, 1, NmbTet, 0, NULL, NULL,
+               GmfIntVec, 5, &TetTab[1][0], &TetTab[ NmbTet ][0]);
+
+   // Close the tet mesh
+   GmfCloseMesh(InpMsh);
+
+   printf("Time for reading: %g seconds\n", GetWallClock() - timer);
+
+
+   /*-----------------------------------*/
+   /* Write the mesh back               */
+   /*-----------------------------------*/
+
+   timer = GetWallClock();
+
+   if(!(OutMsh = GmfOpenMesh("tets_out.meshb", GmfWrite, ver, dim)))
+      return(1);
+
+   // Write the vertices
+   GmfSetKwd(OutMsh, GmfVertices, NmbVer);
+   GmfSetBlock(OutMsh, GmfVertices, 1, NmbVer, 0, NULL, NULL,
+               GmfDoubleVec, 3, &VerTab[1][0], &VerTab[ NmbVer ][0],
+               GmfInt,         &RefTab[1],    &RefTab[ NmbVer ] );
+
+   // Write the tets
+   GmfSetKwd(OutMsh, GmfTetrahedra, NmbTet);
+   GmfSetBlock(InpMsh, GmfTetrahedra, 1, NmbTet, 0, NULL, NULL,
+               GmfIntVec, 5, &TetTab[1][0], &TetTab[ NmbTet ][0]);
+
+   // Do not forget to close the mesh file
+   GmfCloseMesh(OutMsh);
+
+   printf("Time for writing: %g seconds\n", GetWallClock() - timer);
+
+   free(TetTab);
+   free(RefTab);
+   free(VerTab);
+
+   return(0);
+}
diff --git a/utilities/ugrid2mesh.c b/utilities/ugrid2mesh.c
new file mode 100644
index 0000000..c48e5e2
--- /dev/null
+++ b/utilities/ugrid2mesh.c
@@ -0,0 +1,380 @@
+/*
+
+Adrien Loseille, INRIA, 2016
+
+*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <libmeshb7.h>
+
+
+#define max(a,b) ((a)>(b))?(a):(b)
+
+typedef           int  mint;
+typedef long long int lint;
+
+typedef double double3[3];
+typedef  mint     int3[3];
+typedef  mint     int4[4];
+typedef  mint     int5[5];
+typedef  mint     int6[6];
+typedef  mint     int8[8];
+
+int littleEndian()
+{
+  int n = 1;
+  // little endian if true
+  if(*(char *)&n == 1) return 1;
+  else                 return 0;
+}
+
+
+double swapd(double d)
+{
+   double a;
+   unsigned char *dst = (unsigned char *)&a;
+   unsigned char *src = (unsigned char *)&d;
+
+   if ( ! littleEndian() ) return d;
+
+   dst[0] = src[7];
+   dst[1] = src[6];
+   dst[2] = src[5];
+   dst[3] = src[4];
+   dst[4] = src[3];
+   dst[5] = src[2];
+   dst[6] = src[1];
+   dst[7] = src[0];
+
+   return a;
+}
+
+int swapi4(int d)
+{
+   mint a;
+   unsigned char *dst = (unsigned char *)&a;
+   unsigned char *src = (unsigned char *)&d;
+
+   if ( ! littleEndian() ) return d;
+
+   dst[0] = src[3];
+   dst[1] = src[2];
+   dst[2] = src[1];
+   dst[3] = src[0];
+  
+   return a;
+}
+
+lint swapi8(lint d)
+{
+   lint a;
+   unsigned char *dst = (unsigned char *)&a;
+   unsigned char *src = (unsigned char *)&d;
+   
+   if ( ! littleEndian() ) return d;
+
+   dst[0] = src[7];
+   dst[1] = src[6];
+   dst[2] = src[5];
+   dst[3] = src[4];
+   dst[4] = src[3];
+   dst[5] = src[2];
+   dst[6] = src[1];
+   dst[7] = src[0];
+  
+   return a;
+}
+
+
+
+
+
+
+
+mint main(int argc, char *argv[])
+{
+  
+  mint nbrVer=0, nbrTri=0, nbrQua=0, nbrTet = 0, nbrPyr = 0, nbrPri = 0, nbrHex = 0;
+  int3 *tri        = NULL;
+  int4 *qua        = NULL;
+  int4 *tet        = NULL;
+  int5 *pyr        = NULL;
+  int6 *pri        = NULL;
+  int8 *hex        = NULL;
+  
+  mint  *lsurf     = NULL;
+  double3 *crd     = NULL;
+  
+  int    writeSurf = 0;
+  
+  mint iVer,iTri,iQua,iTet,iPyr,iPri,iHex,i,ix;
+  
+  if ( argc < 3 ) {
+    printf("  Usage : ugrid2mesh file.ugrid file.mesh[b] [ -surf -info ]\n");
+    return 0;
+  }
+  
+  /* write surface only */
+  if ( argc >= 4 ) {
+    if( strcmp(argv[3],"-surf") == 0 ) writeSurf = 1;
+  }
+  
+  
+  printf("  *************************************** \n");
+  printf("\n");
+  printf("  UGRID2MESH Conversion, version 1.0 \n");
+  printf("\n");
+  printf("  *************************************** \n");
+  
+  FILE *f = fopen(argv[1],"r");
+  if ( f == NULL ) {
+    printf("  ## ERROR: CANNOT OPEN FILE %s \n", argv[1]);
+    return 0;
+  }
+  
+  //--- should be used to see if it is necessary to swaps the bytes....
+  printf("  ** Local system is in  %s \n",(littleEndian() == 1)?"Little Endian encoding":"Big Endian encoding \n");
+  
+  //printf("  f= %p \n",f);
+  unsigned int size[7];
+  fread(size,sizeof(int),7,f); 
+  
+  //printf(" size %d %d %d %d %d %d %d \n",swapi4(size[0]),swapi4(size[1]),swapi4(size[2]),swapi4(size[3]),swapi4(size[4]),swapi4(size[5]),swapi4(size[6]));
+  
+  printf("  %% %s OPENED \n",argv[1]);
+  printf("  Number of Nodes           %15d\n",swapi4(size[0]));
+  printf("  Number of Triangles       %15d\n",swapi4(size[1]));
+  printf("  Number of Quads           %15d\n",swapi4(size[2]));
+  printf("  Number of Tetrahedra      %15d\n",swapi4(size[3]));
+  printf("  Number of Pyramids        %15d\n",swapi4(size[4]));
+  printf("  Number of Prisms          %15d\n",swapi4(size[5]));
+  printf("  Number of Hexes           %15d\n",swapi4(size[6]));
+  
+  
+  //printf(" size of the file %lg Gb\n",(swapi4(size[0])*3.0*8.0 + swapi4(size[1])*4.*4. + swapi4(size[3])*4.*4.)/(1024*1024*1024));
+  
+  /* exist here if just a query on the mesh */
+  printf(" argc %d \n",argc);
+  
+  if ( argc > 2 ) {
+    if ( strcmp(argv[2],"-info") == 0 ) return 1;
+    if ( argc > 3  ) {
+      if ( strcmp(argv[3],"-info") == 0 ) return 1; 
+    }
+  }
+  
+  printf(" -- allocating vertices \n");
+  nbrVer   = swapi4(size[0]);
+  crd = (double3 *)malloc(sizeof(double)*3*nbrVer);
+  
+  fread(crd,sizeof(double),3*nbrVer,f);
+  
+  printf(" -- reading vertices %d\n",nbrVer);
+  for(iVer=0; iVer<nbrVer; iVer++) {
+    for(i=0; i<3; i++) {
+      double sx  = swapd(crd[iVer][i]);
+      crd[iVer][i] = sx;
+    }
+    //if ( iVer <= 10 )        printf("iVer[%15d] =  %lg %lg %lg \n", iVer+1, crd[iVer][0],crd[iVer][1],crd[iVer][2]);
+    //if ( iVer >= nbrVer-10 ) printf("iVer[%15d] =  %lg %lg %lg \n", iVer+1,crd[iVer][0],crd[iVer][1],crd[iVer][2]);
+    //printf("iVer[%d] =  %lg %lg %lg \n", iVer+1, crd[iVer][0],crd[iVer][1],crd[iVer][2]);
+  }
+    
+  //--- write only surface mesh ?
+  printf(" -- allocating surface link \n");
+  mint  mark  = 1;
+  int  *pmark = malloc(sizeof(mint)*(nbrVer+1));
+  
+  mint npsrf  = 0;
+  mint *lpsrf = malloc(sizeof(mint)*(nbrVer+1));
+
+  memset(pmark,0,sizeof(mint)*(nbrVer+1));
+    
+  mint maxid = 0;
+  
+  printf(" -- reading triangles \n");
+  nbrTri  = swapi4(size[1]);
+  if ( nbrTri > 0 ) {
+    tri   = (int3 *)malloc(sizeof(mint)*3*nbrTri);
+    fread(tri,sizeof(mint),3*nbrTri,f);
+    for(iTri=0; iTri<nbrTri; iTri++) {
+      //printf(" iTri %d \n",iTri);
+      for(i=0; i<3; i++) {
+        ix  = swapi4(tri[iTri][i]);
+        tri[iTri][i] = ix;
+        maxid = max(maxid,ix);
+        if ( pmark[ix] == 0 ) {
+          pmark[ix] = mark;
+          lpsrf[npsrf] = ix;
+          npsrf = npsrf + 1;
+          pmark[ix]    = npsrf;
+        }
+        
+      }
+      //prmintf(" %lg %lg %lg \n",swapd(crd[0]),swapd(crd[1]),swapd(crd[2]));
+    }
+  }
+  
+  printf(" -- reading quads \n");
+  nbrQua  = swapi4(size[2]);
+  if ( nbrQua > 0 ) {
+    qua   = (int4 *)malloc(sizeof(mint)*4*nbrQua);
+    fread(qua,sizeof(mint),4*nbrQua,f);
+    for(iQua=0; iQua<nbrQua; iQua++) {
+      for(i=0; i<4; i++) {
+        ix  = swapi4(qua[iQua][i]);
+        qua[iQua][i] = ix;
+        maxid = max(maxid,ix);
+        if ( pmark[ix] == 0 ) {
+          pmark[ix] = mark;
+          lpsrf[npsrf] = ix;
+          npsrf = npsrf + 1;
+          pmark[ix]    = npsrf;
+        }
+      }
+      //prmintf(" %lg %lg %lg \n",swapd(crd[0]),swapd(crd[1]),swapd(crd[2]));
+    }
+  }
+  
+  printf(" -- setting up surface id \n");
+  lsurf = malloc(sizeof(mint)*(nbrTri+nbrQua));
+  fread(lsurf,sizeof(mint),nbrTri+nbrQua,f);
+  for(iTri=0; iTri<(nbrTri+nbrQua); iTri++) {
+     mint ix      = swapi4(lsurf[iTri]);
+     lsurf[iTri] = ix;
+  }
+  
+  printf(" ... end reading ugrid file \n");
+  
+  printf("  Number of Boundary  Nodes %15d \n",npsrf);
+  
+  int FilVer = 3;
+  int64_t OutMsh = GmfOpenMesh(argv[2], GmfWrite, FilVer, 3);
+  if ( OutMsh <= 0 ) {
+    printf("  ## ERROR: CANNOT CREATE FILE %s \n", argv[2]);
+    return 0;
+  }
+  
+  printf(" %% %s CREATED \n",argv[2]);
+  
+  if ( writeSurf == 1 ) {
+    printf("  outputting surface only \n");
+    GmfSetKwd(OutMsh, GmfVertices,npsrf);
+    for (iVer=0; iVer<npsrf; iVer++) {
+      GmfSetLin(OutMsh, GmfVertices, crd[lpsrf[iVer]-1][0],crd[lpsrf[iVer]-1][1],crd[lpsrf[iVer]-1][2],0); 
+    }
+    GmfSetKwd(OutMsh, GmfTriangles,nbrTri);
+    for (iTri=0; iTri<nbrTri; iTri++) {
+      GmfSetLin(OutMsh, GmfTriangles, pmark[tri[iTri][0]], pmark[tri[iTri][1]], pmark[tri[iTri][2]], lsurf[iTri]); 
+    }
+    GmfSetKwd(OutMsh, GmfQuadrilaterals,nbrQua);
+    for (iQua=0; iQua<nbrQua; iQua++) {
+      GmfSetLin(OutMsh, GmfQuadrilaterals, pmark[qua[iQua][0]], pmark[qua[iQua][1]], pmark[qua[iQua][2]],  pmark[qua[iQua][2]], lsurf[nbrTri+iQua]); 
+    }
+  }
+  else {
+    printf("  outputting volume only \n");
+    GmfSetKwd(OutMsh, GmfVertices, nbrVer);
+    for (iVer=0; iVer<nbrVer; iVer++) {
+      GmfSetLin(OutMsh, GmfVertices, crd[iVer][0],crd[iVer][1],crd[iVer][2],0); 
+    }
+    GmfSetKwd(OutMsh, GmfTriangles,nbrTri);
+    for (iTri=0; iTri<nbrTri; iTri++) {
+      GmfSetLin(OutMsh, GmfTriangles, tri[iTri][0], tri[iTri][1], tri[iTri][2], lsurf[iTri]); 
+    }
+    GmfSetKwd(OutMsh, GmfQuadrilaterals,nbrQua);
+    for (iQua=0; iQua<nbrQua; iQua++) {
+      GmfSetLin(OutMsh, GmfQuadrilaterals, qua[iQua][0], qua[iQua][1], qua[iQua][2],  qua[iQua][2], lsurf[nbrTri+iQua]); 
+    }
+    
+    nbrTet  = swapi4(size[3]);
+    if ( nbrTet > 0 ) {
+      tet   = (int4 *)malloc(sizeof(mint)*4*nbrTet);
+      fread(tet,sizeof(mint),4*nbrTet,f);
+      for(iTet=0; iTet<nbrTet; iTet++) {
+        for(i=0; i<4; i++) {
+          ix  = swapi4(tet[iTet][i]);
+          tet[iTet][i] = ix;
+        }
+      }
+    }
+    GmfSetKwd(OutMsh, GmfTetrahedra,nbrTet);
+    for (iTet=0; iTet<nbrQua; iTet++) {
+      GmfSetLin(OutMsh, GmfTetrahedra, tet[iTet][0], tet[iTet][1], tet[iTet][2],  tet[iTet][3], 0); 
+    }
+    if ( tet ) {free(tet); tet = NULL; }
+    
+    
+    nbrPyr  = swapi4(size[4]);
+    if ( nbrPyr > 0 ) {
+      pyr   = (int5 *)malloc(sizeof(mint)*5*nbrPyr);
+      fread(pyr,sizeof(mint),5*nbrPyr,f);
+      for(iPyr=0; iPyr<nbrPyr; iPyr++) {
+        for(i=0; i<5; i++) {
+          ix  = swapi4(pyr[iPyr][i]);
+          pyr[iPyr][i] = ix;
+        }
+      }
+    }
+    GmfSetKwd(OutMsh, GmfPyramids,nbrPyr);
+    for (iPyr=0; iPyr<nbrQua; iPyr++) {
+      GmfSetLin(OutMsh, GmfPyramids, pyr[iPyr][0], pyr[iPyr][1], pyr[iPyr][2],  pyr[iPyr][3], pyr[iPyr][4] , 0); 
+    }
+    if ( pyr ) {free(pyr); pyr = NULL; }
+    
+    
+    nbrPri  = swapi4(size[5]);
+    if ( nbrPri > 0 ) {
+      pri   = (int6 *)malloc(sizeof(mint)*6*nbrPri);
+      fread(pri,sizeof(mint),6*nbrPri,f);
+      for(iPri=0; iPri<nbrPri; iPri++) {
+        for(i=0; i<6; i++) {
+          ix  = swapi4(pri[iPri][i]);
+          pri[iPri][i] = ix;
+        }
+      }
+    }
+    GmfSetKwd(OutMsh, GmfPrisms,nbrPri);
+    for (iPri=0; iPri<nbrPri; iPri++) {
+      GmfSetLin(OutMsh, GmfPrisms, pri[iPri][0], pri[iPri][1], pri[iPri][2],  pri[iPri][3], pri[iPri][4] , pri[iPri][5], 0); 
+    }
+    if ( pri ) {free(pri); pri = NULL; }
+    
+    
+    nbrHex  = swapi4(size[6]);
+    if ( nbrHex > 0 ) {
+      hex   = (int8 *)malloc(sizeof(mint)*8*nbrHex);
+      fread(hex,sizeof(mint),8*nbrHex,f);
+      for(iHex=0; iHex<nbrHex; iHex++) {
+        for(i=0; i<8; i++) {
+          ix  = swapi4(hex[iHex][i]);
+          hex[iHex][i] = ix;
+        }
+      }
+    }
+    GmfSetKwd(OutMsh, GmfHexahedra,nbrHex);
+    for (iHex=0; iHex<nbrQua; iHex++) {
+      GmfSetLin(OutMsh, GmfHexahedra, hex[iHex][0], hex[iHex][1], hex[iHex][2],  hex[iHex][3], hex[iHex][4], hex[iHex][5], hex[iHex][6],  hex[iHex][7], 0); 
+    }
+    if ( hex ) {free(hex); hex = NULL; }
+    
+  }
+  
+  GmfCloseMesh(OutMsh);
+  
+  
+  if ( crd ) { free(crd); crd = NULL; }
+  if ( tri ) { free(tri); tri = NULL; }
+  if ( qua ) { free(qua); qua = NULL; }
+  if ( tet ) { free(tet); tet = NULL; }
+  if ( pyr ) { free(pyr); pyr = NULL; }
+  if ( pri ) { free(pri); pri = NULL; }
+  if ( hex ) { free(hex); hex = NULL; }
+  
+  fclose(f);
+  
+  
+  return 0;
+}
diff --git a/whatsnew.md b/whatsnew.md
index dbe806b..1fb0960 100644
--- a/whatsnew.md
+++ b/whatsnew.md
@@ -1,12 +1,8 @@
-## Release 7.60
+## Release 7.62
 
-1. Arbitrary polyhedra and polygons storage is now possible with these six new keywords:
-  - GmfBoundaryPolygonHeaders: stores the number and the list nodes of each boundary polygon
-  - GmfBoundaryPolygonVertices: stores the concatenated list of boundary polygons' nodes
-  - GmfInnerPolygonHeaders: stores the number and the list nodes of each inner polygon (polyhedra's faces)
-  - GmfInnerPolygonVertices: stores the concatenated list of inner polygons' nodes
-  - GmfPolyhedraHeaders: stores the number and the list of faces of each polyhedron
-  - GmfPolyhedraFaces: stores the concatenated list of polyhedra's faces
+1. Corrected two bugs:
+  - GmfOpenMesh: could crash with a stack overflow a present some security issues
+  - test\_libmeshb\_pipeline.f: the Fortran version of the user's procedure call was crashing
 
 2. New helpers functions system to easily add specific features related to the libMeshb:
   - See the helper's [readme](utilities/libmeshb7_helpers.md) for more information about the new functions to handle polyhedral meshes.