Skip to content
Snippets Groups Projects
Select Git revision
  • 2c73c7ebf1820678ed4e7cb89587ebfafb405f9e
  • develop default protected
  • feature/variational-hydro
  • origin/stage/bouguettaia
  • feature/gmsh-reader
  • feature/reconstruction
  • save_clemence
  • feature/kinetic-schemes
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

GmshReader.cpp

Blame
  • GmshReader.cpp 28.68 KiB
    #include <GmshReader.hpp>
    
    #include <iostream>
    #include <fstream>
    #include <set>
    #include <rang.hpp>
    
    #include <Connectivity2D.hpp>
    #include <Mesh.hpp>
    #include <map>
    
    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;
    
      __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]= false;
    
      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::__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 << " [IGNORED]";
          } 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
    	  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
          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) {
        std::cerr << "*** using a 3d mesh (NIY)\n";
        std::exit(0);
      } 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;
        }
        const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
        const size_t nb_triangles = __triangles.extent(0);
        for (size_t j=0; j<nb_triangles; ++j) {
          cell_nodes(j,0) = __triangles[j][0];
          cell_nodes(j,1) = __triangles[j][1];
          cell_nodes(j,2) = __triangles[j][2];
        }
        
        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_nodes(jq,0) = __quadrangles[j][0];
          cell_nodes(jq,1) = __quadrangles[j][1];
          cell_nodes(jq,2) = __quadrangles[j][2];
          cell_nodes(jq,3) = __quadrangles[j][3];
        }
        const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
        for (size_t j=0; j<nb_triangles; ++j) {
          cell_nb_nodes[j] = 3;
        }
        for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
          cell_nb_nodes[j] = 4;
        }
        m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
        Connectivity2D& connectivity = *m_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) {
          connectivity.addFaceBoundary(ref_face_list.first,
        				   ref_face_list.second);
        }
    
        const Kokkos::View<const unsigned int**> face_nodes = connectivity.faceNodes();
    
        for (size_t i=0; i<connectivity.numberOfFaceBoundaries(); ++i) {
          const Connectivity2D::FacesBoundary& faces_boundary
    	= connectivity.facesBoundary(i);
    
          const unsigned int ref = faces_boundary.first;
          std::set<unsigned int> node_id_set;
          const Kokkos::View<const unsigned int*> face_ids
    	= faces_boundary.second;
          
          for (unsigned int l=0; l<face_ids.extent(0); ++l) {
    	for (unsigned short r=0; r<2; ++r) {
    	  node_id_set.insert(face_nodes(face_ids[l],r));
    	}
          }
    
          Kokkos::View<unsigned int*> node_id_list("node_ids", node_id_set.size());
          {
    	int r=0;
    	for (auto node_id : node_id_set) {
    	  node_id_list[r] = node_id;
    	  ++r;
    	}
          }
          connectivity.addNodeBoundary(ref, node_id_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 = new MeshType(connectivity, xr);
        MeshType& mesh = *m_mesh;
    
        std::ofstream gnuplot("mesh.gnu");
    
        for (size_t j=0; j<mesh.numberOfCells(); ++j) {
          for (int r=0; r<mesh.connectivity().cellNbNodes()[j]; ++r) {
    	const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,r)];
    	gnuplot << x[0] << ' ' << x[1] << '\n';
          }
          const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,0)];
          gnuplot << x[0] << ' ' << x[1] << "\n\n";
        }
        gnuplot.close();
    
      } else if ((dimension1_mask, elementNumber)>0) {
        std::cerr << "*** using a 1d mesh (NIY)\n";
        std::exit(0);
      } else {
        std::cerr << "*** using a dimension 0 mesh is forbidden!\n";
        std::exit(0);
      }
    }
    
    GmshReader::Keyword
    GmshReader::__nextKeyword()
    {
      GmshReader::Keyword kw;
    
      std::cerr << "warning: " << rang::fg::red << __PRETTY_FUNCTION__ << rang::fg::reset << " keyword validity not checked!\n";
      int retval = 1;
      std::string aKeyword;
      m_fin >> aKeyword;
      if (retval < 0)  {
        kw.second = EndOfFile;
        return kw;
      } else if (retval == 0) {
        kw.second = Unknown;
        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;
      kw = this->__nextKeyword();
      if (kw.second != NODES) {
        // throw ErrorHandler(__FILE__,__LINE__,
        // 		       "reading file '"+m_filename
        // 		       +"': expecting $Nodes, '"+kw.first+"' was found",
        // 		       ErrorHandler::normal);
      }
    
      this->__readVertices();
    
      kw = this->__nextKeyword();
      if (kw.second != ENDNODES) {
        // throw ErrorHandler(__FILE__,__LINE__,
        // 		       "reading file '"+m_filename
        // 		       +"': expecting $EndNodes, '"+kw.first+"' was found",
        // 		       ErrorHandler::normal);
      }
    
      // Getting elements list
      kw = this->__nextKeyword();
      if (kw.second != ELEMENTS) {
        // throw ErrorHandler(__FILE__,__LINE__,
        // 		       "reading file '"+m_filename
        // 		       +"': expecting $Elements, '"+kw.first+"' was found",
        // 		       ErrorHandler::normal);
      }
    
      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);
      }
    }