GmshReader.cpp

#include <GmshReader.hpp>

#include <iostream>
#include <fstream>
#include <set>
#include <rang.hpp>

#include <Connectivity1D.hpp>
#include <Connectivity2D.hpp>
#include <Connectivity3D.hpp>

#include <Mesh.hpp>

#include <RefFaceList.hpp>

#include <map>
#include <regex>
#include <iomanip>

template<typename T>
inline std::string stringify(const T & t)
{
  std::ostringstream oss;
  oss << t;
  return oss.str();
}

template<>
inline std::string stringify<std::string>(const std::string& t)
{
  return t;
}

class ErrorHandler
{
 public:
  enum Type {
    asked,     /**< execution request by the user*/
    compilation, /**< syntax error in a language */
    normal,    /**< normal error due to a bad use of ff3d */
    unexpected /**< Unexpected execution error */
  };

 private:
  const std::string __filename;	/**< The source file name where the error occured */
  const size_t __lineNumber;	/**< The line number where exception was raised */
  const std::string __errorMessage; /**< The reporting message */

  const Type __type;		/**< the type of the error */
 public:
  /**
   * Prints the error message
   *
   */
  virtual void writeErrorMessage();

  /**
   * The copy constructor
   *
   * @param e an handled error
   */
  ErrorHandler(const ErrorHandler& e)
      : __filename(e.__filename),
        __lineNumber(e.__lineNumber),
        __errorMessage(e.__errorMessage),
        __type(e.__type)
  {
    ;
  }

  /**
   * Constructor
   *
   * @param filename the source file name
   * @param lineNumber the source line
   * @param errorMessage the reported message
   * @param type the type of the error
   */
  ErrorHandler(const std::string& filename,
               const size_t& lineNumber,
               const std::string& errorMessage,
               const Type& type);

  /**
   * The destructor
   *
   */
  virtual ~ErrorHandler()
  {
    ;
  }
};
void ErrorHandler::writeErrorMessage()
{
  switch(__type) {
    case asked: {
      std::cerr << "\nremark: exit command explicitly called\n";
    }
    case normal: {
      std::cerr << '\n' << __filename << ':' << __lineNumber
                << ":remark: emitted the following message\n";
      std::cerr << "error: " << __errorMessage << '\n';
      break;
    }
    case compilation: {
      std::cerr << "\nline " << __lineNumber << ':' << __errorMessage << '\n';
      break;
    }
    case unexpected: {
      std::cerr << '\n' << __filename << ':' << __lineNumber << ":\n" << __errorMessage << '\n';
      std::cerr << "\nUNEXPECTED ERROR: this should not occure, please report it\n";
      std::cerr << "\nBUG REPORT: Please send bug reports to:\n"
                << "  ff3d-dev@nongnu.org or freefem@ann.jussieu.fr\n"
                << "or better, use the Bug Tracking System:\n"
                << "  http://savannah.nongnu.org/bugs/?group=ff3d\n";
      break;
    }
    default: {
      std::cerr << __filename << ':' << __lineNumber << ": " << __errorMessage << '\n';
      std::cerr << __FILE__ << ':' << __LINE__ << ":remark:  error type not implemented!\n";
    }
  }
}

ErrorHandler::
ErrorHandler(const std::string& filename,
             const size_t& lineNumber,
             const std::string& errorMessage,
             const Type& type)
    : __filename(filename),
      __lineNumber(lineNumber),
      __errorMessage(errorMessage),
      __type(type)
{
  ;
}


