msh2geo.cc

Go to the documentation of this file.

//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
// 
// =========================================================================
// the msh2geo unix command
// author: Pierre.Saramito@imag.fr
 
namespace rheolef {
} // namespace rheolef
 
#include "rheolef/compiler.h" // after index_set to avoid boost
#include "scatch.icc"
#include "index_set_header.icc"
#include "index_set_body.icc"
#include <array>
#include <cstring>
 
namespace rheolef {
 
template<class T>
struct point_basic: std::array<T,3> {
  typedef std::array<T,3> base;
  point_basic(T x0=T(), T x1=T(), T x2=T()) 
   : base() { 
    base::operator[](0) = x0;
    base::operator[](1) = x1;
    base::operator[](2) = x2;
  }
};
using point = point_basic<double>;
 
namespace edge {
#include "edge.icc"
} // namespace edge
 
namespace triangle {
#include "triangle.icc"
} // namespace triangle
 
namespace quadrangle {
#include "quadrangle.icc"
} // namespace quadrangle
 
namespace tetrahedron {
#include "tetrahedron.icc"
} // namespace tetrahedron
 
namespace prism {
#include "prism.icc"
} // namespace prism
 
namespace hexahedron {
#include "hexahedron.icc"
} // namespace hexahedron
 
} // namespace rheolef
 
#include "reference_element_aux.icc"
 
namespace rheolef {
 
// TODO: move somewhere in reference_element_xxx.cc
const char
_reference_element_name [reference_element__max_variant] = {
        'p',
        'e',
        't',
        'q',
        'T',
        'P',
        'H'
};
} // namespace rheolef
 
namespace rheolef {
 
typedef int orientation_type;
typedef int shift_type;
 
// side (edge & face): element with reduced node list when order > 1, for sides (edges or faces)
struct geo_element_side: std::array<size_t,4> {
  void setname     (char   name)     { _variant = reference_element_variant(name); }
  void setindex    (size_t index)    { _index   = index; }
  char   name()     const { return _reference_element_name [_variant]; }
  size_t index()    const { return _index; }
  size_t variant()  const { return _variant; }
  size_t dimension()const { return reference_element_dimension_by_variant(variant()); }
  size_t n_vertex() const { return reference_element_n_vertex (variant()); }
  size_t size()     const { return n_vertex(); }
  geo_element_side() 
   : array<size_t,4>(), _variant(-1), _index(-1)
   {}
protected:
  size_t _variant;
  size_t _index;
};
// side index with orient and shift for geo_element
struct geo_element_indirect {
  typedef int orientation_type;
  typedef int shift_type;
  void setindex       (size_t index)  { _index  = index; }
  void setorientation (size_t orient) { _orient = orient; }
  void setshift       (size_t shift)  { _shift  = shift; }
  size_t index()       const { return _index; }
  int    orientation() const { return _orient; }
  size_t shift()       const { return _shift; }
  geo_element_indirect()
   : _index(-1), _orient(1), _shift(0) {}
protected:
  size_t _index;
  int    _orient;
  int    _shift;
};
// element with full node list when order > 1
struct geo_element: vector<size_t> {
  void setname     (char   name)     { _name     = name; }
  void setorder    (size_t order)    { _order    = order; }
  void setindex    (size_t index)    { _index    = index; }
  void setgmshtype (size_t gmshtype) { _gmshtype = gmshtype; }
  geo_element_indirect& edge_indirect (size_t i) { return _edge[i]; }
  geo_element_indirect& face_indirect (size_t i) { return _face[i]; }
  char   name()     const { return _name; }
  size_t order()    const { return _order; }
  size_t index()    const { return _index; }
  size_t gmshtype() const { return _gmshtype; }
  size_t variant()  const { return reference_element_variant  (name()); }
  size_t dimension()const { return reference_element_dimension_by_name(name()); }
  size_t n_vertex() const { return reference_element_n_vertex (variant()); }
  const geo_element_indirect& edge_indirect (size_t i) const { return _edge[i]; }
  const geo_element_indirect& face_indirect (size_t i) const { return _face[i]; }
  size_t n_edge() const { return _reference_element_n_edge_by_variant [variant()]; }
  size_t edge (size_t i) const { return _edge[i].index(); }
  size_t edge_local_vertex(size_t iedg, size_t edg_iv) const;
  size_t face (size_t i) const { return _face[i].index(); }
  size_t n_face() const { return _reference_element_n_face_by_variant [variant()]; }
  size_t face_n_vertex(size_t loc_ifac) const;
  size_t face_local_vertex(size_t iedg, size_t edg_iv) const;
  geo_element() 
   : vector<size_t>(), _name('z'), _order(-1), _index(-1), _gmshtype(-1), _edge(), _face()
   {}
protected:
  char   _name; // TODO: store variant instead of name
  size_t _order;
  size_t _index;
  size_t _gmshtype;
  array<geo_element_indirect,12> _edge;
  array<geo_element_indirect,6>  _face;
};
size_t
geo_element::edge_local_vertex (size_t iedg, size_t edg_iv) const
{
  switch (variant()) {
    case reference_element__e: return vector<size_t>::operator[] (edg_iv);
    case reference_element__t: return    triangle::edge [iedg][edg_iv%2];
    case reference_element__q: return  quadrangle::edge [iedg][edg_iv%2];
    case reference_element__T: return tetrahedron::edge [iedg][edg_iv%2];
    case reference_element__P: return       prism::edge [iedg][edg_iv%2];
    case reference_element__H: return  hexahedron::edge [iedg][edg_iv%2];
    default: error_macro ("invalid variant"); return 0;
  }
}
size_t
geo_element::face_local_vertex (size_t loc_ifac, size_t fac_iv) const
{
  switch (variant()) {
    case reference_element__t:
    case reference_element__q: return vector<size_t>::operator[] (fac_iv);
    case reference_element__T: return tetrahedron::face [loc_ifac][fac_iv%3];
    case reference_element__P: return       prism::face [loc_ifac][fac_iv%4];
    case reference_element__H: return  hexahedron::face [loc_ifac][fac_iv%4];
    default: error_macro ("invalid variant"); return 0;
  }
}
size_t
geo_element::face_n_vertex(size_t loc_ifac) const
{
  switch (variant()) {
    case reference_element__t: return 3;
    case reference_element__q: return 4;
    case reference_element__T: return 3;
    case reference_element__H: return 4;
    case reference_element__P: return (loc_ifac < 2) ? 3 : 4;
    default: error_macro ("invalid variant"); return 0;
  }
}
 
} // namespace rheolef
 
