GmshReader.cpp

#include <GmshReader.hpp>
#include <PastisMacros.hpp>

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

#include <CellType.hpp>
#include <Connectivity.hpp>

#include <Mesh.hpp>
#include <MeshData.hpp>

#include <RefFaceList.hpp>
#include <Messenger.hpp>
#include <Partitioner.hpp>

#include <ArrayUtils.hpp>

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

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() const;

  /**
   * 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() const
{
  switch(__type) {
    case asked: {
      perr() << "\nremark: exit command explicitly called\n";
      [[fallthrough]];
    }
    case normal: {
      perr() << '\n' << __filename << ':' << __lineNumber
                << ":remark: emitted the following message\n";
      perr() << "error: " << __errorMessage << '\n';
      break;
    }
    case compilation: {
      perr() << "\nline " << __lineNumber << ':' << __errorMessage << '\n';
      break;
    }
    case unexpected: {
      perr() << '\n' << __filename << ':' << __lineNumber << ":\n" << __errorMessage << '\n';
      perr() << "\nUNEXPECTED ERROR: this should not occure, please report it\n";
      perr() << "\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: {
      perr() << __filename << ':' << __lineNumber << ": " << __errorMessage << '\n';
      perr() << __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)
{
  ;
}

template <size_t Dimension>
class ConnectivityFace;

template<>
class ConnectivityFace<1>
{
 public:
  friend struct Hash;
  struct Hash
  {
    size_t operator()(const ConnectivityFace& f) const;
  };
};

template<>
class ConnectivityFace<2>
{
 public:
  friend struct Hash;
  struct Hash
  {
    size_t operator()(const ConnectivityFace& f) const {
      size_t hash = 0;
      hash ^= std::hash<unsigned int>()(f.m_node0_id);
      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
      return hash;
    }
  };

  unsigned int m_node0_id;
  unsigned int m_node1_id;

  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
  {
    os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
    return os;
  }

  PASTIS_INLINE
  bool operator==(const ConnectivityFace& f) const
  {
    return ((m_node0_id == f.m_node0_id) and
            (m_node1_id == f.m_node1_id));
  }

  PASTIS_INLINE
  bool operator<(const ConnectivityFace& f) const
  {
    return ((m_node0_id<f.m_node0_id) or
            ((m_node0_id == f.m_node0_id) and
             (m_node1_id<f.m_node1_id)));
  }

  PASTIS_INLINE
  ConnectivityFace& operator=(const ConnectivityFace&) = default;

  PASTIS_INLINE
  ConnectivityFace& operator=(ConnectivityFace&&) = default;

  PASTIS_INLINE
  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
  {
    Assert(given_node_id_list.size()==2);
#warning rework this dirty constructor
    const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]);
    m_node0_id = min;
    m_node1_id = max;
  }

  PASTIS_INLINE
  ConnectivityFace(const ConnectivityFace&) = default;

  PASTIS_INLINE
  ConnectivityFace(ConnectivityFace&&) = default;

  PASTIS_INLINE
  ~ConnectivityFace() = default;
};

template <>
class ConnectivityFace<3>
{
 private:
  friend class GmshReader;
  friend struct Hash;
  struct Hash
  {
    size_t operator()(const ConnectivityFace& f) const {
      size_t hash = 0;
      for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
      }
      return hash;
    }
  };

  bool m_reversed;
  std::vector<NodeId::base_type> m_node_id_list;

  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
  {
    for (auto id : f.m_node_id_list) {
      os << id << ' ';
    }
    return os;
  }

  PASTIS_INLINE
  const bool& reversed() const
  {
    return m_reversed;
  }

  PASTIS_INLINE
  const std::vector<unsigned int>& nodeIdList() const
  {
    return m_node_id_list;
  }

  PASTIS_INLINE
  std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
  {
    const auto min_id = std::min_element(node_list.begin(), node_list.end());
    const int shift = std::distance(node_list.begin(), min_id);

    std::vector<unsigned int> rotated_node_list(node_list.size());
    if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
      for (size_t i=0; i<node_list.size(); ++i) {
        rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
        m_reversed = true;
      }
    } else {
      for (size_t i=0; i<node_list.size(); ++i) {
        rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
      }
    }

    return rotated_node_list;
  }

  PASTIS_INLINE
  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
      : m_reversed(false),
        m_node_id_list(_sort(given_node_id_list))
  {
    ;
  }

 public:
  bool operator==(const ConnectivityFace& f) const
  {
    if (m_node_id_list.size() == f.nodeIdList().size()) {
      for (size_t j=0; j<m_node_id_list.size(); ++j) {
        if (m_node_id_list[j] != f.nodeIdList()[j]) {
          return false;
        }
      }
      return true;
    }
    return false;
  }

  PASTIS_INLINE
  bool operator<(const ConnectivityFace& f) const
  {
    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
    for (size_t i=0; i<min_nb_nodes; ++i) {
      if (m_node_id_list[i] <  f.m_node_id_list[i]) return true;
      if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
    }
    return m_node_id_list.size() < f.m_node_id_list.size();
  }

  PASTIS_INLINE
  ConnectivityFace& operator=(const ConnectivityFace&) = default;

  PASTIS_INLINE
  ConnectivityFace& operator=(ConnectivityFace&&) = default;

  PASTIS_INLINE
  ConnectivityFace(const ConnectivityFace&) = default;

  PASTIS_INLINE
  ConnectivityFace(ConnectivityFace&&) = default;


  PASTIS_INLINE
  ConnectivityFace() = delete;

  PASTIS_INLINE
  ~ConnectivityFace() = default;
};

template <int Dimension>
class MeshDispatcher
{
 public:
  using MeshType = Mesh<Connectivity<Dimension>>;

 private:
  const MeshType& m_mesh;
  CellValue<const int> m_cell_new_owner;
  NodeValue<const int> m_node_new_owner;

  using CellListToSendByProc = std::vector<Array<const CellId>>;
  const CellListToSendByProc m_cell_list_to_send_by_proc;

  Array<int> m_nb_cell_to_send_by_proc;
  Array<int> m_nb_cell_to_recv_by_proc;

  CellValue<int> _getCellNewOwner()
  {
    CSRGraph mesh_graph = m_mesh.cellToCellGraph();
    Partitioner P;

    CellValue<int> cell_new_owner(m_mesh.connectivity());
    cell_new_owner = P.partition(mesh_graph);
    return cell_new_owner;
  }

  NodeValue<int> _getNodeNewOwner()
  {
    const auto& node_to_cell_matrix
        = m_mesh.connectivity().nodeToCellMatrix();

    pout() << __FILE__ << ':' << __LINE__ << ": use Min function\n";
#warning could use a better policy
    NodeValue<int> node_new_owner(m_mesh.connectivity());
    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
        const auto& node_to_cell = node_to_cell_matrix[r];
        CellId Jmin = node_to_cell[0];

        for (size_t j=1; j<node_to_cell.size(); ++j) {
          const CellId J = node_to_cell[j];
          if (J<Jmin) {
            Jmin=J;
          }
        }
        node_new_owner[r] = m_cell_new_owner[Jmin];
      });

#warning Add missing synchronize
    return node_new_owner;
  }

  const CellListToSendByProc
  _buildCellListToSend() const
  {
    const auto& node_to_cell_matrix
        = m_mesh.connectivity().nodeToCellMatrix();
    const auto& cell_to_node_matrix
        = m_mesh.connectivity().cellToNodeMatrix();

    std::vector<std::vector<CellId>> cell_vector_to_send_by_proc(parallel::size());
    Array<bool> send_to_rank(parallel::size());
    for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
      send_to_rank.fill(false);
      const auto& cell_to_node = cell_to_node_matrix[j];

      for (size_t R=0; R<cell_to_node.size(); ++R) {
        const NodeId& r = cell_to_node[R];
        const auto& node_to_cell = node_to_cell_matrix[r];
        for (size_t K=0; K<node_to_cell.size(); ++K) {
          const CellId& k = node_to_cell[K];
          send_to_rank[m_cell_new_owner[k]] = true;
        }
      }

      for (size_t k=0; k<send_to_rank.size(); ++k) {
        if (send_to_rank[k]) {
          cell_vector_to_send_by_proc[k].push_back(j);
        }
      }
    }

    std::vector<Array<const CellId>> cell_list_to_send_by_proc(parallel::size());
    for (size_t i=0; i<parallel::size(); ++i) {
      cell_list_to_send_by_proc[i] = convert_to_array(cell_vector_to_send_by_proc[i]);
    }

    return cell_list_to_send_by_proc;
  }

  Array<int> _buildNbCellToSend()
  {
    Array<int> nb_cell_to_send_by_proc(parallel::size());
    for (size_t i=0; i<parallel::size(); ++i) {
      nb_cell_to_send_by_proc[i] = m_cell_list_to_send_by_proc[i].size();
    }
    return nb_cell_to_send_by_proc;
  }

 public:

  PASTIS_INLINE
  const CellValue<const int>& cellNewOwner() const
  {
    return m_cell_new_owner;
  }

  PASTIS_INLINE
  const NodeValue<const int>& nodeNewOwner() const
  {
    return m_node_new_owner;
  }

  template<typename DataType>
  std::vector<Array<std::remove_const_t<DataType>>>
  exchange(const CellValue<DataType>& cell_value) const
  {
    using MutableDataType = std::remove_const_t<DataType>;
    std::vector<Array<DataType>> cell_value_to_send_by_proc(parallel::size());
    for (size_t i=0; i<parallel::size(); ++i) {
      const Array<const CellId>& cell_list = m_cell_list_to_send_by_proc[i];
      Array<MutableDataType> cell_value_list(cell_list.size());
      parallel_for (cell_list.size(), PASTIS_LAMBDA(const CellId& cell_id) {
          cell_value_list[cell_id] = cell_value[cell_list[cell_id]];
        });
      cell_value_to_send_by_proc[i] = cell_value_list;
    }

    std::vector<Array<MutableDataType>> recv_cell_value_by_proc(parallel::size());
    for (size_t i=0; i<parallel::size(); ++i) {
      recv_cell_value_by_proc[i] = Array<MutableDataType>(m_nb_cell_to_recv_by_proc[i]);
    }

    parallel::exchange(cell_value_to_send_by_proc, recv_cell_value_by_proc);
    return recv_cell_value_by_proc;
  }

  [[deprecated]]
  const CellListToSendByProc& cell_list_to_send_by_proc() const
  {
    return m_cell_list_to_send_by_proc;
  }

  MeshDispatcher(const MeshType& mesh)
      : m_mesh(mesh),
        m_cell_new_owner(_getCellNewOwner()),
        m_node_new_owner(_getNodeNewOwner()),
        m_cell_list_to_send_by_proc(_buildCellListToSend()),
        m_nb_cell_to_send_by_proc(_buildNbCellToSend()),
        m_nb_cell_to_recv_by_proc(parallel::allToAll(m_nb_cell_to_send_by_proc))
  {
    ;
  }

  MeshDispatcher(const MeshDispatcher&) = delete;
  ~MeshDispatcher() = default;
};



template <int Dimension>
void GmshReader::_dispatch()
{
  using ConnectivityType = Connectivity<Dimension>;
  using Rd = TinyVector<Dimension>;
  using MeshType = Mesh<ConnectivityType>;

  if (not m_mesh) {
    ConnectivityDescriptor descriptor;
    std::shared_ptr connectivity = std::make_shared<ConnectivityType>(descriptor);
    NodeValue<Rd> xr;
    m_mesh = std::make_shared<MeshType>(connectivity, xr);
  }
  const MeshType& mesh = static_cast<const MeshType&>(*m_mesh);

  MeshDispatcher<Dimension> dispatcher(mesh);

  std::vector<Array<int>> recv_cell_number_by_proc
      = dispatcher.exchange(mesh.connectivity().cellNumber());

  std::vector<Array<CellType>> recv_cell_type_by_proc
      = dispatcher.exchange(mesh.connectivity().cellType());

  std::vector<Array<int>> recv_cell_new_owner_by_proc
      = dispatcher.exchange(dispatcher.cellNewOwner());

  CellValue<int> number_of_node_per_cell(mesh.connectivity());

  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
  parallel_for(mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
      number_of_node_per_cell[j] = cell_to_node_matrix[j].size();
    });

  std::vector<Array<int>> recv_number_of_node_per_cell_by_proc
      = dispatcher.exchange(number_of_node_per_cell);

  const auto& cell_list_to_send_by_proc = dispatcher.cell_list_to_send_by_proc();
  const NodeValue<const int>& node_number = mesh.connectivity().nodeNumber();


  std::vector<Array<int>> cell_node_number_to_send_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
    std::vector<int> node_number_by_cell_vector;
    for (size_t j=0; j<cell_list_to_send_by_proc[i_rank].size(); ++j) {
      const CellId& cell_id = cell_list_to_send_by_proc[i_rank][j];
      const auto& cell_node_list = cell_to_node_matrix[cell_id];
      for (size_t r=0; r<cell_node_list.size(); ++r) {
        const NodeId& node_id = cell_node_list[r];
        node_number_by_cell_vector.push_back(node_number[node_id]);
      }
    }
    cell_node_number_to_send_by_proc[i_rank] = convert_to_array(node_number_by_cell_vector);
  }

  std::vector<Array<int>> recv_cell_node_number_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
    recv_cell_node_number_by_proc[i_rank]
        = Array<int>(Sum(recv_number_of_node_per_cell_by_proc[i_rank]));
  }

  parallel::exchange(cell_node_number_to_send_by_proc, recv_cell_node_number_by_proc);

  const std::unordered_map<int, int> node_number_id_map
      = [&] () {
          std::unordered_map<int, int> node_number_id_map;
          for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
            int cpt=0;
            for (size_t i=0; i<recv_cell_node_number_by_proc[i_rank].size(); ++i) {
              int node_number = recv_cell_node_number_by_proc[i_rank][i];
              auto [iterator, inserted] = node_number_id_map.insert(std::make_pair(node_number, cpt));
              if (inserted) cpt++;
            }
          }
          return node_number_id_map;
        } ();


  ConnectivityDescriptor new_descriptor;

  new_descriptor.node_number_vector.resize(node_number_id_map.size());
  for (const auto& [number, id] : node_number_id_map) {
    new_descriptor.node_number_vector[id] = number;
  }

  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
    int l=0;
    for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
      std::vector<unsigned int> node_vector;
      for (int k=0; k<recv_number_of_node_per_cell_by_proc[i_rank][i]; ++k) {
        const auto& searched_node_id = node_number_id_map.find(recv_cell_node_number_by_proc[i_rank][l++]);
        Assert(searched_node_id != node_number_id_map.end());
        node_vector.push_back(searched_node_id->second);
      }
      new_descriptor.cell_by_node_vector.emplace_back(node_vector);
    }
  }

  const std::unordered_map<int, int> cell_number_id_map
      = [&] () {
          std::unordered_map<int, int> cell_number_id_map;
          for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
            int cpt=0;
            for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
              int cell_number = recv_cell_number_by_proc[i_rank][i];
              auto [iterator, inserted] = cell_number_id_map.insert(std::make_pair(cell_number, cpt));
              if (inserted) cpt++;
            }
          }
          return cell_number_id_map;
        } ();

  new_descriptor.cell_number_vector.resize(cell_number_id_map.size());
  for (const auto& [number, id] : cell_number_id_map) {
    new_descriptor.cell_number_vector[id] = number;
  }

  new_descriptor.cell_type_vector.resize(cell_number_id_map.size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
      int cell_number = recv_cell_number_by_proc[i_rank][i];
      const auto& searched_cell_id = cell_number_id_map.find(cell_number);
      Assert(searched_cell_id != cell_number_id_map.end());
      new_descriptor.cell_type_vector[searched_cell_id->second] = recv_cell_type_by_proc[i_rank][i];
    }
  }

  std::vector<Array<const NodeId>> send_node_id_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    Array<bool> tag(mesh.numberOfNodes());
    tag.fill(false);
    std::vector<NodeId> node_id_vector;
    for (size_t j=0; j<cell_list_to_send_by_proc[i_rank].size(); ++j) {
      const CellId& cell_id = cell_list_to_send_by_proc[i_rank][j];
      const auto& cell_node_list = cell_to_node_matrix[cell_id];
      for (size_t r=0; r<cell_node_list.size(); ++r) {
        const NodeId& node_id = cell_node_list[r];
        if (not tag[node_id]) {
          node_id_vector.push_back(node_id);
          tag[node_id] = true;
        }
      }
    }
    send_node_id_by_proc[i_rank] = convert_to_array(node_id_vector);
  }

  std::vector<Array<const int>> send_node_number_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    Array<int> send_node_number(send_node_id_by_proc[i_rank].size());
    const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
    parallel_for(send_node_number.size(), PASTIS_LAMBDA(const size_t& j){
        send_node_number[j] = node_number[send_node_id[j]];
      });
    send_node_number_by_proc[i_rank] = send_node_number;
  }

  Array<unsigned int> nb_node_to_send_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    nb_node_to_send_by_proc[i_rank] = send_node_id_by_proc[i_rank].size();
  }
  Array<const unsigned int> nb_node_to_recv_by_proc
      = parallel::allToAll(nb_node_to_send_by_proc);

  std::vector<Array<int>> recv_node_number_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank < parallel::size(); ++i_rank) {
    recv_node_number_by_proc[i_rank] = Array<int>(nb_node_to_recv_by_proc[i_rank]);
  }
  parallel::exchange(send_node_number_by_proc, recv_node_number_by_proc);

  std::vector<Array<const NodeId>> recv_node_id_correspondance_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<nb_node_to_recv_by_proc.size(); ++i_rank) {
    Array<NodeId> node_id_correspondace(nb_node_to_recv_by_proc[i_rank]);
    for (size_t l=0; l<nb_node_to_recv_by_proc[i_rank]; ++l) {
      const int& node_number = recv_node_number_by_proc[i_rank][l];
      const auto& searched_node_id = node_number_id_map.find(node_number);
      Assert(searched_node_id != node_number_id_map.end());
      node_id_correspondace[l] = searched_node_id->second;
    }
    recv_node_id_correspondance_by_proc[i_rank] = node_id_correspondace;
  }

  using ConnectivityType = Connectivity<Dimension>;

  new_descriptor.cell_owner_vector.resize(cell_number_id_map.size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
      int cell_number = recv_cell_number_by_proc[i_rank][i];
      const auto& searched_cell_id = cell_number_id_map.find(cell_number);
      Assert(searched_cell_id != cell_number_id_map.end());
      new_descriptor.cell_owner_vector[searched_cell_id->second] = recv_cell_new_owner_by_proc[i_rank][i];
    }
  }


  const NodeValue<const int>& new_node_owner = dispatcher.nodeNewOwner();
  std::vector<Array<const int>> send_node_owner_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    Array<int> send_node_owner(nb_node_to_send_by_proc[i_rank]);
    const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
    parallel_for(send_node_id.size(), PASTIS_LAMBDA(const size_t& r) {
        const NodeId& node_id = send_node_id[r];
        send_node_owner[r] = new_node_owner[node_id];
      });
    send_node_owner_by_proc[i_rank] = send_node_owner;
  }

  std::vector<Array<int>> recv_node_owner_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    recv_node_owner_by_proc[i_rank] = Array<int>(nb_node_to_recv_by_proc[i_rank]);
  }
  parallel::exchange(send_node_owner_by_proc, recv_node_owner_by_proc);

  new_descriptor.node_owner_vector.resize(new_descriptor.node_number_vector.size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    for (size_t r=0; r<recv_node_owner_by_proc[i_rank].size(); ++r) {
      const NodeId& node_id = recv_node_id_correspondance_by_proc[i_rank][r];
      new_descriptor.node_owner_vector[node_id] = recv_node_owner_by_proc[i_rank][r];
    }
  }

  std::shared_ptr p_connectivity
      = std::make_shared<ConnectivityType>(new_descriptor);

  const NodeValue<const Rd>& xr = mesh.xr();
  std::vector<Array<const Rd>> send_node_coord_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    Array<Rd> send_node_coord(nb_node_to_send_by_proc[i_rank]);
    const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
    parallel_for(send_node_id.size(), PASTIS_LAMBDA(const size_t& r) {
        const NodeId& node_id = send_node_id[r];
        send_node_coord[r] = xr[node_id];
      });
    send_node_coord_by_proc[i_rank] = send_node_coord;
  }

  std::vector<Array<Rd>> recv_node_coord_by_proc(parallel::size());
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    recv_node_coord_by_proc[i_rank] = Array<Rd>(nb_node_to_recv_by_proc[i_rank]);
  }
  parallel::exchange(send_node_coord_by_proc, recv_node_coord_by_proc);

  NodeValue<Rd> new_xr(*p_connectivity);
  for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
    for (size_t r=0; r<recv_node_coord_by_proc[i_rank].size(); ++r) {
      const NodeId& node_id = recv_node_id_correspondance_by_proc[i_rank][r];
      new_xr[node_id] = recv_node_coord_by_proc[i_rank][r];
    }
  }

  m_mesh = std::make_shared<MeshType>(p_connectivity, new_xr);
}


GmshReader::GmshReader(const std::string& filename)
    : m_filename(filename)
{
  if (parallel::rank() == 0) {
    try {
      m_fin.open(m_filename);
      if (not m_fin) {
        perr() << "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;

      pout() << "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 '"+std::to_string(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 '"+std::to_string(fileType)+"'",
                               ErrorHandler::normal);
          }

          int dataSize = this->_getInteger();
          if (dataSize != sizeof(double)) {
            throw ErrorHandler(__FILE__,__LINE__,
                               "Data size not supported '"+std::to_string(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);
            // 	}
            // }

            //       pout() << "- Binary format: ";
            // #ifdef    WORDS_BIGENDIAN
            //       if (not __convertEndian) {
            // 	pout() << "big endian\n";
            //       } else {
            // 	pout() << "little endian\n";
            //       }
            // #else  // WORDS_BIGENDIAN
            //       if (not __convertEndian) {
            // 	pout() << "little endian\n";
            //       } else {
            // 	pout() << "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(const ErrorHandler& e) {
      e.writeErrorMessage();
      std::exit(0);
    }
  }
  pout() << std::flush;
  if (parallel::size() > 1) {
    pout() << "Sequential mesh read! Need to be dispatched\n" << std::flush;

    const int mesh_dimension
        = [&]() {
            int mesh_dimension = -1; // unknown mesh dimension
            if (m_mesh) {
              mesh_dimension = m_mesh->meshDimension();
            }

            Array<int> dimensions = parallel::allGather(mesh_dimension);
            std::set<int> dimension_set;
            for (size_t i=0; i<dimensions.size(); ++i) {
              const int i_dimension = dimensions[i];
              if (i_dimension != -1) {
                dimension_set.insert(i_dimension);
              }
            }
            if (dimension_set.size() != 1) {
              std::cerr << "error dimensions of read mesh parts differ!\n";
              std::exit(1);
            }

            return *begin(dimension_set);
          }();
    switch (mesh_dimension) {
      case 1: {
        this->_dispatch<1>();
        break;
      }
      case 2: {
        this->_dispatch<2>();
        break;
      }
      case 3: {
        this->_dispatch<3>();
        break;
      }
      default: {
        perr() << "unexpected dimension " << mesh_dimension << '\n';
        std::exit(1);
      }
    }
  }
}

void GmshReader::__readVertices()
{
  const int numberOfVerices = this->_getInteger();
  pout() << "- 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 = Array<R3>(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();
//   pout() << "- 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 '"+std::to_string(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();
  pout() << "- Number Of Elements: " << numberOfElements << '\n';
  if (numberOfElements<0) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "reading file '"+m_filename
                       +"': number of elements is negative",
                       ErrorHandler::normal);
  }

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

  if (not __binary) {
    for (int i=0; i<numberOfElements; ++i) {
      __elementNumber[i] = this->_getInteger();
      const int elementType = this->_getInteger()-1;

      if ((elementType < 0) or (elementType > 14)) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "reading file '"+m_filename
                           +"': unknown element type '"+std::to_string(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 '"+std::to_string(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()) {
      perr() << "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::__compute2DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
{
  // const auto& cell_to_node_matrix
  //     = this->_getMatrix(ItemType::cell, ItemType::node);
  using Face = ConnectivityFace<2>;

  // In 2D faces are simply define
  using CellFaceId = std::pair<CellId, unsigned short>;
  std::map<Face, std::vector<CellFaceId>> face_cells_map;
  for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
    const auto& cell_nodes = descriptor.cell_by_node_vector[j];
    for (unsigned short r=0; r<cell_nodes.size(); ++r) {
      NodeId node0_id = cell_nodes[r];
      NodeId node1_id = cell_nodes[(r+1)%cell_nodes.size()];
      if (node1_id<node0_id) {
        std::swap(node0_id, node1_id);
      }
      face_cells_map[Face({node0_id, node1_id})].push_back(std::make_pair(j, r));
    }
  }

  {
    descriptor.face_to_node_vector.resize(face_cells_map.size());
    int l=0;
    for (const auto& face_info : face_cells_map) {
      const Face& face = face_info.first;
      descriptor.face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
      ++l;
    }
  }

  {
    descriptor.face_to_cell_vector.resize(face_cells_map.size());
    int l=0;
    for (const auto& face_cells_vector : face_cells_map) {
      const auto& [face, cell_info_vector] = face_cells_vector;
      for (const auto& cell_info : cell_info_vector) {
        descriptor.face_to_cell_vector[l].push_back(cell_info.first);
      }
      ++l;
    }
  }

  {
    descriptor.face_number_vector.resize(face_cells_map.size());
    for (size_t l=0; l<face_cells_map.size(); ++l) {
      descriptor.face_number_vector[l] = l;
    }
    perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
  }
}

void
GmshReader::__compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
{
  using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
  using Face = ConnectivityFace<3>;

  // const auto& cell_to_node_matrix
  //     = this->_getMatrix(ItemType::cell, ItemType::node);

  Array<unsigned short> cell_nb_faces(descriptor.cell_by_node_vector.size());
  std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
  for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
    const auto& cell_nodes = descriptor.cell_by_node_vector[j];

    switch (descriptor.cell_type_vector[j]) {
      case CellType::Tetrahedron: {
        cell_nb_faces[j] = 4;
        // face 0
        Face f0({cell_nodes[1],
                 cell_nodes[2],
                 cell_nodes[3]});
        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));

        // face 1
        Face f1({cell_nodes[0],
                 cell_nodes[3],
                 cell_nodes[2]});
        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));

        // face 2
        Face f2({cell_nodes[0],
                 cell_nodes[1],
                 cell_nodes[3]});
        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));

        // face 3
        Face f3({cell_nodes[0],
                 cell_nodes[2],
                 cell_nodes[1]});
        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
        break;
      }
      case CellType::Hexahedron: {
        // face 0
        Face f0({cell_nodes[3],
                 cell_nodes[2],
                 cell_nodes[1],
                 cell_nodes[0]});
        face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));

        // face 1
        Face f1({cell_nodes[4],
                 cell_nodes[5],
                 cell_nodes[6],
                 cell_nodes[7]});
        face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));

        // face 2
        Face f2({cell_nodes[0],
                 cell_nodes[4],
                 cell_nodes[7],
                 cell_nodes[3]});
        face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));

        // face 3
        Face f3({cell_nodes[1],
                 cell_nodes[2],
                 cell_nodes[6],
                 cell_nodes[5]});
        face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));

        // face 4
        Face f4({cell_nodes[0],
                 cell_nodes[1],
                 cell_nodes[5],
                 cell_nodes[4]});
        face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));

        // face 5
        Face f5({cell_nodes[3],
                 cell_nodes[7],
                 cell_nodes[6],
                 cell_nodes[2]});
        face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));

        cell_nb_faces[j] = 6;
        break;
      }
      default: {
        perr() << "unexpected cell type!\n";
        std::exit(0);
      }
    }
  }

  {
    descriptor.cell_to_face_vector.resize(descriptor.cell_by_node_vector.size());
    for (CellId j=0; j<descriptor.cell_to_face_vector.size(); ++j) {
      descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
    }
    FaceId l=0;
    for (const auto& face_cells_vector : face_cells_map) {
      const auto& cells_vector = face_cells_vector.second;
      for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
        const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
        descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
      }
      ++l;
    }
  }

  {
    descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_by_node_vector.size());
    for (CellId j=0; j<descriptor.cell_face_is_reversed_vector.size(); ++j) {
      descriptor.cell_face_is_reversed_vector[j].resize(cell_nb_faces[j]);
    }
    for (const auto& face_cells_vector : face_cells_map) {
      const auto& cells_vector = face_cells_vector.second;
      for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
        const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
        descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
      }
    }
  }

  {
    descriptor.face_to_node_vector.resize(face_cells_map.size());
    int l=0;
    for (const auto& face_info : face_cells_map) {
      const Face& face = face_info.first;
      descriptor.face_to_node_vector[l] = face.nodeIdList();
      ++l;
    }
  }

  {
    descriptor.face_number_vector.resize(face_cells_map.size());
    for (size_t l=0; l<face_cells_map.size(); ++l) {
      descriptor.face_number_vector[l] = l;
    }
    perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
  }
}


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

        switch (i) {
          // Supported entities
          case 0: { // edges
            __edges = Array<Edge>(elementNumber[i]);
            __edges_ref.resize(elementNumber[i]);
            __edges_number.resize(elementNumber[i]);
            break;
          }
          case 1: { // triangles
            __triangles = Array<Triangle>(elementNumber[i]);
            __triangles_ref.resize(elementNumber[i]);
            __triangles_number.resize(elementNumber[i]);
            break;
          }
          case  2: { // quadrangles
            __quadrangles = Array<Quadrangle>(elementNumber[i]);
            __quadrangles_ref.resize(elementNumber[i]);
            __quadrangles_number.resize(elementNumber[i]);
            break;
          }
          case 3: { // tetrahedra
            __tetrahedra = Array<Tetrahedron>(elementNumber[i]);
            __tetrahedra_ref.resize(elementNumber[i]);
            __tetrahedra_number.resize(elementNumber[i]);
            break;
          }
          case  4: {// hexahedra
            __hexahedra = Array<Hexahedron>(elementNumber[i]);
            __hexahedra_ref.resize(elementNumber[i]);
            __hexahedra_number.resize(elementNumber[i]);
            break;
          }
          case 14: {// point
            __points = Array<Point>(elementNumber[i]);
            __points_ref.resize(elementNumber[i]);
            __points_number.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);
          }
        }
      }
      pout() << '\n';
    }
  }

  pout() << "- 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);
  }

#warning should use an unordered_map
  // 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 "+std::to_string(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 "+std::to_string(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __triangles[triangleNumber]
            = Triangle(a, b, c);
        __triangles_ref[triangleNumber] = __references[i];
        __triangles_number[triangleNumber] = __elementNumber[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 "+std::to_string(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __quadrangles[quadrilateralNumber] = Quadrangle(a,b,c,d);
        __quadrangles_ref[quadrilateralNumber] = __references[i];
        __quadrangles_number[quadrilateralNumber] = __elementNumber[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 "+std::to_string(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __tetrahedra[tetrahedronNumber] = Tetrahedron(a, b, c, d);
        __tetrahedra_ref[tetrahedronNumber] = __references[i];
        __tetrahedra_number[tetrahedronNumber] = __elementNumber[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 "+std::to_string(i)
                             +" [bad vertices definition]",
                             ErrorHandler::normal);
        }
        __hexahedra[hexahedronNumber]
            = Hexahedron(a, b, c, d,
                         e, f, g, h);
        __hexahedra_ref[hexahedronNumber] = __references[i];
        __hexahedra_number[hexahedronNumber] = __elementNumber[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];
        __points_number[point_number] = __elementNumber[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;

  pout() << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
  pout() << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
  pout() << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
  pout() << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
  if ((dimension3_mask, elementNumber)>0) {
    const size_t nb_cells = (dimension3_mask, elementNumber);

    ConnectivityDescriptor descriptor;

    descriptor.node_number_vector = __verticesNumbers;
    descriptor.cell_type_vector.resize(nb_cells);
    descriptor.cell_number_vector.resize(nb_cells);
    descriptor.cell_by_node_vector.resize(nb_cells);

    const size_t nb_tetrahedra = __tetrahedra.size();
    for (size_t j=0; j<nb_tetrahedra; ++j) {
      descriptor.cell_by_node_vector[j].resize(4);
      for (int r=0; r<4; ++r) {
        descriptor.cell_by_node_vector[j][r] = __tetrahedra[j][r];
      }
      descriptor.cell_type_vector[j] = CellType::Tetrahedron;
      descriptor.cell_number_vector[j] = __tetrahedra_number[j];
    }
    const size_t nb_hexahedra = __hexahedra.size();
    for (size_t j=0; j<nb_hexahedra; ++j) {
      const size_t jh = nb_tetrahedra+j;
      descriptor.cell_by_node_vector[jh].resize(8);
      for (int r=0; r<8; ++r) {
        descriptor.cell_by_node_vector[jh][r] = __hexahedra[j][r];
      }
      descriptor.cell_type_vector[jh] = CellType::Hexahedron;
      descriptor.cell_number_vector[jh] = __hexahedra_number[j];
    }

    __compute3DCellFaceAndFaceNodeConnectivities(descriptor);

    descriptor.cell_owner_vector.resize(nb_cells);
    std::fill(descriptor.cell_owner_vector.begin(),
              descriptor.cell_owner_vector.end(),
              parallel::rank());

    descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
    std::fill(descriptor.face_owner_vector.begin(),
              descriptor.face_owner_vector.end(),
              parallel::rank());

    descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
    std::fill(descriptor.node_owner_vector.begin(),
              descriptor.node_owner_vector.end(),
              parallel::rank());

    std::shared_ptr p_connectivity = std::make_shared<Connectivity3D>(descriptor);
    Connectivity3D& connectivity = *p_connectivity;

    using Face = ConnectivityFace<3>;
    const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
        = [&]  {
            std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
            const auto& face_node_list = connectivity.faceToNodeMatrix();
            for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
              const auto& node_list = face_node_list[l];
              std::vector<unsigned int> node_vector(node_list.size());
              for (size_t r=0; r<node_list.size(); ++r) {
                node_vector[r] = node_list[r];
              }
              face_to_id_map[Face(node_vector)] = l;
            }
            return face_to_id_map;
          } ();

    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
    for (unsigned int f=0; f<__triangles.size(); ++f) {
      const unsigned int face_number
          = [&]{
              auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]}));
              if (i == face_to_id_map.end()) {
                std::cerr << "face not found!\n";
                std::terminate();
              }
              return i->second;
            }();

      const unsigned int& ref = __triangles_ref[f];
      ref_faces_map[ref].push_back(face_number);
    }
    for (unsigned int f=0; f<__quadrangles.size(); ++f) {
      const unsigned int face_number
          = [&]{
              auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
                                                 __quadrangles[f][2], __quadrangles[f][3]}));
              if (i == face_to_id_map.end()) {
                std::cerr << "face not found!\n";
                std::terminate();
              }
              return i->second;
            }();

      const unsigned int& ref = __quadrangles_ref[f];
      ref_faces_map[ref].push_back(face_number);
    }
    for (const auto& ref_face_list : ref_faces_map) {
      Array<FaceId> 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));
    }

    using MeshType = Mesh<Connectivity3D>;
    using Rd = TinyVector<3, double>;

    NodeValue<Rd> xr(connectivity);
    for (NodeId i=0; i<__vertices.size(); ++i) {
      xr[i][0] = __vertices[i][0];
      xr[i][1] = __vertices[i][1];
      xr[i][2] = __vertices[i][2];
    }
    m_mesh = std::make_shared<MeshType>(p_connectivity, xr);

  } else if ((dimension2_mask, elementNumber)>0) {
    const size_t nb_cells = (dimension2_mask, elementNumber);

    ConnectivityDescriptor descriptor;

    descriptor.node_number_vector = __verticesNumbers;
    descriptor.cell_type_vector.resize(nb_cells);
    descriptor.cell_number_vector.resize(nb_cells);
    descriptor.cell_by_node_vector.resize(nb_cells);

    const size_t nb_triangles = __triangles.size();
    for (size_t j=0; j<nb_triangles; ++j) {
      descriptor.cell_by_node_vector[j].resize(3);
      for (int r=0; r<3; ++r) {
        descriptor.cell_by_node_vector[j][r] = __triangles[j][r];
      }
      descriptor.cell_type_vector[j] = CellType::Triangle;
      descriptor.cell_number_vector[j] = __triangles_number[j];
    }

    const size_t nb_quadrangles = __quadrangles.size();
    for (size_t j=0; j<nb_quadrangles; ++j) {
      const size_t jq = j+nb_triangles;
      descriptor.cell_by_node_vector[jq].resize(4);
      for (int r=0; r<4; ++r) {
        descriptor.cell_by_node_vector[jq][r] = __quadrangles[j][r];
      }
      descriptor.cell_type_vector[jq] = CellType::Quadrangle;
      descriptor.cell_number_vector[jq] = __quadrangles_number[j];
    }

    this->__compute2DCellFaceAndFaceNodeConnectivities(descriptor);

    descriptor.cell_owner_vector.resize(nb_cells);
    std::fill(descriptor.cell_owner_vector.begin(),
              descriptor.cell_owner_vector.end(),
              parallel::rank());

    descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
    std::fill(descriptor.face_owner_vector.begin(),
              descriptor.face_owner_vector.end(),
              parallel::rank());

    descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
    std::fill(descriptor.node_owner_vector.begin(),
              descriptor.node_owner_vector.end(),
              parallel::rank());

    std::shared_ptr p_connectivity = std::make_shared<Connectivity2D>(descriptor);
    Connectivity2D& connectivity = *p_connectivity;

    using Face = ConnectivityFace<2>;
    const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
        = [&]  {
            std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
            const auto& face_node_list = connectivity.faceToNodeMatrix();
            for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
              const auto& node_list = face_node_list[l];
              face_to_id_map[Face({node_list[0], node_list[1]})] = l;
            }
            return face_to_id_map;
          } ();

    std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
    for (unsigned int e=0; e<__edges.size(); ++e) {
      const unsigned int edge_number
          = [&]{
              auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}));
              if (i == face_to_id_map.end()) {
                std::cerr << "face " << __edges[e][0] << " not found!\n";
                std::terminate();
              }
              return i->second;
            }();
      const unsigned int& ref = __edges_ref[e];
      ref_faces_map[ref].push_back(edge_number);
    }

    for (const auto& ref_face_list : ref_faces_map) {
      Array<FaceId> 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.size(); ++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) {
      Array<NodeId> 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));
    }

    using MeshType = Mesh<Connectivity2D>;
    using Rd = TinyVector<2, double>;

    NodeValue<Rd> xr(connectivity);
    for (NodeId i=0; i<__vertices.size(); ++i) {
      xr[i][0] = __vertices[i][0];
      xr[i][1] = __vertices[i][1];
    }
    m_mesh = std::make_shared<MeshType>(p_connectivity, xr);

  } else if ((dimension1_mask, elementNumber)>0) {
    const size_t nb_cells = (dimension1_mask, elementNumber);

    ConnectivityDescriptor descriptor;

    descriptor.node_number_vector = __verticesNumbers;
    descriptor.cell_type_vector.resize(nb_cells);
    descriptor.cell_number_vector.resize(nb_cells);
    descriptor.cell_by_node_vector.resize(nb_cells);

    for (size_t j=0; j<nb_cells; ++j) {
      descriptor.cell_by_node_vector[j].resize(2);
      for (int r=0; r<2; ++r) {
        descriptor.cell_by_node_vector[j][r] = __edges[j][r];
      }
      descriptor.cell_type_vector[j] = CellType::Line;
      descriptor.cell_number_vector[j] =  __edges_number[j];
    }

    descriptor.cell_owner_vector.resize(nb_cells);
    std::fill(descriptor.cell_owner_vector.begin(),
              descriptor.cell_owner_vector.end(),
              parallel::rank());

    descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
    std::fill(descriptor.node_owner_vector.begin(),
              descriptor.node_owner_vector.end(),
              parallel::rank());


    std::shared_ptr p_connectivity = std::make_shared<Connectivity1D>(descriptor);
    Connectivity1D& connectivity = *p_connectivity;

    std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
    for (unsigned int r=0; r<__points.size(); ++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) {
      Array<NodeId> 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));
    }

    using MeshType = Mesh<Connectivity1D>;
    using Rd = TinyVector<1, double>;

    NodeValue<Rd> xr(connectivity);
    for (NodeId i=0; i<__vertices.size(); ++i) {
      xr[i][0] = __vertices[i][0];
    }

    m_mesh = std::make_shared<MeshType>(p_connectivity, xr);

  } else {
    perr() << "*** 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()
{
  pout() << "- 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);
      }
    }
  }
}