GmshReader::GmshReader(const std::string& filename)
    : m_filename(filename)
{
  try {
    m_fin.open(m_filename);
    if (not m_fin) {
      std::cerr << "cannot read file '" << m_filename << "'\n";
      std::exit(0);
    }

    // gmsh 2.2 format keywords
    __keywordList["$MeshFormat"]       = MESHFORMAT;
    __keywordList["$EndMeshFormat"]    = ENDMESHFORMAT;
    __keywordList["$Nodes"]            = NODES;
    __keywordList["$EndNodes"]         = ENDNODES;
    __keywordList["$Elements"]         = ELEMENTS;
    __keywordList["$EndElements"]      = ENDELEMENTS;
    __keywordList["$PhysicalNames"]    = PHYSICALNAMES;
    __keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;

    __numberOfPrimitiveNodes.resize(16);
    __numberOfPrimitiveNodes[ 0] =  2; // edge
    __numberOfPrimitiveNodes[ 1] =  3; // triangle
    __numberOfPrimitiveNodes[ 2] =  4; // quadrangle
    __numberOfPrimitiveNodes[ 3] =  4; // Tetrahedron
    __numberOfPrimitiveNodes[ 4] =  8; // Hexaredron
    __numberOfPrimitiveNodes[ 5] =  6; // Prism
    __numberOfPrimitiveNodes[ 6] =  5; // Pyramid
    __numberOfPrimitiveNodes[ 7] =  3; // second order edge
    __numberOfPrimitiveNodes[ 8] =  6; // second order triangle
    __numberOfPrimitiveNodes[ 9] =  9; // second order quadrangle
    __numberOfPrimitiveNodes[10] = 10; // second order tetrahedron
    __numberOfPrimitiveNodes[11] = 27; // second order hexahedron
    __numberOfPrimitiveNodes[12] = 18; // second order prism
    __numberOfPrimitiveNodes[13] = 14; // second order pyramid
    __numberOfPrimitiveNodes[14] =  1; // point

    __primitivesNames[0]     = "edges";
    __supportedPrimitives[0] = true;
    __primitivesNames[1]     = "triangles";
    __supportedPrimitives[1] = true;
    __primitivesNames[2]     = "quadrangles";
    __supportedPrimitives[2] = true;
    __primitivesNames[3]     = "tetrahedra";
    __supportedPrimitives[3] = true;
    __primitivesNames[4]     = "hexahedra";
    __supportedPrimitives[4] = true;
    __primitivesNames[5]     = "prisms";
    __supportedPrimitives[5] = false;
    __primitivesNames[6]     = "pyramids";
    __supportedPrimitives[6] = false;
    __primitivesNames[7]     = "second order edges";
    __supportedPrimitives[7] = false;
    __primitivesNames[8]     = "second order triangles";
    __supportedPrimitives[8] = false;
    __primitivesNames[9]     = "second order quadrangles";
    __supportedPrimitives[9] = false;
    __primitivesNames[10]    = "second order tetrahedra";
    __supportedPrimitives[10]= false;
    __primitivesNames[11]    = "second order hexahedra";
    __supportedPrimitives[11]= false;
    __primitivesNames[12]    = "second order prisms";
    __supportedPrimitives[12]= false;
    __primitivesNames[13]    = "second order pyramids";
    __supportedPrimitives[13]= false;
    __primitivesNames[14]    = "point";
    __supportedPrimitives[14]= true;

    std::cout << "Reading file '" << m_filename << "'\n";

    // Getting vertices list
    GmshReader::Keyword kw = this->__nextKeyword();
    switch(kw.second) {
      // case NOD: {
      //   this->__readGmsh1();
      //   break;
      // }
      case MESHFORMAT: {
        double fileVersion = this->_getReal();
        if (fileVersion != 2.2) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "Cannot read Gmsh format '"+stringify(fileVersion)+"'",
                             ErrorHandler::normal);
        }
        int fileType = this->_getInteger();
        __binary = (fileType == 1);
        if ((fileType < 0) or (fileType > 1)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "Cannot read Gmsh file type '"+stringify(fileType)+"'",
                             ErrorHandler::normal);
        }

        int dataSize = this->_getInteger();
        if (dataSize != sizeof(double)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "Data size not supported '"+stringify(dataSize)+"'",
                             ErrorHandler::normal);
        }

        if (__binary) {
          //       int one=0;

          //       fseek(__ifh,1L,SEEK_CUR);
          //       fread(reinterpret_cast<char*>(&one),sizeof(int),1,__ifh);

          //       if (one == 1) {
          // 	__convertEndian = false;
          //       } else {
          // 	invertEndianess(one);

          // 	if (one == 1) {
          // 	  __convertEndian = true;
          // 	} else {
          // 	  throw ErrorHandler(__FILE__,__LINE__,
          // 	  		     "Cannot determine endianess",
          // 	  		     ErrorHandler::normal);
          // 	}
          // }

          //       std::cout << "- Binary format: ";
          // #ifdef    WORDS_BIGENDIAN
          //       if (not __convertEndian) {
          // 	std::cout << "big endian\n";
          //       } else {
          // 	std::cout << "little endian\n";
          //       }
          // #else  // WORDS_BIGENDIAN
          //       if (not __convertEndian) {
          // 	std::cout << "little endian\n";
          //       } else {
          // 	std::cout << "big endian\n";
          //       }
          // #endif // WORDS_BIGENDIAN
        }

        kw = this->__nextKeyword();
        if (kw.second != ENDMESHFORMAT) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': expecting $EndMeshFormat, '"+kw.first+"' was found",
                             ErrorHandler::normal);
        }

        this->__readGmshFormat2_2();

        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "cannot determine format version of '"+m_filename+"'",
                           ErrorHandler::normal);
      }
    }

    this->__proceedData();
    // this->__createMesh();
  }
  catch(ErrorHandler e) {
    e.writeErrorMessage();
    std::exit(0);
  }
}