#include "msh2geo_defs.icc"
#include "msh2geo_node_renum.icc"
#include "geo_element_aux.icc"
 
using namespace std;
using namespace rheolef;
// necessite de numeroter les aretes 
// puis de le stoker dans ts les triangles
// le num arete et son orient
// => doit etre appeler apres la num des aretes
void
numbering_Pk_dis_idof (
  size_t                                                order,
  const array<size_t,reference_element__max_variant>&   size_by_variant,
  const array<size_t,reference_element__max_variant+1>& first_by_variant,
  const array<size_t,reference_element__max_variant>&   loc_ndof_by_variant,
  const geo_element&                                    K, 
  vector<size_t>&                                       dis_idof1)
{
  dis_idof1.resize (reference_element_n_node (K.variant(), order));
  std::fill (dis_idof1.begin(), dis_idof1.end(), std::numeric_limits<size_t>::max());
 
  for (size_t subgeo_variant = 0, variant = K.variant(); subgeo_variant <= variant; subgeo_variant++) {
    size_t subgeo_dim   = reference_element_dimension_by_variant (subgeo_variant);
    size_t loc_sub_ndof = loc_ndof_by_variant [subgeo_variant];
    for (size_t first_loc_idof = reference_element_first_inod_by_variant (variant, order, subgeo_variant),
                 last_loc_idof = reference_element_last_inod_by_variant (variant, order, subgeo_variant), 
                loc_idof = first_loc_idof;
                loc_idof <  last_loc_idof; loc_idof++) {
      // 1) local loc_igev on subgeo
      size_t loc_igev   = (loc_idof - first_loc_idof) / loc_sub_ndof;
      size_t loc_igev_j = (loc_idof - first_loc_idof) % loc_sub_ndof;
      // 2) then compute ige; computation depends upon the subgeo dimension:
      size_t ige;
      switch (subgeo_dim) {
        case 0: {
          // convert node numbering to vertex numbering, for geo order > 1
          size_t loc_inod = loc_idof; // only one dof per vertex
          ige = K [loc_inod];
          break;
        }
        case 1: {
          loc_igev_j = geo_element_fix_edge_indirect (K, loc_igev, loc_igev_j, order);
          size_t loc_ige = loc_igev;
          ige = K.edge_indirect(loc_ige).index();
          break;
        }
        case 2: {
          size_t loc_ige = loc_igev;
          if (subgeo_variant == reference_element__t) {
             loc_igev_j = geo_element_fix_triangle_indirect (K, loc_igev, loc_igev_j, order);
          } else {
             size_t loc_ntri = (K.variant() == reference_element__P) ? 2 : 0;
             loc_ige += loc_ntri;
             loc_igev_j = geo_element_fix_quadrangle_indirect (K, loc_ige, loc_igev_j, order);
          }
          ige = K.face(loc_ige);
          break;
        }
        case 3: {
          ige = K.index();
          break;
        }
        default: {
          error_macro ("unexpected subgeo_dim="<<subgeo_dim);
        }
      }
      // 3) then compute igev, by variant
      size_t igev = ige;
      for (size_t prev_subgeo_variant = reference_element_first_variant_by_dimension(subgeo_dim);
                     prev_subgeo_variant < subgeo_variant; 
                     prev_subgeo_variant++) {
        size_t nprev = size_by_variant [prev_subgeo_variant];
        assert_macro (ige >= nprev, "invalid index");
        igev -= nprev;
      } 
      size_t dis_igev = igev;
      // 4) finally compute dis_idof
      size_t dis_idof = 0;
      for (size_t prev_subgeo_variant = 0;
                     prev_subgeo_variant < subgeo_variant; 
                     prev_subgeo_variant++) {
        dis_idof += size_by_variant     [prev_subgeo_variant]
                   *loc_ndof_by_variant [prev_subgeo_variant];
      }
      dis_idof +=   dis_igev
                   *loc_ndof_by_variant [subgeo_variant]
                   + loc_igev_j;
      assert_macro (loc_idof < dis_idof1.size(), "msh2geo: invalid index");
      dis_idof1 [loc_idof] = dis_idof;
    }
  }
}
void
put_domain_noupgrade (
  ostream&                         out,
  size_t                           dim,
  size_t                           dom_dim,
  const map<size_t,list<size_t> >& domain_map,
        map<size_t, string>&       phys,
  const vector<geo_element>&       element)
{
  if (dim == dom_dim && domain_map.size() == 1) return;
  for (map<size_t,list<size_t> >::const_iterator
          first = domain_map.begin(),
          last  = domain_map.end(); first != last; first++) {
    size_t           dom_idx = (*first).first;
    const list<size_t>& dom     = (*first).second;
    string dom_name = phys[dom_idx];
    if (dom_name == "") dom_name = "unnamed" + std::to_string(dom_idx);
    out << "domain" << endl
        << dom_name << endl
        << "1 " << dom_dim << " " << dom.size() << endl;
    for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
                variant <  reference_element_last_variant_by_dimension(dom_dim); variant++) {
      for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
        size_t ie = *fe;
        if (element[ie].variant() != variant) continue;
        out << element[ie].name() << "\t";
        if (element[ie].order() > 1) {
          out << "p" << element[ie].order() << " ";
        }
        for (size_t j = 0; j < element[ie].size(); j++) {
          out << element[ie][j] << " ";
        }
        out << endl;
      }
    }
    out << endl;
  }
}
void
put_domain_upgrade (
  ostream&                           out,
  size_t                             dim,
  size_t                             dom_dim,
  const map<size_t,list<size_t> >&   domain_map,
        map<size_t, string>&         phys,
  const vector<geo_element>&         element,
  const vector<geo_element_side>&    edge,
  const vector<index_set>&           edge_ball,
  const vector<geo_element_side>&    face,
  const vector<index_set>&           face_ball)
{
  if (dim == dom_dim && domain_map.size() == 1) return;
  for (map<size_t,list<size_t> >::const_iterator
          first = domain_map.begin(),
          last  = domain_map.end(); first != last; first++) {
    size_t           dom_idx = (*first).first;
    const list<size_t>& dom     = (*first).second;
    string dom_name = phys[dom_idx];
    if (dom_name == "") dom_name = "unnamed" + std::to_string(dom_idx);
    out << "domain" << endl
        << dom_name << endl
        << "2 " << dom_dim << " " << dom.size() << endl;
    for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
                variant <  reference_element_last_variant_by_dimension(dom_dim); variant++) {
      for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
        size_t ie = *fe;
        if (element[ie].variant() != variant) continue;
        size_t isid = 0;
        orientation_type orient = 1;
        if (element[ie].name() == 'p') {
          isid   = element[ie][0];
          orient = 1;
        } else if (element[ie].name() == 'e') {
          size_t inod0 = element[ie][0];
          size_t inod1 = element[ie][1];
          index_set iedge_set = edge_ball[inod0];
          iedge_set.inplace_intersection (edge_ball[inod1]);
          check_macro (iedge_set.size() == 1,
            "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
          isid = *(iedge_set.begin());
          orient = (inod0 == edge[isid][0]) ? 1 : -1;
        } else if (element[ie].dimension() == 2) {
          size_t inod0 = element[ie][0];
          index_set iface_set = face_ball[inod0];
          for (size_t j = 1; j < element[ie].n_vertex(); ++j) {
            size_t inodj = element[ie][j];
            iface_set.inplace_intersection (face_ball[inodj]);
          }
          check_macro (iface_set.size() == 1,
            "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
          isid = *(iface_set.begin());
          shift_type shift;
          if (element[ie].name() == 't') {
            geo_element_get_orientation_and_shift (
              face[isid],
              element[ie][0],
              element[ie][1],
              element[ie][2],
              orient, shift);
          } else {
            geo_element_get_orientation_and_shift (
              face[isid],
              element[ie][0],
              element[ie][1],
              element[ie][2],
              element[ie][3],
              orient, shift);
          }
        } else { // 3d domain
          isid   = element[ie].index();
          orient = 1;
        }
        out << orient*int(isid) << endl;
      }
    }
    out << endl;
  }
}
void msh2geo (istream& in, ostream& out, string sys_coord_name, bool do_upgrade)
{
  // ----------------------------------
  // 1. input gmsh
  // ----------------------------------
  // 1.0. preambule
  check_macro (scatch(in,"$MeshFormat",true),
        "msh2geo: input stream does not contains a gmsh mesh file ($MeshFormat not found).");
  double gmsh_fmt_version;
  size_t file_type, float_data_size;
  in >> gmsh_fmt_version >> file_type >> float_data_size;
  check_macro (gmsh_fmt_version <= 2.2,
    "gmsh format version " << gmsh_fmt_version << " founded ; expect version 2.2 (HINT: use gmsh -format msh2)");
  check_macro (file_type == 0, "msh2geo: unsupported gmsh non-ascii format");
  check_macro (scatch(in,"$EndMeshFormat",true), "msh2geo: gmsh input error: $EndMeshFormat not found.");
  //
  // 1.1 optional domain names
  //
  bool full_match = true;
  bool partial_match = !full_match;
  check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
  size_t nphys;
  map<size_t, string> phys;
  string label;
  in >> label;
  size_t n_names = 0;
  if (label == "PhysicalNames") {
    in  >> nphys;
    for (size_t i = 0; i < nphys; i++) {
      string name;
      size_t dom_dim, dom_idx;
      in  >> dom_dim >> dom_idx;
      // get name:
      char c; 
      in >> std::ws >>  c;
      if (c != '"') name.push_back(c);
      do {
        in.get(c);
        name.push_back(c);
      } while (c != '"' && c != '\n');
      // strip '"' in name
      size_t start = 0, end = name.length();
      if (name[start] == '"') start++;
      if (name[end-1] == '"') end--;
      name = name.substr(start,end-start);
      // rename spaces and tabs
      for (size_t i = 0; i < name.size(); i++) {
        if (name[i] == ' ' || name[i] == '\t') name[i] = '_';
      }
      phys[dom_idx] = name.substr(start,end-start);
    }
    check_macro (scatch(in,"$EndPhysicalNames",true), "msh2geo: gmsh input error ($EndPhysicalNames not found).");
    check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
    in >> label;
  }
  //
  // 1.2. nodes
  //
  if (label != "Nodes") {
    check_macro (scatch(in,"$Nodes",true), "msh2geo: $Nodes not found in gmsh file");
  }
  size_t nnode;
  in  >> nnode;
  vector<point> node (nnode);
  double infty = numeric_limits<double>::max();
  point xmin ( infty,  infty,  infty);
  point xmax (-infty, -infty, -infty);
  for (size_t k = 0; k < node.size(); k++) {
    size_t dummy;
    in  >> dummy;
    for (size_t i = 0; i < 3; ++i) {
      in >> node[k][i];
    }
    for (size_t j = 0 ; j < 3; j++) {
      xmin[j] = min(node[k][j], xmin[j]);
      xmax[j] = max(node[k][j], xmax[j]);
    }
  }
  //
  // dimension is deduced from bounding box
  //
  size_t dim = 3;
  if (xmax[2] == xmin[2]) {
    dim = 2;
    if (xmax[1] == xmin[1]) {
      dim = 1;
      if (xmax[0] == xmin[0]) dim = 0;
    }
  }
  //
  // 1.3. elements
  //
  check_macro (scatch(in,"$Elements",true), "msh2geo: $Elements not found in gmsh file");
  size_t n_element;
  in  >> n_element;
  vector<geo_element> element (n_element);
  vector<size_t>  node_subgeo_variant (node.size(), std::numeric_limits<size_t>::max());
  array<map<size_t, list<size_t> >,4>  domain_map; // domain_map[dim][idom][ie]
  size_t map_dim = 0;
  size_t order = 0;
  for (size_t i = 0; i < element.size(); i++) {
    size_t id, dummy, gmshtype;
    in  >> id >> gmshtype;
    size_t n_tag_gmsh;
    in >> n_tag_gmsh;
    size_t domain_idx = 0;
    for (size_t j = 0 ; j < n_tag_gmsh; j++) {
        // the first tag is the physical domain index
        // the second tag is the object index, defined for all elements
        // the third is zero (in all examples)
        size_t tag_dummy;
        in >> tag_dummy;
        if (j == 0) {
          domain_idx = tag_dummy;
        }
    }
    check_macro (gmshtype < gmshtype_max,
        "msh2geo: element #" << id << ": unexpected gmsh type '" << gmshtype << "'");
    check_macro (gmsh_table[gmshtype].supported,
        "msh2geo: element #" << id << ": unsupported gmsh type '" << gmshtype << "'");
 
    element[i].setname (gmsh_table[gmshtype].name);
    element[i].setorder(gmsh_table[gmshtype].order);
    element[i].setgmshtype(gmshtype);
    size_t nv    = gmsh_table[gmshtype].nv;
    size_t nn    = gmsh_table[gmshtype].nn_tot;
    size_t dim_i   = element[i].dimension();
    size_t variant = element[i].variant();
    map_dim = max(map_dim,dim_i);
    if (order == 0) {
      order = element[i].order();
    } else {
      check_macro (order == element[i].order(), "msh2geo: unexpected order="<<element[i].order());
    }
    element[i].resize(nn);
    for (size_t j = 0; j < nn; j++) {
      in >> element[i][j];
      element[i][j]--;
    }
    for (size_t subgeo_variant = 0; subgeo_variant <= variant; subgeo_variant++) { // set node dimension
      for (size_t loc_inod = reference_element_first_inod_by_variant(variant,element[i].order(),subgeo_variant), 
                  loc_nnod =  reference_element_last_inod_by_variant(variant,element[i].order(),subgeo_variant);
                     loc_inod < loc_nnod; loc_inod++) {
        node_subgeo_variant [element[i][loc_inod]] = subgeo_variant;
      }
    }
    // from gmsh to rheolef/geo local node renumbering: element[i][j] are modified
    msh2geo_node_renum (element[i], element[i].name(), element[i].order()); 
    if (domain_idx != 0) {
      domain_map[dim_i][domain_idx].push_back(i);
    }
  }
  array<size_t,reference_element__max_variant> loc_ndof_by_variant;
  reference_element_init_local_nnode_by_variant (order, loc_ndof_by_variant);
  // set dis_ie: by increasing variant order, for map_dim element only
  size_t last_index = 0;
  for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
              variant <  reference_element_last_variant_by_dimension(map_dim); variant++) {
    for (size_t ie = 0; ie < element.size(); ++ie) {
      if (element[ie].variant() != variant) continue;
      element[ie].setindex (last_index);
      switch (element[ie].dimension()) {
        case 0: element[ie].resize(1); element[ie][0] = last_index;  break;
        case 1: element[ie].edge_indirect(0).setindex(last_index); break;
        case 2: element[ie].face_indirect(0).setindex(last_index); break;
        default: break;
      }
      last_index++;
    }
  }
  // -------------------------------------------------
  // 2. node reordering, first pass
  // -------------------------------------------------
  // permut: first vertex, then edge, faces and inernal nodes, by increasing subgeo_dim 
  vector<size_t> old2new_inode (node.size(), std::numeric_limits<size_t>::max());
  if (true || !do_upgrade) {
    // 2.1. when no upgrade
    size_t new_inode = 0;
    for (size_t subgeo_variant = 0; subgeo_variant < reference_element_last_variant_by_dimension(map_dim); subgeo_variant++) {
      for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
        if (node_subgeo_variant [old_inode] != subgeo_variant) continue;
        old2new_inode[old_inode] = new_inode++;
      }
    }
  } else {
    // 2.2. when no upgrade: will be renumbered after side creation
    for (size_t inode = 0; inode < node.size(); inode++) {
      old2new_inode[inode] = inode;
    }
  }
  for (size_t inode = 0; inode < node.size(); inode++) {
    check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
  }
  // inv permut
  vector<size_t> new2old_inode (node.size(), std::numeric_limits<size_t>::max());
  for (size_t inode = 0; inode < node.size(); inode++) {
    new2old_inode[old2new_inode[inode]] = inode;
  }
  for (size_t inode = 0; inode < node.size(); inode++) {
    check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
  }
  // for convenience, apply the new numbering to nodes and elements
  {
    vector<point> new_node (node.size());
    for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
      new_node [old2new_inode[old_inode]] = node [old_inode];
    }
    for (size_t inode = 0; inode < node.size(); inode++) {
      node [inode] = new_node [inode];
    }
  }
  for (size_t ie = 0; ie < element.size(); ++ie) {
    for (size_t i = 0; i < element[ie].size(); ++i) {
      element[ie][i] = old2new_inode[element[ie][i]];
    }
  }
  // ------------------------------------------------------------------------
  // 3.1) compute all edges
  // ------------------------------------------------------------------------
  vector<index_set> edge_ball (node.size());
  vector<geo_element_side> edge;
  if (do_upgrade && map_dim >= 2) {
    list  <geo_element_side> edge_list;
    size_t nedg = 0;
    for (size_t ie = 0; ie < element.size(); ++ie) {
      if (element[ie].dimension() != map_dim) continue;
      for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
        size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
        size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
        size_t inod0 = element[ie][iv0];
        size_t inod1 = element[ie][iv1];
        index_set iedge_set = edge_ball[inod0];
        iedge_set.inplace_intersection (edge_ball[inod1]);
        check_macro (iedge_set.size() <= 1,
          "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
        if (iedge_set.size() == 1) continue; // edge already exists
        edge_ball[inod0].insert (nedg);
        edge_ball[inod1].insert (nedg);
        geo_element_side new_edge;
        new_edge.setname('e');
        new_edge[0] = inod0;
        new_edge[1] = inod1;
        new_edge.setindex(nedg);
        edge_list.push_back (new_edge);
        nedg++;
      }
    }
    edge.resize (edge_list.size());
    std::copy (edge_list.begin(), edge_list.end(), edge.begin());
  }
  // propagate edge numbering inside elements
  if (do_upgrade && map_dim >= 2) {
    for (size_t ie = 0; ie < element.size(); ++ie) {
      if (element[ie].dimension() != map_dim) continue;
      for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
        size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
        size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
        size_t inod0 = element[ie][iv0];
        size_t inod1 = element[ie][iv1];
        index_set iedge_set = edge_ball[inod0];
        iedge_set.inplace_intersection (edge_ball[inod1]);
        check_macro (iedge_set.size() == 1,
          "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
        size_t iedge = *(iedge_set.begin());
        element[ie].edge_indirect(iedg).setindex            (iedge);
        element[ie].edge_indirect(iedg).setorientation((edge[iedge][0] == inod0) ? 1 : -1);
      }
    }
  }
  // ------------------------------------------------------------------------
  // 3.2) compute all faces
  // ------------------------------------------------------------------------
  vector<index_set> face_ball (node.size());
  vector<geo_element_side> face;
  if (do_upgrade && map_dim >= 3) {
    list  <geo_element_side> face_list;
    size_t nfac = 0;
    for (size_t ie = 0; ie < element.size(); ++ie) {
      if (element[ie].dimension() != map_dim) continue;
      for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
        size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
        size_t inod0 = element[ie][iv0];
        index_set iface_set = face_ball[inod0];
        for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
          size_t ivj   = element[ie].face_local_vertex(loc_ifac, j);
          size_t inodj = element[ie][ivj];
          iface_set.inplace_intersection (face_ball[inodj]);
        }
        check_macro (iface_set.size() <= 1,
          "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
        if (iface_set.size() == 1) continue; // face already exists
        geo_element_side new_face;
        char face_name = (element[ie].face_n_vertex(loc_ifac) == 3) ? 't' : 'q';
        new_face.setname(face_name);
        new_face.setindex(nfac);
        for (size_t j = 0; j < element[ie].face_n_vertex(loc_ifac); ++j) {
          size_t ivj   = element[ie].face_local_vertex(loc_ifac, j);
          size_t inodj = element[ie][ivj];
          face_ball[inodj].insert (nfac);
          new_face[j] = inodj;
        }
        face_list.push_back (new_face);
        nfac++;
      }
    }
    face.resize (face_list.size());
    std::copy (face_list.begin(), face_list.end(), face.begin());
  }
  // propagate face numbering inside elements
  if (do_upgrade && map_dim >= 3) {
    for (size_t ie = 0; ie < element.size(); ++ie) {
      if (element[ie].dimension() != map_dim) continue;
      for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
        size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
        size_t inod0 = element[ie][iv0];
        index_set iface_set = face_ball[inod0];
        for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
          size_t ivj   = element[ie].face_local_vertex(loc_ifac, j);
          size_t inodj = element[ie][ivj];
          iface_set.inplace_intersection (face_ball[inodj]);
        }
        check_macro (iface_set.size() == 1,
          "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
        size_t iface = *(iface_set.begin());
        orientation_type orient;
        shift_type shift;
        if (element[ie].face_n_vertex(loc_ifac) == 3) {
          geo_element_get_orientation_and_shift (
            face[iface],
            element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
            element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
            element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
            orient, shift);
        } else {
          geo_element_get_orientation_and_shift (
            face[iface],
            element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
            element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
            element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
            element[ie][element[ie].face_local_vertex(loc_ifac, 3)],
            orient, shift);
        }
        element[ie].face_indirect(loc_ifac).setindex       (iface);
        element[ie].face_indirect(loc_ifac).setorientation (orient);
        element[ie].face_indirect(loc_ifac).setshift       (shift);
      }
    }
  }
  // ------------------------------------------------------------------------
  // 3.4) count all elements
  // ------------------------------------------------------------------------
  array<size_t,reference_element__max_variant>    size_by_variant;
  array<size_t,reference_element__max_variant+1> first_by_variant;
   size_by_variant.fill (0);
  first_by_variant.fill (0);
  if (true) {
    size_t n_vertex = 0;
    for (size_t inode = 0; inode < node.size(); inode++) {
      if (node_subgeo_variant [inode] == reference_element__p)  n_vertex++;
    }
    size_by_variant [reference_element__p] = n_vertex;
    if (map_dim >= 2 && edge.size() != 0) {
      size_by_variant [reference_element__e] = edge.size();
    }
    if (map_dim >= 3 && face.size() != 0) {
      for (size_t loc_ifac = 0; loc_ifac < face.size(); ++loc_ifac) {
        size_by_variant [face[loc_ifac].variant()]++;
      }
    }
    for (size_t ie = 0; ie < element.size(); ++ie) {
      if (element[ie].dimension() != map_dim) continue;
      size_t variant = element[ie].variant();
      size_by_variant [variant]++;
    }
    for (size_t dim = 0; dim <= 3; ++dim) {
      for (size_t variant = reference_element_first_variant_by_dimension(dim);
                  variant <  reference_element_last_variant_by_dimension(dim); variant++) {
        first_by_variant [variant+1] = size_by_variant [variant];
      }
    }
  }
  // -------------------------------------------------
  // 4. node reordering, second pass
  // -------------------------------------------------
  // need edges & faces with index & orient
  if (do_upgrade) {
    std::fill (old2new_inode.begin(), old2new_inode.end(), std::numeric_limits<size_t>::max());
    std::fill (new2old_inode.begin(), new2old_inode.end(), std::numeric_limits<size_t>::max());
    for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
                variant <  reference_element_last_variant_by_dimension(map_dim); variant++) {
      for (size_t ie = 0; ie < element.size(); ++ie) {
        if (element[ie].variant() != variant) continue;
        const geo_element& old_K = element[ie];
        std::vector<size_t> new_K;
        numbering_Pk_dis_idof (
          order, size_by_variant, first_by_variant, loc_ndof_by_variant, old_K, new_K);
        for (size_t loc_inod = 0; loc_inod < old_K.size(); loc_inod++) {
          size_t old_inod = old_K [loc_inod];
          size_t new_inod = new_K [loc_inod];
          if (old2new_inode [old_inod] != std::numeric_limits<size_t>::max()) continue; // already done
          old2new_inode [old_inod] = new_inod;
        }
      }
    }
  } else {
    // when no upgrade:
    for (size_t inode = 0; inode < node.size(); inode++) {
      old2new_inode[inode] = inode;
    }
  }
  for (size_t inode = 0; inode < node.size(); inode++) {
    check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
  }
  // inv permut
  for (size_t inode = 0; inode < node.size(); inode++) {
    new2old_inode[old2new_inode[inode]] = inode;
  }
  for (size_t inode = 0; inode < node.size(); inode++) {
    check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
  }