void GmshReader::__readVertices()
{
  const int numberOfVerices = this->_getInteger();
  std::cout << "- Number Of Vertices: " << numberOfVerices << '\n';
  if (numberOfVerices<0) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "reading file '"+this->m_filename
                       +"': number of vertices is negative",
                       ErrorHandler::normal);
  }

  __verticesNumbers.resize(numberOfVerices);
  __vertices = Kokkos::View<R3*>("vertices",numberOfVerices);

  if (not __binary) {
    for (int i=0; i<numberOfVerices; ++i) {
      __verticesNumbers[i] = this->_getInteger();
      const double x = this->_getReal();
      const double y = this->_getReal();
      const double z = this->_getReal();
      __vertices[i] = TinyVector<3, double>(x,y,z);
    }
  } else {
    // fseek(m_fin.file()__ifh,1L,SEEK_CUR);
    // for (int i=0; i<numberOfVerices; ++i) {
    //   int number = 0;
    //   fread(reinterpret_cast<char*>(&number),sizeof(int),1,__ifh);
    //   __verticesNumbers[i] = number;
    //   TinyVector<3,double> x;
    //   fread(reinterpret_cast<char*>(&(x[0])),sizeof(double),3,__ifh);
    //   for (size_t j=0; j<3; ++j) {
    // 	__vertices[i][j] = x[j];
    //   }
    // }

    // if (__convertEndian) {
    //   for (int i=0; i<numberOfVerices; ++i) {
    // 	invertEndianess(__verticesNumbers[i]);
    // 	for (size_t j=0; j<3; ++j) {
    // 	  invertEndianess(vertices[i][j]);
    // 	}
    //   }
    // }
  }
}

// void
// GmshReader::__readElements1()
// {
//   const int numberOfElements = this->_getInteger();
//   std::cout << "- Number Of Elements: " << numberOfElements << '\n';
//   if (numberOfElements<0) {
//     throw ErrorHandler(__FILE__,__LINE__,
// 		       "reading file '"+m_filename
// 		       +"': number of elements is negative",
// 		       ErrorHandler::normal);
//   }

//   __elementType.resize(numberOfElements);
//   __references.resize(numberOfElements);
//   __elementVertices.resize(numberOfElements);

//   for (int i=0; i<numberOfElements; ++i) {
//     this->_getInteger(); // drops element number
//     const int elementType = this->_getInteger()-1;

//     if ((elementType < 0) or (elementType > 14)) {
//       throw ErrorHandler(__FILE__,__LINE__,
// 			 "reading file '"+m_filename
// 			 +"': unknown element type '"+stringify(elementType)+"'",
// 			 ErrorHandler::normal);
//     }
//     __elementType[i] = elementType;
//     __references[i] = this->_getInteger(); // physical reference
//     this->_getInteger(); // drops entity reference

//     const int numberOfVertices = this->_getInteger();
//     if (numberOfVertices != __numberOfPrimitiveNodes[elementType]) {
//       std::stringstream errorMsg;
//       errorMsg << "reading file '" <<m_filename
// 	       << "':number of vertices (" << numberOfVertices
// 	       << ") does not match expected ("
// 	       << __numberOfPrimitiveNodes[elementType] << ")" << std::ends;

//       throw ErrorHandler(__FILE__,__LINE__,
// 			 errorMsg.str(),
// 			 ErrorHandler::normal);
//     }
//     __elementVertices[i].resize(numberOfVertices);
//     for (int j=0; j<numberOfVertices; ++j) {
//       __elementVertices[i][j] = this->_getInteger();
//     }
//   }
// }

void
GmshReader::__readElements2_2()
{
  const int numberOfElements = this->_getInteger();
  std::cout << "- Number Of Elements: " << numberOfElements << '\n';
  if (numberOfElements<0) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "reading file '"+m_filename
                       +"': number of elements is negative",
                       ErrorHandler::normal);
  }

  __elementType.resize(numberOfElements);
  __references.resize(numberOfElements);
  __elementVertices.resize(numberOfElements);

  if (not __binary) {
    for (int i=0; i<numberOfElements; ++i) {
      this->_getInteger(); // drops element number
      const int elementType = this->_getInteger()-1;

      if ((elementType < 0) or (elementType > 14)) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "reading file '"+m_filename
                           +"': unknown element type '"+stringify(elementType)+"'",
                           ErrorHandler::normal);
      }
      __elementType[i] = elementType;
      const int numberOfTags =  this->_getInteger();
      __references[i] = this->_getInteger(); // physical reference
      for (int tag=1; tag<numberOfTags; ++tag) {
        this->_getInteger(); // drops remaining tags
      }

      const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
      __elementVertices[i].resize(numberOfVertices);
      for (int j=0; j<numberOfVertices; ++j) {
        __elementVertices[i][j] = this->_getInteger();
      }
    }
  } else {
    // fseek(__ifh,1L,SEEK_CUR);
    // int i=0;
    // for (;i<numberOfElements;) {
    //   int elementType = 0;
    //   fread(reinterpret_cast<char*>(&elementType),sizeof(int),1,__ifh);
    //   if (__convertEndian) {
    // 	invertEndianess(elementType);
    //   }
    //   --elementType;

    //   if ((elementType < 0) or (elementType > 14)) {
    // 	 throw ErrorHandler(__FILE__,__LINE__,
    // 	 		   "reading file '"+m_filename
    // 	 		   +"': unknown element type '"+stringify(elementType)+"'",
    // 	 		   ErrorHandler::normal);
    //   }

    //   int elementTypeNumber = 0;
    //   fread(reinterpret_cast<char*>(&elementTypeNumber),sizeof(int),1,__ifh);
    //   if (__convertEndian) {
    // 	invertEndianess(elementTypeNumber);
    //   }

    //   int numberOfTags = 0;
    //   fread(reinterpret_cast<char*>(&numberOfTags),sizeof(int),1,__ifh);
    //   if (__convertEndian) {
    // 	invertEndianess(numberOfTags);
    //   }

    //   for (int k=0; k<elementTypeNumber; ++k) {
    // 	__elementType[i] = elementType;

    // 	int elementNumber = 0;
    // 	fread(reinterpret_cast<char*>(&elementNumber),sizeof(int),1,__ifh);

    // 	int reference = 0;
    // 	fread(reinterpret_cast<char*>(&reference),sizeof(int),1,__ifh);

    // 	__references[i] = reference; // physical reference
    // 	for (int tag=1; tag<numberOfTags; ++tag) {
    // 	  int t;
    // 	  fread(reinterpret_cast<char*>(&t),sizeof(int),1,__ifh);
    // 	}

    // 	const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
    // 	__elementVertices[i].resize(numberOfVertices);
    // 	fread(reinterpret_cast<char*>(&(__elementVertices[i][0])),sizeof(int),numberOfVertices,__ifh);
    // 	++i;
    //   }
    // }

    // if (__convertEndian) {
    //   for (size_t i=0; i<__references.size(); ++i) {
    // 	invertEndianess(__references[i]);
    //   }
    //   for (size_t i=0; i<__elementVertices.size(); ++i) {
    // 	for (size_t j=0; j<__elementVertices[i].size(); ++i) {
    // 	  invertEndianess(__elementVertices[i][j]);
    // 	}
    //   }
    // }
  }
}

void
GmshReader::
__readPhysicalNames2_2()
{
  const int number_of_names = this->_getInteger();
  for (int i=0; i<number_of_names; ++i) {
    const int physical_dimension = this->_getInteger();
    const int physical_number = this->_getInteger();
    std::string physical_name;
    m_fin >> physical_name;
    physical_name = std::regex_replace(physical_name, std::regex("(\")"), "");

    PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
    auto searched_physical_ref_id = m_physical_ref_map.find(physical_number);
    if (searched_physical_ref_id != m_physical_ref_map.end()) {
      std::cerr << "Physical reference id '" << physical_ref_id
                << "' already defined as '" << searched_physical_ref_id->second << "'!";
      std::exit(0);
    }
    m_physical_ref_map[physical_number] = physical_ref_id;
  }
}