#ifdef TODO
#endif // TODO
  // ----------------------------------
  // 5. output geo
  // ---------------------------------- 
  // 5.1. output the mesh
  size_t version = 4;
  static const char* geo_variant_name [reference_element__max_variant] = {
    "points",
    "edges",
    "triangles",
    "quadrangles",
    "tetrahedra",
    "prisms",
    "hexahedra"
  };
  out << setprecision(numeric_limits<double>::digits10) 
      << "#!geo" << endl
      << endl
      << "mesh" << endl
      << version << endl
      << "header" << endl
      << " dimension\t" << dim << endl;
  if (sys_coord_name != "cartesian") {
      out << " coordinate_system " << sys_coord_name << endl;
  }
  if (order != 1) {
      out << " order\t\t" << order << endl;
  }
  out << " nodes\t\t" << node.size() << endl;
  
  if (map_dim > 0) {
    for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
                variant <  reference_element_last_variant_by_dimension(map_dim); variant++) {
      if (size_by_variant[variant] > 0) {
        out << " " << geo_variant_name [variant] << "\t" << size_by_variant[variant] << endl;
      }
    }
  }
  if (map_dim >= 3 &&          size_by_variant[reference_element__t] != 0) {
    out << " triangles\t"   << size_by_variant[reference_element__t] << endl;
  }
  if (map_dim >= 3 &&          size_by_variant[reference_element__q] != 0) {
    out << " quadrangles\t" << size_by_variant[reference_element__q] << endl;
  }
  if (map_dim >= 2 &&          size_by_variant[reference_element__e] != 0) {
    out << " edges\t"       << size_by_variant[reference_element__e] << endl;
  }
  out << "end header" << endl
      << endl;
  // nodes:
  for (size_t inode = 0; inode < node.size(); inode++) {
    for (size_t i = 0; i < dim; ++i) {
      out << node[new2old_inode[inode]][i];
      if (i+1 != dim) out << " ";
    }
    out << endl;
  }
  out << endl;
  // elements:
  for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
              variant <  reference_element_last_variant_by_dimension(map_dim); variant++) {
    for (size_t ie = 0; ie < element.size(); ie++) {
      if (element[ie].variant() != variant) continue;
      if (!do_upgrade) {
        if (element[ie].name() != 'e' || element[ie].order() > 1) {
          out << element[ie].name() << "\t";
        }
        if (element[ie].order() > 1) {
          out << "p" << element[ie].order() << " ";
        }
        for (size_t iloc = 0, nloc = element[ie].size(); iloc < nloc; iloc++) {
          out << element[ie][iloc];
          if (iloc+1 != nloc) out << " ";
        }
      } else {
        // upgrade: truncate inodes from pk to p1
        out << element[ie].name() << "\t";
        for (size_t iloc = 0, nloc = element[ie].n_vertex(); iloc < nloc; iloc++) {
          out << element[ie][iloc];
          if (iloc+1 != nloc) out << " ";
        }
      }
      out << endl;
    }
  }
  out << endl;
  // faces: no-empty when upgrade & map_dim >= 3
  for (size_t variant = reference_element_first_variant_by_dimension(2);
              variant <  reference_element_last_variant_by_dimension(2); variant++) {
    for (size_t ie = 0; ie < face.size(); ie++) {
      if (face[ie].variant() != variant) continue;
      out << face[ie].name() << "\t";
      for (size_t iloc = 0, nloc = face[ie].n_vertex(); iloc < nloc; iloc++) {
        out << face[ie][iloc];
        if (iloc+1 != nloc) out << " ";
      }
      out << endl;
    }
  }
  // edges: no-empty when upgrade & map_dim >= 2
  for (size_t ie = 0, ne = edge.size(); ie < ne; ++ie) {
      out << "e\t" << old2new_inode[edge[ie][0]]
          << " "   << old2new_inode[edge[ie][1]] << endl;
  }
  out << endl;
  //
  // 5.2. output domains
  //
  for (size_t d = 0; d <= 3; ++d) {
    if (!do_upgrade) {
      put_domain_noupgrade (out, map_dim, d, domain_map[d],  phys, element);
    } else {
      put_domain_upgrade   (out, map_dim, d, domain_map[d],  phys, element,
                                edge, edge_ball, face, face_ball);
    }
  }
}
void usage()
{
  cerr << "msh2geo: usage:" << endl
       << "msh2geo [-[no]upgrade][-rz|-zr]   in.msh > out.geo" << endl
       << "msh2geo [-[no]upgrade][-rz|-zr] < in.msh > out.geo" << endl;
  exit (1);
}
int main (int argc, char**argv) {
  // ----------------------------
  // scan the command line
  // ----------------------------
  string input_filename = "";
  std::string sys_coord_name = "cartesian";
  bool do_upgrade = true; // upgrade: put sides & truncate pk inodes to p1
  for (int i = 1; i < argc; i++) {
         if (strcmp (argv[i], "-rz") == 0)        sys_coord_name = "rz";
    else if (strcmp (argv[i], "-zr") == 0)        sys_coord_name = "zr";
    else if (strcmp (argv[i], "-cartesian") == 0) sys_coord_name = "cartesian";
    else if (strcmp (argv[i],   "-upgrade") == 0) do_upgrade = true;
    else if (strcmp (argv[i], "-noupgrade") == 0) do_upgrade = false;
    else if (argv[i][0] == '-') {
            cerr << "field: unknown option `" << argv[i] << "'" << endl;
            usage();
    } else {
      input_filename = argv[i];
    }
  }
  if (input_filename == "") {
    msh2geo (cin, cout, sys_coord_name, do_upgrade);
  } else {
    ifstream fin (input_filename.c_str());
    if (!fin.good()) { 
      cerr << "msh2geo: unable to read file \""<<input_filename<<"\"" << endl; exit (1);
    }
    msh2geo (fin, cout, sys_coord_name, do_upgrade);
  }
}