void
GmshReader::__proceedData()
{
  TinyVector<15,int> elementNumber = zero;
  std::vector<std::set<size_t> > elementReferences(15);

  for (size_t i=0; i<__elementType.size(); ++i) {
    const int elementType = __elementType[i];
    elementNumber[elementType]++;
    elementReferences[elementType].insert(__references[i]);
  }

  for (size_t i=0; i<elementNumber.dimension(); ++i) {
    if (elementNumber[i]>0) {
      std::cout << "  - Number of "
                << __primitivesNames[i]
                << ": " << elementNumber[i];
      if (not(__supportedPrimitives[i])) {
        std::cout << " [" << rang::fg::yellow << "IGNORED" << rang::style::reset << "]";
      } else {
        std::cout << " references: ";
        for(std::set<size_t>::const_iterator iref
                = elementReferences[i].begin();
            iref != elementReferences[i].end(); ++iref) {
          if (iref != elementReferences[i].begin()) std::cout << ',';
          std::cout << *iref;
        }

        switch (i) {
          // Supported entities
          case 0: { // edges
            __edges = Kokkos::View<Edge*>("edges",
                                          elementNumber[i]);
            __edges_ref.resize(elementNumber[i]);
            break;
          }
          case 1: { // triangles
            __triangles = Kokkos::View<Triangle*>("triangles",
                                                  elementNumber[i]);
            __triangles_ref.resize(elementNumber[i]);
            break;
          }
          case  2: { // quadrangles
            __quadrangles = Kokkos::View<Quadrangle*>("quadrangles",
                                                      elementNumber[i]);
            __quadrangles_ref.resize(elementNumber[i]);
            break;
          }
          case 3: { // tetrahedra
            __tetrahedra = Kokkos::View<Tetrahedron*>("tetrahedra",
                                                      elementNumber[i]);
            __tetrahedra_ref.resize(elementNumber[i]);
            break;
          }
          case  4: {// hexahedra
            __hexahedra = Kokkos::View<Hexahedron*>("hexahedra",
                                                    elementNumber[i]);
            __hexahedra_ref.resize(elementNumber[i]);
          }
            // Ignored entities
          case 14: {// point
            __points = Kokkos::View<Point*>("points",
                                            elementNumber[i]);
            __points_ref.resize(elementNumber[i]);
            break;
          }
            // Unsupported entities
          case  5: // prism
          case  6: // pyramid
          case  7: // second order edge
          case  8: // second order triangle
          case  9: // second order quadrangle
          case 10: // second order tetrahedron
          case 11: // second order hexahedron
          case 12: // second order prism
          case 13: // second order pyramid
          default: {
            throw ErrorHandler(__FILE__,__LINE__,
                               "reading file '"+m_filename
                               +"': ff3d cannot read this kind of element",
                               ErrorHandler::normal);
          }
        }
      }
      std::cout << '\n';
    }
  }

  std::cout << "- Building correspondance table\n";
  int minNumber = std::numeric_limits<int>::max();
  int maxNumber = std::numeric_limits<int>::min();
  for (size_t i=0; i<__verticesNumbers.size(); ++i) {
    const int& vertexNumber = __verticesNumbers[i];
    minNumber = std::min(minNumber,vertexNumber);
    maxNumber = std::max(maxNumber,vertexNumber);
  }

  if (minNumber<0) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "reading file '"+m_filename
                       +"': ff3d does not implement negative vertices number",
                       ErrorHandler::normal);
  }

  // A value of -1 means that the vertex is unknown
  __verticesCorrepondance.resize(maxNumber+1,-1);

  for (size_t i=0; i<__verticesNumbers.size(); ++i) {

    __verticesCorrepondance[__verticesNumbers[i]] = i;
  }

  // reset element number to count them while filling structures
  elementNumber = zero;

  // Element structures filling
  for (size_t i=0; i<__elementType.size(); ++i) {
    std::vector<int>& elementVertices = __elementVertices[i];
    switch (__elementType[i]) {
      // Supported entities
      case  0: { // edge
        int& edgeNumber = elementNumber[0];
        const int a = __verticesCorrepondance[elementVertices[0]];
        const int b = __verticesCorrepondance[elementVertices[1]];
        if ((a<0) or (b<0)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': error reading element "+stringify(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __edges[edgeNumber]
            = Edge(a, b);
        __edges_ref[edgeNumber] = __references[i];
        edgeNumber++; // one more edge
        break;
      }
      case 1: { // triangles
        int& triangleNumber = elementNumber[1];

        const int a = __verticesCorrepondance[elementVertices[0]];
        const int b = __verticesCorrepondance[elementVertices[1]];
        const int c = __verticesCorrepondance[elementVertices[2]];
        if ((a<0) or (b<0) or (c<0)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': error reading element "+stringify(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __triangles[triangleNumber]
            = Triangle(a, b, c);
        __triangles_ref[triangleNumber] = __references[i];
        triangleNumber++; // one more triangle
        break;
      }
      case 2: { // quadrangle
        int& quadrilateralNumber = elementNumber[2];

        const int a = __verticesCorrepondance[elementVertices[0]];
        const int b = __verticesCorrepondance[elementVertices[1]];
        const int c = __verticesCorrepondance[elementVertices[2]];
        const int d = __verticesCorrepondance[elementVertices[3]];
        if ((a<0) or (b<0) or (c<0) or (d<0)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': error reading element "+stringify(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __quadrangles[quadrilateralNumber] = Quadrangle(a,b,c,d);
        __quadrangles_ref[quadrilateralNumber] = __references[i];
        quadrilateralNumber++; // one more quadrangle
        break;
      }
      case 3: { // tetrahedra
        int& tetrahedronNumber = elementNumber[3];

        const int a = __verticesCorrepondance[elementVertices[0]];
        const int b = __verticesCorrepondance[elementVertices[1]];
        const int c = __verticesCorrepondance[elementVertices[2]];
        const int d = __verticesCorrepondance[elementVertices[3]];
        if ((a<0) or (b<0) or (c<0) or (d<0)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': error reading element "+stringify(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __tetrahedra[tetrahedronNumber] = Tetrahedron(a, b, c, d);
        __tetrahedra_ref[tetrahedronNumber] = __references[i];
        tetrahedronNumber++; // one more tetrahedron
        break;
      }
      case  4: { // hexaredron
        int& hexahedronNumber = elementNumber[4];

        const int a = __verticesCorrepondance[elementVertices[0]];
        const int b = __verticesCorrepondance[elementVertices[1]];
        const int c = __verticesCorrepondance[elementVertices[2]];
        const int d = __verticesCorrepondance[elementVertices[3]];
        const int e = __verticesCorrepondance[elementVertices[4]];
        const int f = __verticesCorrepondance[elementVertices[5]];
        const int g = __verticesCorrepondance[elementVertices[6]];
        const int h = __verticesCorrepondance[elementVertices[7]];
        if ((a<0) or (b<0) or (c<0) or (d<0) or (e<0) or (f<0) or (g<0) or (h<0)) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': error reading element "+stringify(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __hexahedra[hexahedronNumber]
            = Hexahedron(a, b, c, d,
                         e, f, g, h);
        __hexahedra_ref[hexahedronNumber] = __references[i];
        hexahedronNumber++; // one more hexahedron
        break;
      }
        // Unsupported entities
      case 14:{// point
        int& point_number = elementNumber[14];
        const int a = __verticesCorrepondance[elementVertices[0]];
        __points[point_number] = a;
        __points_ref[point_number] = __references[i];
        point_number++;
        break;
      }
      case  5: // prism
      case  6: // pyramid
      case  7: // second order edge
      case  8: // second order triangle
      case  9: // second order quadrangle
      case 10: // second order tetrahedron
      case 11: // second order hexahedron
      case 12: // second order prism
      case 13:{// second order pyramid
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "reading file '"+m_filename
                           +"': ff3d cannot read this kind of element",
                           ErrorHandler::normal);
      }
    }
  }

  TinyVector<15,int> dimension0_mask = zero;
  dimension0_mask[14]=1;
  TinyVector<15,int> dimension1_mask = zero;
  dimension1_mask[0]=1;
  dimension1_mask[7]=1;
  TinyVector<15,int> dimension2_mask = zero;
  dimension2_mask[1]=1;
  dimension2_mask[2]=1;
  dimension2_mask[8]=1;
  dimension2_mask[9]=1;
  TinyVector<15,int> dimension3_mask = zero;
  dimension3_mask[3]=1;
  dimension3_mask[4]=1;
  dimension3_mask[5]=1;
  dimension3_mask[6]=1;
  dimension3_mask[10]=1;
  dimension3_mask[11]=1;
  dimension3_mask[12]=1;
  dimension3_mask[13]=1;

  std::cout << "- dimension 0 entities : " << (dimension0_mask, elementNumber) << '\n';
  std::cout << "- dimension 1 entities : " << (dimension1_mask, elementNumber) << '\n';
  std::cout << "- dimension 2 entities : " << (dimension2_mask, elementNumber) << '\n';
  std::cout << "- dimension 3 entities : " << (dimension3_mask, elementNumber) << '\n';
  if ((dimension3_mask, elementNumber)>0) {
    const size_t nb_cells = (dimension3_mask, elementNumber);
    size_t max_nb_node_per_cell=4;
    if (elementNumber[4] > 0) {
      max_nb_node_per_cell = 8;
    }

    std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
    const size_t nb_tetrahedra = __tetrahedra.extent(0);
    for (size_t j=0; j<nb_tetrahedra; ++j) {
      cell_by_node_vector[j].resize(4);
      for (int r=0; r<4; ++r) {
        cell_by_node_vector[j][r] = __tetrahedra[j][r];
      }
    }
    const size_t nb_hexahedra = __hexahedra.extent(0);
    for (size_t j=0; j<nb_hexahedra; ++j) {
      const size_t jh = nb_tetrahedra+j;
      cell_by_node_vector[jh].resize(8);
      for (int r=0; r<8; ++r) {
        cell_by_node_vector[jh][r] = __hexahedra[j][r];
      }
    }

    std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector));
    Connectivity3D& connectivity = *p_connectivity;

    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
    for (unsigned int f=0; f<__triangles.extent(0); ++f) {
      const unsigned int face_number
          = connectivity.getFaceNumber({__triangles[f][0], __triangles[f][1], __triangles[f][2]});
      const unsigned int& ref = __triangles_ref[f];
      ref_faces_map[ref].push_back(face_number);
    }
    for (unsigned int f=0; f<__quadrangles.extent(0); ++f) {
      const unsigned int face_number
          = connectivity.getFaceNumber({__quadrangles[f][0], __quadrangles[f][1],
                                        __quadrangles[f][2], __quadrangles[f][3]});
      const unsigned int& ref = __quadrangles_ref[f];
      ref_faces_map[ref].push_back(face_number);
    }
    for (const auto& ref_face_list : ref_faces_map) {
      Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
      for (size_t j=0; j<ref_face_list.second.size(); ++j) {
        face_list[j]=ref_face_list.second[j];
      }
      const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
      connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
    }

    typedef Mesh<Connectivity3D> MeshType;
    typedef TinyVector<3, double> Rd;

    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
    for (size_t i=0; i<__vertices.extent(0); ++i) {
      xr[i][0] = __vertices[i][0];
      xr[i][1] = __vertices[i][1];
      xr[i][2] = __vertices[i][2];
    }
    m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));

  } else if ((dimension2_mask, elementNumber)>0) {
    const size_t nb_cells = (dimension2_mask, elementNumber);
    size_t max_nb_node_per_cell=3;
    if (elementNumber[2] > 0) {
      max_nb_node_per_cell = 4;
    }

    std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
    const size_t nb_triangles = __triangles.extent(0);
    for (size_t j=0; j<nb_triangles; ++j) {
      cell_by_node_vector[j].resize(3);
      for (int r=0; r<3; ++r) {
        cell_by_node_vector[j][r] = __triangles[j][r];
      }
    }

    const size_t nb_quadrangles = __quadrangles.extent(0);
    for (size_t j=0; j<nb_quadrangles; ++j) {
      const size_t jq = j+nb_triangles;
      cell_by_node_vector[jq].resize(4);
      for (int r=0; r<4; ++r) {
        cell_by_node_vector[jq][r] = __quadrangles[j][r];
      }
    }

    std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector));
    Connectivity2D& connectivity = *p_connectivity;

    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
    for (unsigned int e=0; e<__edges.extent(0); ++e) {
      const unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
      const unsigned int& ref = __edges_ref[e];
      ref_faces_map[ref].push_back(edge_number);
    }

    for (const auto& ref_face_list : ref_faces_map) {
      Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
      for (size_t j=0; j<ref_face_list.second.size(); ++j) {
        face_list[j]=ref_face_list.second[j];
      }
      const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
      connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
    }

    std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
    for (unsigned int r=0; r<__points.extent(0); ++r) {
      const unsigned int point_number = __points[r];
      const unsigned int& ref = __points_ref[r];
      ref_points_map[ref].push_back(point_number);
    }

    for (const auto& ref_point_list : ref_points_map) {
      Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
      for (size_t j=0; j<ref_point_list.second.size(); ++j) {
        point_list[j]=ref_point_list.second[j];
      }
      const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
      connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
    }

    typedef Mesh<Connectivity2D> MeshType;
    typedef TinyVector<2, double> Rd;

    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
    for (size_t i=0; i<__vertices.extent(0); ++i) {
      xr[i][0] = __vertices[i][0];
      xr[i][1] = __vertices[i][1];
    }
    m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));

  } else if ((dimension1_mask, elementNumber)>0) {

    const size_t nb_cells = (dimension1_mask, elementNumber);

    std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
    for (size_t j=0; j<nb_cells; ++j) {
      cell_by_node_vector[j].resize(2);
      for (int r=0; r<2; ++r) {
        cell_by_node_vector[j][r] = __edges[j][r];
      }
    }

    std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector));
    Connectivity1D& connectivity = *p_connectivity;

    std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
    for (unsigned int r=0; r<__points.extent(0); ++r) {
      const unsigned int point_number = __points[r];
      const unsigned int& ref = __points_ref[r];
      ref_points_map[ref].push_back(point_number);
    }

    for (const auto& ref_point_list : ref_points_map) {
      Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
      for (size_t j=0; j<ref_point_list.second.size(); ++j) {
        point_list[j]=ref_point_list.second[j];
      }
      const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
      connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
    }

    typedef Mesh<Connectivity1D> MeshType;
    typedef TinyVector<1, double> Rd;

    Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
    for (size_t i=0; i<__vertices.extent(0); ++i) {
      xr[i][0] = __vertices[i][0];
    }

    m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));

  } else {
    std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
    std::exit(0);
  }
}

GmshReader::Keyword
GmshReader::__nextKeyword()
{
  GmshReader::Keyword kw;

  std::string aKeyword;
  m_fin >> aKeyword;
  if (not m_fin) {
    kw.second = EndOfFile;
    return kw;
  }

  KeywordList::iterator i = __keywordList.find(aKeyword.c_str());

  if (i != __keywordList.end()) {
    kw.first = (*i).first;
    kw.second = (*i).second;
    return kw;
  }

  throw ErrorHandler(__FILE__,__LINE__,
                     "reading file '"+m_filename
                     +"': unknown keyword '"+aKeyword+"'",
                     ErrorHandler::normal);

  kw.first = aKeyword;
  kw.second = Unknown;
  return kw;
}

void GmshReader::
__readGmshFormat2_2()
{
  std::cout << "- Reading Gmsh format 2.2\n";
  GmshReader::Keyword kw = std::make_pair("", Unknown);
  while (kw.second != EndOfFile) {
    kw = this->__nextKeyword();
    switch (kw.second) {
      case NODES: {
        this->__readVertices();
        if (this->__nextKeyword().second != ENDNODES) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': expecting $EndNodes, '"+kw.first+"' was found",
                             ErrorHandler::normal);
        }
        break;
      }
      case ELEMENTS: {
        this->__readElements2_2();
        kw = this->__nextKeyword();
        if (kw.second != ENDELEMENTS) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': expecting $EndElements, '"+kw.first+"' was found",
                             ErrorHandler::normal);
        }
        break;
      }
      case PHYSICALNAMES: {
        this->__readPhysicalNames2_2();
        if (this->__nextKeyword().second != ENDPHYSICALNAMES) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "reading file '"+m_filename
                             +"': expecting $EndNodes, '"+kw.first+"' was found",
                             ErrorHandler::normal);
        }
        break;
      }

      case EndOfFile: {
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "reading file '"+m_filename
                           +"': unexpected '"+kw.first+"'",
                           ErrorHandler::normal);
      }
    }
  }
}