<filename>sfepy/discrete/fem/meshio.py
from __future__ import print_function
from __future__ import absolute_import
import sys
from copy import copy
import logging
import numpy as nm

from sfepy.base.base import (complex_types, dict_from_keys_init,
                             assert_, is_derived_class, ordered_iteritems,
                             insert_static_method, output, get_default,
                             get_default_attr, Struct, basestr)
from sfepy.base.ioutils import (skip_read_line, look_ahead_line, read_token,
                                read_array, pt, enc, dec,
                                edit_filename,
                                read_from_hdf5, write_to_hdf5,
                                HDF5ContextManager, get_or_create_hdf5_group)
import os.path as op
import six
from six.moves import range
import meshio as meshiolib
try:
    from meshio import CellBlock as meshio_Cells  # for meshio >= 4.0.3
except:
    from meshio import Cells as meshio_Cells  # for 4.0.3 > meshio > 4.0.0

_supported_formats = {
    # format name: IO class, suffix, modes[, variants]
    # modes: r = read, w = write, c = test cell groups, v = test vertex groups
    'abaqus': ('meshio', None, 'cv'),
    'exodus': ('meshio', None, 'v'),
    'gmsh': ('meshio', None, 'cv', ['gmsh4-binary', 'gmsh4-ascii',
                                    'gmsh2-binary', 'gmsh2-ascii']),
    'medit': ('meshio', None, 'cv'),
    'nastran': ('meshio', None, 'cv'),
    'vtk': ('meshio', None, 'cv', ['vtk-binnary', 'vtk-ascii']),
    'vtu': ('meshio', None, 'cv'),
    'med': ('meshio', None, 'cv'),
    'xdmf': ('meshio', None, 'cv'),
    'tetgen': ('meshio', None, ''),
    'ansys': ('ansys_cdb', '.cdb', 'r'),
    'hdf5': ('hdf5', '.h5', 'rwcv'),
    'hdf5-xdmf': ('hdf5-xdmf', '.h5x', 'rwcv'),
    'xyz': ('xyz', '.xyz', 'rw'),
    'comsol': ('comsol', '.txt', 'r'),
    'hmascii': ('hmascii', '.hmascii', 'r'),
    'gambit': ('gambit', '.neu', 'r'),
    'mesh3d': ('mesh3d', '.mesh3d', 'r'),
}

def update_supported_formats(formats):
    from meshio._helpers import reader_map, _writer_map,\
        extension_to_filetype

    f2e = {}
    for k, v in extension_to_filetype.items():
        f2e.setdefault(v, []).append(k)

    out = {}
    for format, info in formats.items():
        io, ext, _flag = info[:3]
        variants = info[3] if len(info) >= 4 else []
        for f in [format] + variants:
            if io == 'meshio':
                flag = _flag[:]
                if ext is None:
                    ext = f2e[format]
                if f in _writer_map:
                    flag = 'w' + flag
                if format in reader_map:
                    flag = 'r' + flag
                if not f == format:
                    flag = '*' + flag  # label format variants
            else:
                flag = _flag

            if not isinstance(ext, list):
                ext = [ext]

            out[f] = (io, ext, flag)

    return out

supported_formats = update_supported_formats(_supported_formats)
del _supported_formats


def output_mesh_formats(mode='r'):
    for key, vals in ordered_iteritems(supported_formats):
        if mode in vals[2]:
            output('%s (%s)' % (key,
                vals[1] if len(vals[1]) > 1 else vals[1][0]))


def split_conns_mat_ids(conns_in):
    """
    Split connectivities (columns except the last ones in `conns_in`) from cell
    groups (the last columns of `conns_in`).
    """
    conns, mat_ids = [], []
    for conn in conns_in:
        conn = nm.asarray(conn, dtype=nm.int32)
        conns.append(conn[:, :-1])
        mat_ids.append(conn[:, -1])

    return conns, mat_ids


def convert_complex_output(out_in):
    """
    Convert complex values in the output dictionary `out_in` to pairs of
    real and imaginary parts.
    """
    out = {}
    for key, val in six.iteritems(out_in):

        if val.data.dtype in  complex_types:
            rval = copy(val)
            rval.data = val.data.real
            out['real.%s' % key] = rval

            ival = copy(val)
            ival.data = val.data.imag
            out['imag.%s' % key] = ival

        else:
            out[key] = val

    return out


class MeshIO(Struct):
    """
    The abstract class for importing and exporting meshes.

    Read the docstring of the Mesh() class. Basically all you need to do is to
    implement the read() method::

        def read(self, mesh, **kwargs):
            nodes = ...
            ngroups = ...
            conns = ...
            mat_ids = ...
            descs = ...
            mesh._set_io_data(nodes, ngroups, conns, mat_ids, descs)
            return mesh

    See the Mesh class' docstring how the nodes, ngroups, conns, mat_ids and
    descs should look like. You just need to read them from your specific
    format from disk.

    To write a mesh to disk, just implement the write() method and use the
    information from the mesh instance (e.g. nodes, conns, mat_ids and descs)
    to construct your specific format.

    Optionally, subclasses can implement read_data() to read also computation
    results. This concerns mainly the subclasses with implemented write()
    supporting the 'out' kwarg.

    The default implementation od read_last_step() just returns 0. It should be
    reimplemented in subclasses capable of storing several steps.
    """
    format = None
    call_msg = 'called an abstract MeshIO instance!'

    def __init__(self, filename, **kwargs):
        Struct.__init__(self, filename=filename, **kwargs)
        self.set_float_format()

    def get_filename_trunk(self):
        if isinstance(self.filename, basestr):
            trunk = op.splitext(self.filename)[0]
        else:
            trunk = 'from_descriptor'

        return trunk

    def read_last_step(self):
        """The default implementation: just return 0 as the last step."""
        return 0

    def read_times(self, filename=None):
        """
        Read true time step data from individual time steps.

        Returns
        -------
        steps : array
            The time steps.
        times : array
            The times of the time steps.
        nts : array
            The normalized times of the time steps, in [0, 1].

        Notes
        -----
        The default implementation returns empty arrays.
        """
        aux = nm.array([0.0], dtype=nm.float64)
        return aux.astype(nm.int32), aux, aux

    def read(self, mesh, omit_facets=False, **kwargs):
        raise ValueError(MeshIO.call_msg)

    def write(self, filename, mesh, **kwargs):
        raise ValueError(MeshIO.call_msg)

    def read_data(self, step, filename=None, cache=None):
        raise ValueError(MeshIO.call_msg)

    def set_float_format(self, format=None):
        self.float_format = get_default(format, '%e')

    def get_vector_format(self, dim):
        return ' '.join([self.float_format] * dim)

class UserMeshIO(MeshIO):
    """
    Special MeshIO subclass that enables reading and writing a mesh using a
    user-supplied function.
    """
    format = 'function'

    def __init__(self, filename, **kwargs):
        assert_(hasattr(filename, '__call__'))
        self.function = filename

        MeshIO.__init__(self, filename='function:%s' % self.function.__name__,
                        **kwargs)

    def get_filename_trunk(self):
        return self.filename

    def read(self, mesh, *args, **kwargs):
        aux = self.function(mesh, mode='read')
        if aux is not None:
            mesh = aux

        self.filename = mesh.name

        return mesh

    def write(self, filename, mesh, *args, **kwargs):
        self.function(mesh, mode='write')


def _suppress_meshio_warnings(f):
    def __suppress_meshio_warnings(*args, **kwargs):
        logger = logging.getLogger()
        level = logger.level
        logger.setLevel(logging.ERROR)
        out = f(*args, **kwargs)
        logger.setLevel(level)

        return out

    return __suppress_meshio_warnings

def _decorate_all(module, decorator):
    import types
    for name in dir(module):
        obj = getattr(module, name)
        if isinstance(obj, types.FunctionType):
            setattr(module, name, decorator(obj))

_decorate_all(meshiolib, _suppress_meshio_warnings)
del _decorate_all, _suppress_meshio_warnings


class MeshioLibIO(MeshIO):
    format = 'meshio'

    cell_types = {
        ('hexahedron', 3): '3_8',
        ('tetra', 3): '3_4',
        ('triangle', 3): '2_3',
        ('triangle', 2): '2_3',
        ('quad', 3): '2_4',
        ('quad', 2): '2_4',
        ('line', 3): '3_2',
        ('line', 2): '2_2',
        ('line', 1): '1_2',
    }

    def __init__(self, filename, file_format=None, **kwargs):
        MeshIO.__init__(self, filename=filename, **kwargs)
        from meshio._helpers import _filetype_from_path
        import pathlib

        if file_format is None:
            file_format = _filetype_from_path(pathlib.Path(filename))

        self.file_format = file_format

    def read_bounding_box(self, ret_dim=False):
        m = meshiolib.read(self.filename, file_format=self.file_format)

        bbox = nm.vstack([nm.amin(m.points, 0), nm.amax(m.points, 0)])

        if ret_dim:
            return bbox, m.points.shape[1]
        else:
            return bbox

    @staticmethod
    def _get_dimension(points):
        dim = nm.sum(nm.max(points, axis=0)
                     - nm.min(points, axis=0) > 1e-15)
        return dim

    def read_dimension(self):
        m = meshiolib.read(self.filename, file_format=self.file_format)
        dim = self._get_dimension(m.points)

        return dim

    def read(self, mesh, omit_facets=False, **kwargs):
        m = meshiolib.read(self.filename, file_format=self.file_format)
        dim = self._get_dimension(m.points)

        ngkey = None
        for k in m.point_data.keys():
            if k == 'node_groups' or k.endswith(':ref'):
                ngkey = k
                break

        if ngkey is not None:
            ngroups = nm.asarray(m.point_data[ngkey]).flatten()
        elif hasattr(m, 'point_sets') and len(m.point_sets) > 0:
            ngroups = nm.zeros((len(m.points),), dtype=nm.int32)
            keys = list(m.point_sets.keys())
            keys.sort()
            try:
                ngrps = [int(ii) for ii in keys]
            except:
                ngrps = nm.arange(len(keys)) + 1

            for ik, k in enumerate(keys):
                ngroups[m.point_sets[k]] = ngrps[ik]
        else:
            ngroups = None

        cells, cgroups, cell_types = [], [], []

        # meshio.__version__ > 3.3.2
        cgkey = None
        for k in list(m.cell_data.keys()):
            if k == 'mat_id' or k.endswith(':ref'):
                cgkey = k
                break

        if cgkey is not None:
            cgdata = m.cell_data[cgkey]
        elif len(m.cell_sets) > 0:
            cgdata = []
            keys = list(m.cell_sets.keys())
            keys.sort()
            try:
                cgrps = [int(ii) for ii in keys]
            except:
                cgrps = nm.arange(len(keys)) + 1

            for ic, c in enumerate(m.cells):
                cgdata0 = nm.zeros((len(c.data),), dtype=nm.int32)
                for ik, k in enumerate(keys):
                    cgdata0[m.cell_sets[k][ic]] = cgrps[ik]
                cgdata.append(cgdata0)
        else:
            cgdata = None

        for ic, c in enumerate(m.cells):
            cells.append(c.data)
            cell_types.append(self.cell_types[(c.type, dim)])

            if cgdata is not None:
                cgroups.append(nm.asarray(cgdata[ic]).flatten())
            else:
                cgroups.append(nm.ones((len(c.data),), dtype=nm.int32))

        mesh._set_io_data(m.points[:,:dim], ngroups,
                          cells, cgroups, cell_types)

        output('number of vertices: %d' % m.points.shape[0])
        output('number of cells:')
        for ii, k in enumerate(cell_types):
            output('  %s: %d' % (k, cells[ii].shape[0]))

        return mesh

    def write(self, filename, mesh, out=None, ts=None, **kwargs):
        inv_cell_types = {v: k for k, v in self.cell_types.items()}

        coors, ngroups, conns, _, descs = mesh._get_io_data()

        out = {} if out is None else out

        point_data = {k: v.data for k, v in out.items() if v.mode == 'vertex'}
        cell_data_keys = [k for k, v in out.items() if v.mode == 'cell']

        if self.file_format in ['vtk', 'vtu']:
            ngkey = 'node_groups'
            cgkey = 'mat_id'
        else:
            ngkey = '%s:ref' % self.file_format
            cgkey = '%s:ref' % self.file_format

        point_data[ngkey] = ngroups
        point_sets = {str(k): nm.where(ngroups == k)[0]
            for k in nm.unique(ngroups)}

        cmesh = mesh.cmesh
        cell_groups = cmesh.cell_groups
        cgrps = nm.unique(cell_groups)

        # meshio.__version__ > 3.3.2
        cells = []
        cgroups = [ ]
        cell_data = {k: [] for k in cell_data_keys}
        cell_sets = {str(k): [] for k in cgrps}
        for ii, desc in enumerate(descs):
            cells.append(meshio_Cells(type=inv_cell_types[desc][0],
                                      data=conns[ii]))
            cidxs = nm.where(cmesh.cell_types == cmesh.key_to_index[desc])
            cidxs = cidxs[0].astype(nm.uint32)

            cgroups.append(cell_groups[cidxs])
            for k in cell_data_keys:
                cell_data[k].append(out[k].data[cidxs, 0, :, 0])

            for k in cgrps:
                idxs = nm.where(cell_groups[cidxs] == k)[0]
                cell_sets[str(k)].append(cidxs[idxs])

        cell_data[cgkey] = cgroups

        meshiolib.write_points_cells(filename, coors, cells,
                                     point_data=point_data,
                                     point_sets=point_sets,
                                     cell_data=cell_data,
                                     cell_sets=cell_sets,
                                     file_format=self.file_format)

    def read_data(self, step, filename=None, cache=None):
        m = meshiolib.read(self.filename, file_format=self.file_format)
        dim = self._get_dimension(m.points)

        def _fix_shape(data):
            if data.ndim == 2:
                data = data[:, :dim]

            elif data.ndim == 3:
                data = data[:, None, ...]

            return data

        out = {}
        for key, data in m.point_data.items():
            aux = _fix_shape(data).astype(nm.float64)
            out[key] = Struct(name=key, mode='vertex', data=aux)

        for key, data in m.cell_data.items():
            aux = _fix_shape(data[0]).astype(nm.float64)
            out[key] = Struct(name=key, mode='cell', data=aux)

        return out

class ComsolMeshIO(MeshIO):
    format = 'comsol'

    def _read_commented_int(self):
        return int(skip_read_line(self.fd).split('#')[0])

    def _skip_comment(self):
        read_token(self.fd)
        self.fd.readline()

    def read(self, mesh, **kwargs):

        self.fd = fd = open(self.filename, 'r')
        mode = 'header'

        coors = conns = None
        while 1:
            if mode == 'header':
                line = skip_read_line(fd)

                n_tags = self._read_commented_int()
                for ii in range(n_tags):
                    skip_read_line(fd)
                n_types = self._read_commented_int()
                for ii in range(n_types):
                    skip_read_line(fd)

                skip_read_line(fd)
                assert_(skip_read_line(fd).split()[1] == 'Mesh')
                skip_read_line(fd)
                dim = self._read_commented_int()
                assert_((dim == 2) or (dim == 3))
                n_nod = self._read_commented_int()
                i0 = self._read_commented_int()
                mode = 'points'

            elif mode == 'points':
                self._skip_comment()
                coors = read_array(fd, n_nod, dim, nm.float64)
                mode = 'cells'

            elif mode == 'cells':

                n_types = self._read_commented_int()
                conns = []
                descs = []
                mat_ids = []
                for it in range(n_types):
                    t_name = skip_read_line(fd).split()[1]
                    n_ep = self._read_commented_int()
                    n_el = self._read_commented_int()

                    self._skip_comment()
                    aux = read_array(fd, n_el, n_ep, nm.int32)
                    if t_name == 'tri':
                        conns.append(aux)
                        descs.append('2_3')
                        is_conn = True
                    elif t_name == 'quad':
                        # Rearrange element node order to match SfePy.
                        aux = aux[:,(0,1,3,2)]
                        conns.append(aux)
                        descs.append('2_4')
                        is_conn = True
                    elif t_name == 'hex':
                        # Rearrange element node order to match SfePy.
                        aux = aux[:,(0,1,3,2,4,5,7,6)]
                        conns.append(aux)
                        descs.append('3_8')
                        is_conn = True
                    elif t_name == 'tet':
                        conns.append(aux)
                        descs.append('3_4')
                        is_conn = True
                    else:
                        is_conn = False

                    # Skip parameters.
                    n_pv = self._read_commented_int()
                    n_par = self._read_commented_int()
                    for ii in range(n_par):
                        skip_read_line(fd)

                    n_domain = self._read_commented_int()
                    assert_(n_domain == n_el)
                    if is_conn:
                        self._skip_comment()
                        mat_id = read_array(fd, n_domain, 1, nm.int32)
                        mat_ids.append(mat_id.squeeze())
                    else:
                        for ii in range(n_domain):
                            skip_read_line(fd)

                    # Skip up/down pairs.
                    n_ud = self._read_commented_int()
                    for ii in range(n_ud):
                        skip_read_line(fd)
                break

        fd.close()
        self.fd = None

        mesh._set_io_data(coors, None, conns, mat_ids, descs)

        return mesh

    def write(self, filename, mesh, out=None, **kwargs):

        def write_elements(fd, ig, conn, mat_ids, type_name,
                           npe, format, norder, nm_params):
            fd.write("# Type #%d\n\n" % ig)
            fd.write("%s # type name\n\n\n" % type_name)
            fd.write("%d # number of nodes per element\n" % npe)
            fd.write("%d # number of elements\n" % conn.shape[0])
            fd.write("# Elements\n")
            for ii in range(conn.shape[0]):
                nn = conn[ii] # Zero based
                fd.write(format % tuple(nn[norder]))
            fd.write("\n%d # number of parameter values per element\n"
                     % nm_params)
            # Top level always 0?
            fd.write("0 # number of parameters\n")
            fd.write("# Parameters\n\n")
            fd.write("%d # number of domains\n"
                     % sum([mi.shape[0] for mi in mat_ids]))
            fd.write("# Domains\n")
            for mi in mat_ids:
                # Domains in comsol have to be > 0
                if (mi <= 0).any():
                    mi += mi.min() + 1
                for dom in mi:
                    fd.write("%d\n" % abs(dom))
            fd.write("\n0 # number of up/down pairs\n")
            fd.write("# Up/down\n")

        fd = open(filename, 'w')

        coors, ngroups, conns, mat_ids, desc = mesh._get_io_data()

        n_nod, dim = coors.shape

        # Header
        fd.write("# Created by SfePy\n\n\n")
        fd.write("# Major & minor version\n")
        fd.write("0 1\n")
        fd.write("1 # number of tags\n")
        fd.write("# Tags\n")
        fd.write("2 m1\n")
        fd.write("1 # number of types\n")
        fd.write("# Types\n")
        fd.write("3 obj\n\n")

        # Record
        fd.write("# --------- Object 0 ----------\n\n")
        fd.write("0 0 1\n") # version unused serializable
        fd.write("4 Mesh # class\n")
        fd.write("1 # version\n")
        fd.write("%d # sdim\n" % dim)
        fd.write("%d # number of mesh points\n" % n_nod)
        fd.write("0 # lowest mesh point index\n\n") # Always zero in SfePy

        fd.write("# Mesh point coordinates\n")

        format = self.get_vector_format(dim) + '\n'
        for ii in range(n_nod):
            nn = tuple(coors[ii])
            fd.write(format % tuple(nn))

        fd.write("\n%d # number of element types\n\n\n" % len(conns))

        for ig, conn in enumerate(conns):
            if (desc[ig] == "2_4"):
                write_elements(fd, ig, conn, mat_ids,
                               "4 quad", 4, "%d %d %d %d\n", [0, 1, 3, 2], 8)

            elif (desc[ig] == "2_3"):
                # TODO: Verify number of parameters for tri element
                write_elements(fd, ig, conn, mat_ids,
                               "3 tri", 3, "%d %d %d\n", [0, 1, 2], 4)

            elif (desc[ig] == "3_4"):
                # TODO: Verify number of parameters for tet element
                write_elements(fd, ig, conn, mat_ids,
                               "3 tet", 4, "%d %d %d %d\n", [0, 1, 2, 3], 16)

            elif (desc[ig] == "3_8"):
                write_elements(fd, ig, conn, mat_ids,
                               "3 hex", 8, "%d %d %d %d %d %d %d %d\n",
                               [0, 1, 3, 2, 4, 5, 7, 6], 24)

            else:
                raise ValueError('unknown element type! (%s)' % desc[ig])

        fd.close()

        if out is not None:
            for key, val in six.iteritems(out):
                raise NotImplementedError


class HDF5MeshIO(MeshIO):
    format = "hdf5"

    import string
    _all = ''.join(map(chr, list(range(256))))
    _letters = string.ascii_letters + string.digits + '_'
    _rubbish = ''.join([ch for ch in set(_all) - set(_letters)])
    if sys.version_info[0] >= 3:
        _tr = str.maketrans(_rubbish, '_' * len(_rubbish))
    else:
        _tr = string.maketrans(_rubbish, '_' * len(_rubbish))

    @staticmethod
    def read_mesh_from_hdf5(filename, group=None, mesh=None):
        """
        Read the mesh from a HDF5 file.

        filename: str or tables.File
            The HDF5 file to read the mesh from.
        group: tables.group.Group or str, optional
            The HDF5 file group to read the mesh from.
            If None, the root group is used.
        mesh: sfepy.dicrete.fem.Mesh or None
            If None, the new mesh is created and returned, otherwise
            content of this argument is replaced by the read mesh.

        Returns
        -------
        sfepy.dicrete.fem.Mesh
            readed mesh
        """
        with HDF5ContextManager(filename, mode='r') as fd:
            if group is None:
                group = fd.root
            elif not isinstance(group, pt.group.Group):
                group = fd.get_node(group)

            set_shape_info = mesh is None
            if mesh is None:
                from .mesh import Mesh
                mesh = Mesh('mesh')

            mesh.name = dec(group.name.read())
            coors = group.coors.read()
            ngroups = group.ngroups.read()

            n_gr = group.n_gr.read()

            conns = []
            descs = []
            mat_ids = []
            for ig in range(n_gr):
                gr_name = 'group%d' % ig
                conn_group = group._f_get_child(gr_name)
                conns.append(conn_group.conn.read())
                mat_ids.append(conn_group.mat_id.read())
                descs.append(dec(conn_group.desc.read()))

            nodal_bcs = {}
            try:
                node_sets_groups = group.node_sets

            except:
                pass

            else:
                for group in node_sets_groups:
                    key = dec(group.key.read())
                    nods = group.nods.read()
                    nodal_bcs[key] = nods

            mesh._set_io_data(coors, ngroups, conns, mat_ids, descs,
                              nodal_bcs=nodal_bcs)

            if set_shape_info:
                mesh._set_shape_info()
        return mesh

    @staticmethod
    def write_mesh_to_hdf5(filename, group, mesh, force_3d=False):
        """
        Write mesh to a hdf5 file.

        filename: str or tables.File
            The HDF5 file to write the mesh to.
        group: tables.group.Group or None or str
            The HDF5 file group to write the mesh to.
            If None, the root group is used.
            The group can be given as a path from root, e.g. /path/to/mesh
        mesh: sfepy.dicrete.fem.Mesh
            The mesh to write.
        """
        with HDF5ContextManager(filename, mode='w') as fd:
            if group is None:
                group = fd.root
            elif not isinstance(group, pt.group.Group):
                group = get_or_create_hdf5_group(fd, group)

            coors, ngroups, conns, mat_ids, descs = mesh._get_io_data()
            if force_3d and coors.shape[1] == 2:
                coors = nm.hstack([coors, nm.zeros((len(coors), 1))])
            fd.create_array(group, 'name', enc(mesh.name), 'name')
            fd.create_array(group, 'coors', coors, 'coors')
            fd.create_array(group, 'ngroups', ngroups, 'ngroups')
            fd.create_array(group, 'n_gr', len(conns), 'n_gr')
            for ig, conn in enumerate(conns):
                conn_group = fd.create_group(group, 'group%d' % ig,
                                             'connectivity group')
                fd.create_array(conn_group, 'conn', conn, 'connectivity')
                fd.create_array(conn_group, 'mat_id', mat_ids[ig],
                                'material id')
                fd.create_array(conn_group, 'desc', enc(descs[ig]),
                                'element Type')

            node_sets_groups = fd.create_group(group, 'node_sets',
                                               'node sets groups')
            ii = 0
            for key, nods in six.iteritems(mesh.nodal_bcs):
                group = fd.create_group(node_sets_groups, 'group%d' % ii,
                                        'node sets group')
                fd.create_array(group, 'key', enc(key), 'key')
                fd.create_array(group, 'nods', nods, 'nods')
                ii += 1

    def read_dimension(self, ret_fd=False):
        fd = pt.open_file(self.filename, mode="r")

        dim = fd.root.mesh.coors.shape[1]

        if ret_fd:
            return dim, fd

        else:
            fd.close()
            return dim

    def read_bounding_box(self, ret_fd=False, ret_dim=False):
        fd = pt.open_file(self.filename, mode="r")

        mesh_group = fd.root.mesh

        coors = mesh_group.coors.read()
        bbox = nm.vstack((nm.amin(coors, 0),
                          nm.amax(coors, 0)))

        if ret_dim:
            dim = coors.shape[1]

            if ret_fd:
                return bbox, dim, fd

            else:
                fd.close()
                return bbox, dim

        else:
            if ret_fd:
                return bbox, fd

            else:
                fd.close()
                return bbox

    def read(self, mesh=None, **kwargs):
        return self.read_mesh_from_hdf5(self.filename, '/mesh', mesh=mesh)

    @staticmethod
    def write_xdmf_file(filename, **kwargs):
        def get_path(node):
            _, fname = op.split(filename)
            return '%s:%s' % (fname, node._v_pathname)

        def get_data_dim(shape):
            if len(shape) == 4:
               return shape[2]
            if len(shape) == 2:
               return shape[1]
            else:
               return 1

        def data_item(data):
            dtype = data.dtype
            if nm.issubdtype(dtype, nm.integer):
                data_type = 'Int'
            elif nm.issubdtype(dtype, nm.floating):
                data_type = 'Float'
            else:
                raise ValueError('wrong data type! (%s)' % dtype)

            dim = get_data_dim(data.shape)
            sh = (data.shape[0], dim)
            ditem = et.Element('DataItem',
                               attrib={'DataType': data_type,
                                       'Dimensions': '%d %d' % sh,
                                       'Format': 'HDF',
                                       'Precision': str(dtype.itemsize)})
            ditem.text = get_path(data)

            return ditem

        def attr_item(data, center=None, name=None):
            if isinstance(data, pt.Group):
                if center is None and 'mode' in data:
                    mode = data.mode.read().decode('ascii')
                    if mode == 'custom':
                        return None
                    center = {'vertex': 'Node',
                              'cell': 'Cell'}[mode]
                if name is None and 'dname' in data:
                    name = data.dname.read().decode('ascii').lstrip('_')
                data = data.data

            dim = get_data_dim(data.shape)
            atype = {1: 'Scalar', 3: 'Vector', 6: 'Tensor6', 9: 'Tensor'}[dim]

            aitem = et.Element('Attribute',
                               attrib={'AttributeType': atype,
                                       'Center': center,
                                       'Name': name})
            aitem.append(data_item(data))

            return aitem

        import xml.etree.ElementTree as et
        from xml.dom import minidom

        topology_table = {
            '2_2': 'Line',
            '3_2': 'Line',
            '2_3': 'Triangle',
            '2_4': 'Quadrilateral',
            '3_4': 'Tetrahedron',
            '3_8': 'Hexahedron',
        }

        et_root = et.Element('Xdmf', attrib={'Version': '3.0'})
        if 'extra_data' in kwargs:
            for k, v in kwargs['extra_data'].items():
                d = et.SubElement(et_root, 'Information', attrib={'Name': k})
                d.text = str(v)

        et_domain = et.SubElement(et_root, 'Domain')

        with HDF5ContextManager(filename, mode='r') as fd:
            root = fd.root
            mesh = root.mesh
            name = mesh.name.read().decode('ascii')
            et_mesh = []
            geom = et.Element('Geometry', attrib={'GeometryType': 'XYZ'})
            geom.append(data_item(mesh.coors))
            et_mesh.append(geom)
            et_mesh.append(attr_item(mesh.ngroups, 'Node', 'node_groups'))

            n_gr = mesh.n_gr.read()
            for ig in range(n_gr):
                gr_name = 'group%d' % ig
                conn_group = mesh._f_get_child(gr_name)
                nc, nnd = conn_group.conn.shape
                ttype = topology_table[dec(conn_group.desc.read())]
                et_conn = et.Element('Topology',
                                     attrib={'NumberOfElements': str(nc),
                                             'TopologyType': ttype})
                et_conn.append(data_item(conn_group.conn))
                et_mesh.append(et_conn)
                et_mesh.append(attr_item(conn_group.mat_id, 'Cell', 'mat_id'))

            steps = [k for k in root if k._v_name.startswith('step')]
            et_ts = et.SubElement(et_domain, 'Grid',
                                    attrib={'Name': 'TimeSeries',
                                            'GridType': 'Collection',
                                            'CollectionType': 'Temporal'})

            for step in steps:
                istep = int(step._v_name[4:])
                et_grid = et.SubElement(et_ts, 'Grid',
                                        attrib={'Name': 'grid%d' % istep,
                                                'GridType': 'Uniform'})
                et.SubElement(et_grid, 'Time', attrib={'Value': '%d' % istep})
                et_grid.extend(et_mesh)

                for val in filter(lambda x: x._v_name.startswith('__'), step):
                    aitem = attr_item(val)
                    if aitem is not None:
                        et_grid.append(aitem)

        out = minidom.parseString(et.tostring(et_root)).toprettyxml(indent="  ")
        xdmf_filename = op.splitext(filename)[0] + '.xdmf'
        with open(xdmf_filename, 'w') as f:
            f.write(out[(out.find('\n') + 1):])

    def write(self, filename, mesh, out=None, ts=None, cache=None,
              xdmf=False, **kwargs):
        def expand_data_3d(data):
            expand_tab = {
                2: (3, [0, 1]),
                3: (6, [0, 1, 3]),
                4: (9, [0, 1, 3, 4]),
            }

            n, dim = data.shape
            if dim in expand_tab:
                dim3, order = expand_tab[dim]
                out = nm.zeros((n, dim3), dtype=data.dtype)
                out[:, order] = data
                return out
            else:
                return data

        from time import asctime

        if pt is None:
            raise ValueError('pytables not imported!')

        step = get_default_attr(ts, 'step', 0)
        if (step == 0) or not op.exists(filename):
            # A new file.
            with pt.open_file(filename, mode="w",
                              title="SfePy output file") as fd:
                mesh_group = fd.create_group('/', 'mesh', 'mesh')
                self.write_mesh_to_hdf5(fd, mesh_group, mesh, force_3d=xdmf)

                if ts is not None:
                    ts_group = fd.create_group('/', 'ts', 'time stepper')
                    fd.create_array(ts_group, 't0', ts.t0, 'initial time')
                    fd.create_array(ts_group, 't1', ts.t1, 'final time' )
                    fd.create_array(ts_group, 'dt', ts.dt, 'time step')
                    fd.create_array(ts_group, 'n_step', ts.n_step, 'n_step')

                tstat_group = fd.create_group('/', 'tstat',
                                              'global time statistics')
                fd.create_array(tstat_group, 'created', enc(asctime()),
                                'file creation time')
                fd.create_array(tstat_group, 'finished', enc('.' * 24),
                                'file closing time')

                fd.create_array(fd.root, 'last_step',
                                nm.array([step], dtype=nm.int32),
                                'last saved step')

        if out is not None:
            if ts is None:
                step, time, nt  = 0, 0.0, 0.0
            else:
                step, time, nt = ts.step, ts.time, ts.nt

            # Existing file.
            fd = pt.open_file(filename, mode="r+")

            step_group_name = 'step%d' % step
            if step_group_name in fd.root:
                raise ValueError('step %d is already saved in "%s" file!'
                                 ' Possible help: remove the old file or'
                                 ' start saving from the initial time.'
                                 % (step, filename))
            step_group = fd.create_group('/', step_group_name, 'time step data')

            ts_group = fd.create_group(step_group, 'ts', 'time stepper')
            fd.create_array(ts_group, 'step', step, 'step')
            fd.create_array(ts_group, 't', time, 'time')
            fd.create_array(ts_group, 'nt', nt, 'normalized time')

            name_dict = {}
            for key, val in six.iteritems(out):
                if xdmf and mesh.coors.shape[1] == 2:
                    data = expand_data_3d(val.data)
                else:
                    data = val.data

                group_name = '__' + key.translate(self._tr)
                data_group = fd.create_group(step_group, group_name,
                                             '%s data' % key)
                fd.create_array(data_group, 'dname', enc(key), 'data name')
                fd.create_array(data_group, 'mode', enc(val.mode), 'mode')
                name = val.get('name', 'output_data')
                fd.create_array(data_group, 'name', enc(name), 'object name')
                if val.mode == 'custom':
                    write_to_hdf5(fd, data_group, 'data', data,
                                  cache=cache,
                                  unpack_markers=getattr(val, 'unpack_markers',
                                                         False))
                    continue

                shape = val.get('shape', data.shape)
                dofs = val.get('dofs', None)
                if dofs is None:
                    dofs = [''] * nm.squeeze(shape)[-1]
                var_name = val.get('var_name', '')

                fd.create_array(data_group, 'data', data, 'data')
                fd.create_array(data_group, 'dofs', [enc(ic) for ic in dofs],
                                'dofs')
                fd.create_array(data_group, 'shape', shape, 'shape')
                fd.create_array(data_group, 'var_name',
                                enc(var_name), 'object parent name')
                if val.mode == 'full':
                    fd.create_array(data_group, 'field_name',
                                    enc(val.field_name), 'field name')

                name_dict[key] = group_name

            step_group._v_attrs.name_dict = name_dict
            fd.root.last_step[0] = step

            fd.remove_node(fd.root.tstat.finished)
            fd.create_array(fd.root.tstat, 'finished', enc(asctime()),
                            'file closing time')
            fd.close()

        if xdmf:
            self.write_xdmf_file(filename, **kwargs)

    def read_last_step(self, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")
        last_step = fd.root.last_step[0]
        fd.close()
        return last_step

    def read_time_stepper(self, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        try:
            ts_group = fd.root.ts
            out =  (ts_group.t0.read(), ts_group.t1.read(),
                    ts_group.dt.read(), ts_group.n_step.read())

        except:
            raise ValueError('no time stepper found!')

        finally:
            fd.close()

        return out

    def _get_step_group_names(self, fd):
        return sorted([name for name in fd.root._v_groups.keys()
                       if name.startswith('step')],
                      key=lambda name: int(name[4:]))

    def read_times(self, filename=None):
        """
        Read true time step data from individual time steps.

        Returns
        -------
        steps : array
            The time steps.
        times : array
            The times of the time steps.
        nts : array
            The normalized times of the time steps, in [0, 1].
        """
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode='r')

        steps = []
        times = []
        nts = []
        for gr_name in self._get_step_group_names(fd):
            ts_group = fd.get_node(fd.root, gr_name + '/ts')

            steps.append(ts_group.step.read())
            times.append(ts_group.t.read())
            nts.append(ts_group.nt.read())
        fd.close()

        steps = nm.asarray(steps, dtype=nm.int32)
        times = nm.asarray(times, dtype=nm.float64)
        nts = nm.asarray(nts, dtype=nm.float64)

        return steps, times, nts

    def _get_step_group(self, step, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        if step is None:
            step = int(self._get_step_group_names(fd)[0][4:])

        gr_name = 'step%d' % step
        try:
            step_group = fd.get_node(fd.root, gr_name)
        except:
            output('step %d data not found - premature end of file?' % step)
            fd.close()
            return None, None

        return fd, step_group

    def read_data(self, step, filename=None, cache=None):
        fd, step_group = self._get_step_group(step, filename=filename)
        if fd is None: return None

        out = {}
        for data_group in step_group:
            try:
                key = dec(data_group.dname.read())

            except pt.exceptions.NoSuchNodeError:
                continue

            mode = dec(data_group.mode.read())
            if mode == 'custom':
               out[key] = read_from_hdf5(fd, data_group.data, cache=cache)
               continue

            name = dec(data_group.name.read())
            data = data_group.data.read()
            dofs = tuple([dec(ic) for ic in data_group.dofs.read()])
            try:
                shape = tuple(int(ii) for ii in data_group.shape.read())

            except pt.exceptions.NoSuchNodeError:
                shape = data.shape

            if mode == 'full':
                field_name = dec(data_group.field_name.read())

            else:
                field_name = None

            out[key] = Struct(name=name, mode=mode, data=data,
                              dofs=dofs, shape=shape, field_name=field_name)

            if out[key].dofs == (-1,):
                out[key].dofs = None

        fd.close()

        return out

    def read_data_header(self, dname, step=None, filename=None):
        fd, step_group = self._get_step_group(step, filename=filename)
        if fd is None: return None

        groups = step_group._v_groups
        for name, data_group in six.iteritems(groups):
            try:
                key = dec(data_group.dname.read())

            except pt.exceptions.NoSuchNodeError:
                continue

            if key == dname:
                mode = dec(data_group.mode.read())
                fd.close()
                return mode, name

        fd.close()
        raise KeyError('non-existent data: %s' % dname)

    def read_time_history(self, node_name, indx, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        th = dict_from_keys_init(indx, list)
        for gr_name in self._get_step_group_names(fd):
            step_group = fd.get_node(fd.root, gr_name)
            data = step_group._f_get_child(node_name).data

            for ii in indx:
                th[ii].append(nm.array(data[ii]))

        fd.close()

        for key, val in six.iteritems(th):
            aux = nm.array(val)
            if aux.ndim == 4: # cell data.
                aux = aux[:,0,:,0]
            th[key] = aux

        return th

    def read_variables_time_history(self, var_names, ts, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        assert_((fd.root.last_step[0] + 1) == ts.n_step)

        ths = dict_from_keys_init(var_names, list)

        arr = nm.asarray
        for step in range(ts.n_step):
            gr_name = 'step%d' % step
            step_group = fd.get_node(fd.root, gr_name)
            name_dict = step_group._v_attrs.name_dict
            for var_name in var_names:
                data = step_group._f_get_child(name_dict[var_name]).data
                ths[var_name].append(arr(data.read()))

        fd.close()

        return ths


class HDF5XdmfMeshIO(HDF5MeshIO):
    format="hdf5-xdmf"

    def write(self, filename, mesh, out=None, ts=None, cache=None, **kwargs):
        HDF5MeshIO.write(self, filename, mesh, out=out, ts=ts, cache=cache,
                         xdmf=True, **kwargs)


class Mesh3DMeshIO(MeshIO):
    format = "mesh3d"

    def read(self, mesh, **kwargs):
        f = open(self.filename)
        # read the whole file:
        vertices = self._read_section(f, integer=False)
        tetras = self._read_section(f)
        hexes = self._read_section(f)
        prisms = self._read_section(f)
        tris = self._read_section(f)
        quads = self._read_section(f)

        # substract 1 from all elements, because we count from 0:
        conns = []
        mat_ids = []
        descs = []
        if len(tetras) > 0:
            conns.append(tetras - 1)
            mat_ids.append([0]*len(tetras))
            descs.append("3_4")
        if len(hexes) > 0:
            conns.append(hexes - 1)
            mat_ids.append([0]*len(hexes))
            descs.append("3_8")
        mesh._set_io_data(vertices, None, conns, mat_ids, descs)
        return mesh

    def read_dimension(self):
        return 3

    def _read_line(self, f):
        """
        Reads one non empty line (if it's a comment, it skips it).
        """
        l = f.readline().strip()
        while l == "" or l[0] == "#": # comment or an empty line
            l = f.readline().strip()
        return l

    def _read_section(self, f, integer=True):
        """
        Reads one section from the mesh3d file.

        integer ... if True, all numbers are passed to int(), otherwise to
            float(), before returning

        Some examples how a section can look like:

        2
        1 2 5 4 7 8 11 10
        2 3 6 5 8 9 12 11

        or

        5
        1 2 3 4     1
        1 2 6 5     1
        2 3 7 6     1
        3 4 8 7     1
        4 1 5 8     1

        or

        0

        """
        if integer:
            dtype=int
        else:
            dtype=float
        l = self._read_line(f)
        N = int(l)
        rows = []
        for i in range(N):
            l = self._read_line(f)
            row = nm.fromstring(l, sep=" ", dtype=dtype)
            rows.append(row)
        return nm.array(rows)

def mesh_from_groups(mesh, ids, coors, ngroups,
                     tris, mat_tris, quads, mat_quads,
                     tetras, mat_tetras, hexas, mat_hexas, remap=None):
    ids = nm.asarray(ids, dtype=nm.int32)
    coors = nm.asarray(coors, dtype=nm.float64)

    if remap is None:
        n_nod = coors.shape[0]
        remap = nm.zeros((ids.max()+1,), dtype=nm.int32)
        remap[ids] = nm.arange(n_nod, dtype=nm.int32)

    tris = remap[nm.array(tris, dtype=nm.int32)]
    quads = remap[nm.array(quads, dtype=nm.int32)]
    tetras = remap[nm.array(tetras, dtype=nm.int32)]
    hexas = remap[nm.array(hexas, dtype=nm.int32)]

    conns = [tris, quads, tetras, hexas]
    mat_ids = [nm.array(ar, dtype=nm.int32)
               for ar in [mat_tris, mat_quads, mat_tetras, mat_hexas]]
    descs = ['2_3', '2_4', '3_4', '3_8']

    # Remove empty groups.
    conns, mat_ids, descs = zip(*[(conns[ig], mat_ids[ig], descs[ig])
                                  for ig in range(4)
                                  if conns[ig].shape[0] > 0])

    mesh._set_io_data(coors, ngroups, conns, mat_ids, descs)
    return mesh


class HypermeshAsciiMeshIO(MeshIO):
    format = 'hmascii'

    def read(self, mesh, **kwargs):
        fd = open(self.filename, 'r')

        ids = []
        coors = []
        tetras = []
        mat_tetras = []
        hexas = []
        mat_hexas = []
        quads = []
        mat_quads = []
        trias = []
        mat_trias = []
        mat_id = 0

        for line in fd:
            if line and (line[0] == '*'):
                if line[1:10] == 'component':
                    line = line.strip()[11:-1].split(',')
                    mat_id = int(line[0])
                if line[1:5] == 'node':
                    line = line.strip()[6:-1].split(',')
                    ids.append(int(line[0]))
                    coors.append([float(coor) for coor in line[1:4]])

                elif line[1:7] == 'tetra4':
                    line = line.strip()[8:-1].split(',')
                    mat_tetras.append(mat_id)
                    tetras.append([int(ic) for ic in line[2:6]])

                elif line[1:6] == 'hexa8':
                    line = line.strip()[7:-1].split(',')
                    mat_hexas.append(mat_id)
                    hexas.append([int(ic) for ic in line[2:10]])

                elif line[1:6] == 'quad4':
                    line = line.strip()[7:-1].split(',')
                    mat_quads.append(mat_id)
                    quads.append([int(ic) for ic in line[2:6]])

                elif line[1:6] == 'tria3':
                    line = line.strip()[7:-1].split(',')
                    mat_trias.append(mat_id)
                    trias.append([int(ic) for ic in line[2:5]])
        fd.close()

        mesh = mesh_from_groups(mesh, ids, coors, None,
                                trias, mat_trias, quads, mat_quads,
                                tetras, mat_tetras, hexas, mat_hexas)

        return mesh

    def read_dimension(self):
        return 3

    def write(self, filename, mesh, out=None, **kwargs):
        raise NotImplementedError


class NEUMeshIO(MeshIO):
    format = 'gambit'

    def read_dimension(self, ret_fd=False):

        fd = open(self.filename, 'r')

        row = fd.readline().split()
        while 1:
            if not row: break
            if len(row) == 0: continue

            if (row[0] == 'NUMNP'):
                row = fd.readline().split()
                n_nod, n_el, dim = row[0], row[1], int(row[4])
                break

        if ret_fd:
            return dim, fd
        else:
            fd.close()
            return dim

    def read(self, mesh, **kwargs):

        el = {'3_8' : [], '3_4' : [], '2_4' : [], '2_3' : []}
        nod = []

        conns_in = []
        descs = []

        group_ids = []
        group_n_els = []
        groups = []
        nodal_bcs = {}

        fd = open(self.filename, 'r')

        row = fd.readline()
        while 1:
            if not row: break
            row = row.split()

            if len(row) == 0:
                row = fd.readline()
                continue

            if (row[0] == 'NUMNP'):
                row = fd.readline().split()
                n_nod, n_el, dim = int(row[0]), int(row[1]), int(row[4])

            elif (row[0] == 'NODAL'):
                row = fd.readline().split()
                while not(row[0] == 'ENDOFSECTION'):
                    nod.append(row[1:])
                    row = fd.readline().split()

            elif (row[0] == 'ELEMENTS/CELLS'):
                row = fd.readline().split()
                while not(row[0] == 'ENDOFSECTION'):
                    elid = [row[0]]
                    gtype = int(row[1])
                    if gtype == 6:
                        el['3_4'].append(row[3:]+elid)
                    elif gtype == 4:
                        rr = row[3:]
                        if (len(rr) < 8):
                            rr.extend(fd.readline().split())
                        el['3_8'].append(rr+elid)
                    elif gtype == 3:
                        el['2_3'].append(row[3:]+elid)
                    elif gtype == 2:
                        el['2_4'].append(row[3:]+elid)
                    row = fd.readline().split()

            elif (row[0] == 'GROUP:'):
                group_ids.append(row[1])
                g_n_el = int(row[3])
                group_n_els.append(g_n_el)
                name = fd.readline().strip()

                els = []
                row = fd.readline().split()
                row = fd.readline().split()
                while not(row[0] == 'ENDOFSECTION'):
                    els.extend(row)
                    row = fd.readline().split()
                if g_n_el != len(els):
                    msg = 'wrong number of group elements! (%d == %d)'\
                          % (n_el, len(els))
                    raise ValueError(msg)
                groups.append(els)

            elif (row[0] == 'BOUNDARY'):
                row = fd.readline().split()
                key = row[0]
                num = int(row[2])
                inod = read_array(fd, num, None, nm.int32) - 1
                nodal_bcs[key] = inod.squeeze()

                row = fd.readline().split()
                assert_(row[0] == 'ENDOFSECTION')

            row = fd.readline()

        fd.close()

        if int(n_el) != sum(group_n_els):
            print('wrong total number of group elements! (%d == %d)'\
                  % (int(n_el), len(group_n_els)))

        mat_ids = nm.zeros(n_el, dtype=nm.int32)
        for ii, els in enumerate(groups):
            els = nm.array(els, dtype=nm.int32)
            mat_ids[els - 1] = group_ids[ii]

        for elem in el.keys():
            if len(el[elem]) > 0:
                els = nm.array(el[elem], dtype=nm.int32)
                els[:, :-1] -= 1
                els[:, -1] = mat_ids[els[:, -1]-1]

                if elem == '3_8':
                    els = els[:, [0, 1, 3, 2, 4, 5, 7, 6, 8]]

                conns_in.append(els)
                descs.append(elem)

        nod = nm.array(nod, nm.float64)

        conns, mat_ids = split_conns_mat_ids(conns_in)
        mesh._set_io_data(nod, None, conns, mat_ids, descs, nodal_bcs=nodal_bcs)

        return mesh

    def write(self, filename, mesh, out=None, **kwargs):
        raise NotImplementedError

class ANSYSCDBMeshIO(MeshIO):
    format = 'ansys_cdb'

    @staticmethod
    def guess(filename):
        fd = open(filename, 'r')

        for ii in range(1000):
            row = fd.readline()
            if not row: break
            if len(row) == 0: continue

            row = row.split(',')
            kw = row[0].lower()

            if (kw == 'nblock'):
                ok = True
                break

        else:
            ok = False

        fd.close()

        return ok

    @staticmethod
    def make_format(format, nchar=1000):
        idx = []
        dtype = []
        start = 0

        for iform in format:
            ret = iform.partition('i')
            if not ret[1]:
                ret = iform.partition('e')
            if not ret[1]:
                raise ValueError
            aux = ret[2].partition('.')
            step = int(aux[0])
            for j in range(int(ret[0])):
                if (start + step) > nchar:
                    break
                idx.append((start, start+step))
                start += step
                dtype.append(ret[1])

        return idx, dtype

    def write(self, filename, mesh, out=None, **kwargs):
        raise NotImplementedError

    def read_bounding_box(self):
        raise NotImplementedError

    def read_dimension(self, ret_fd=False):
        return 3

    def read(self, mesh, **kwargs):
        ids = []
        coors = []
        tetras = []
        hexas = []
        qtetras = []
        qhexas = []
        nodal_bcs = {}

        fd = open(self.filename, 'r')

        while True:
            row = fd.readline()
            if not row: break
            if len(row) == 0: continue

            row = row.split(',')
            kw = row[0].lower()

            if (kw == 'nblock'):
                # Solid keyword -> 3, otherwise 1 is the starting coors index.
                ic = 3 if len(row) == 3 else 1
                fmt = fd.readline()
                fmt = fmt.strip()[1:-1].split(',')
                row = look_ahead_line(fd)
                nchar = len(row)
                idx, dtype = self.make_format(fmt, nchar)
                ii0, ii1 = idx[0]
                while True:
                    row = fd.readline()
                    if ((row[0] == '!') or (row[:2] == '-1')
                        or len(row) != nchar):
                        break

                    line = [float(row[i0:i1]) for i0, i1 in idx[ic:]]

                    ids.append(int(row[ii0:ii1]))
                    coors.append(line)

            elif (kw == 'eblock'):
                if (len(row) <= 2) or row[2].strip().lower() != 'solid':
                    continue

                fmt = fd.readline()
                fmt = [fmt.strip()[1:-1]]
                row = look_ahead_line(fd)
                nchar = len(row)
                idx, dtype = self.make_format(fmt, nchar)

                imi0, imi1 = idx[0] # Material id.
                inn0, inn1 = idx[8] # Number of nodes in line.
                ien0, ien1 = idx[10] # Element number.
                ic0 = 11
                while True:
                    row = fd.readline()
                    if ((row[0] == '!') or (row[:2] == '-1')
                        or (len(row) != nchar)):
                        break

                    line = [int(row[imi0:imi1])]
                    n_nod = int(row[inn0:inn1])

                    line.extend(int(row[i0:i1])
                                for i0, i1 in idx[ic0 : ic0 + n_nod])
                    if n_nod == 4:
                        tetras.append(line)

                    elif n_nod == 8:
                        hexas.append(line)

                    elif n_nod == 10:
                        row = fd.readline()
                        line.extend(int(row[i0:i1])
                                    for i0, i1 in idx[:2])
                        qtetras.append(line)

                    elif n_nod == 20:
                        row = fd.readline()
                        line.extend(int(row[i0:i1])
                                    for i0, i1 in idx[:12])
                        qhexas.append(line)

                    else:
                        raise ValueError('unsupported element type! (%d nodes)'
                                         % n_nod)

            elif kw == 'cmblock':
                if row[2].lower() != 'node': # Only node sets support.
                    continue

                n_nod = int(row[3].split('!')[0])
                fd.readline() # Format line not needed.

                nods = read_array(fd, n_nod, 1, nm.int32)
                nodal_bcs[row[1].strip()] = nods.ravel()

        fd.close()

        coors = nm.array(coors, dtype=nm.float64)

        tetras = nm.array(tetras, dtype=nm.int32)
        if len(tetras):
            mat_ids_tetras = tetras[:, 0]
            tetras = tetras[:, 1:]

        else:
            tetras.shape = (0, 4)
            mat_ids_tetras = nm.array([])

        hexas = nm.array(hexas, dtype=nm.int32)
        if len(hexas):
            mat_ids_hexas = hexas[:, 0]
            hexas = hexas[:, 1:]

        else:
            hexas.shape = (0, 8)
            mat_ids_hexas = nm.array([])

        if len(qtetras):
            qtetras = nm.array(qtetras, dtype=nm.int32)
            tetras.shape = (max(0, tetras.shape[0]), 4)
            tetras = nm.r_[tetras, qtetras[:, 1:5]]
            mat_ids_tetras = nm.r_[mat_ids_tetras, qtetras[:, 0]]

        if len(qhexas):
            qhexas = nm.array(qhexas, dtype=nm.int32)
            hexas.shape = (max(0, hexas.shape[0]), 8)
            hexas = nm.r_[hexas, qhexas[:, 1:9]]
            mat_ids_hexas = nm.r_[mat_ids_hexas, qhexas[:, 0]]

        if len(qtetras) or len(qhexas):
            ii = nm.union1d(tetras.ravel(), hexas.ravel())
            n_nod = len(ii)

            remap = nm.zeros((ii.max()+1,), dtype=nm.int32)
            remap[ii] = nm.arange(n_nod, dtype=nm.int32)

            ic = nm.searchsorted(ids, ii)
            coors = coors[ic]

        else:
            n_nod = coors.shape[0]
            remap = nm.zeros((nm.array(ids).max() + 1,), dtype=nm.int32)
            remap[ids] = nm.arange(n_nod, dtype=nm.int32)

        # Convert tetras as degenerate hexas to true tetras.
        ii = nm.where((hexas[:, 2] == hexas[:, 3])
                      & (hexas[:, 4] == hexas[:, 5])
                      & (hexas[:, 4] == hexas[:, 6])
                      & (hexas[:, 4] == hexas[:, 7]))[0]

        if len(ii) == len(hexas):
            tetras = nm.r_[tetras, hexas[ii[:, None], [0, 1, 2, 4]]]
            mat_ids_tetras = nm.r_[mat_ids_tetras, mat_ids_hexas[ii]]

            hexas = nm.delete(hexas, ii, axis=0)
            mat_ids_hexas = nm.delete(mat_ids_hexas, ii)

        else:
            output('WARNING: mesh "%s" has both tetrahedra and hexahedra!'
                   % mesh.name)

        ngroups = nm.zeros(len(coors), dtype=nm.int32)

        mesh = mesh_from_groups(mesh, ids, coors, ngroups,
                                [], [], [], [],
                                tetras, mat_ids_tetras,
                                hexas, mat_ids_hexas, remap=remap)

        mesh.nodal_bcs = {}
        for key, nods in six.iteritems(nodal_bcs):
            nods = nods[nods < len(remap)]
            mesh.nodal_bcs[key] = remap[nods]

        return mesh

class XYZMeshIO(MeshIO):
    """
    Trivial XYZ format working only with coordinates (in a .XYZ file) and the
    connectivity stored in another file with the same base name and .IEN
    suffix.
    """
    format = 'xyz'

    def _read_coors(self):
        coors = nm.loadtxt(self.filename)
        if (coors[:, -1] == 0).all():
            coors = coors[:, :-1].copy()

        return coors

    def read_dimension(self, ret_fd=False):
        coors = self._read_coors()
        dim = coors.shape[1]

        if ret_fd:
            fd = open(self.filename, 'r')
            return dim, fd

        else:
            return dim

    def read_bounding_box(self, ret_fd=False, ret_dim=False):
        coors = self._read_coors()
        bbox = nm.vstack((nm.amin(coors, 0), nm.amax(coors, 0)))

        if ret_fd: fd = open(self.filename, 'r')
        if ret_dim:
            dim = coors.shape[1]
            if ret_fd:
                return bbox, dim, fd
            else:
                return bbox, dim
        else:
            if ret_fd:
                return bbox, fd
            else:
                return bbox

    def read(self, mesh, omit_facets=False, **kwargs):
        coors = self._read_coors()
        n_nod, dim = coors.shape

        conn_ext = '.IEN' if op.splitext(self.filename)[1].isupper() else '.ien'
        conn = nm.loadtxt(edit_filename(self.filename, new_ext=conn_ext),
                          dtype=nm.int32) - 1
        desc = '%d_%d' % (dim, conn.shape[1])

        mesh._set_io_data(coors, nm.zeros(n_nod, dtype=nm.int32),
                          [conn], [nm.zeros(conn.shape[0], dtype=nm.int32)],
                          [desc])
        return mesh

    def write(self, filename, mesh, out=None, **kwargs):
        coors, ngroups, conns, mat_ids, desc = mesh._get_io_data()
        n_nod, dim = coors.shape

        zz = nm.zeros((n_nod, 3-dim))
        nm.savetxt(filename, nm.c_[coors, zz])

        conn_ext = '.IEN' if op.splitext(filename)[1].isupper() else '.ien'
        nm.savetxt(edit_filename(self.filename, new_ext=conn_ext),
                   conns[0] + 1)

        if out is not None:
            raise NotImplementedError


var_dict = list(vars().items())
io_table = {}

for key, var in var_dict:
    try:
        if is_derived_class(var, MeshIO):
            io_table[var.format] = var
    except TypeError:
        pass
del var_dict

def any_from_filename(filename, prefix_dir=None, file_format=None, mode='r'):
    """
    Create a MeshIO instance according to the kind of `filename`.

    Parameters
    ----------
    filename : str, function or MeshIO subclass instance
        The name of the mesh file. It can be also a user-supplied function
        accepting two arguments: `mesh`, `mode`, where `mesh` is a Mesh
        instance and `mode` is one of 'read','write', or a MeshIO subclass
        instance.
    prefix_dir : str
        The directory name to prepend to `filename`.

    Returns
    -------
    io : MeshIO subclass instance
        The MeshIO subclass instance corresponding to the kind of `filename`.
    """
    if not isinstance(filename, basestr):
        if isinstance(filename, MeshIO):
            return filename

        else:
            return UserMeshIO(filename)

    if prefix_dir is not None:
        filename = op.normpath(op.join(prefix_dir, filename))

    kwargs = {}
    if file_format is not None:
        if file_format in supported_formats:
            io_class = supported_formats[file_format][0]
            kwargs['file_format'] = file_format
        else:
            raise ValueError('unknown mesh format! (%s)' % file_format)
    else:
        ext2io = {e: (v[0], k) for k, v in supported_formats.items()
            for e in v[1] if '*' not in v[2]}
        ext = op.splitext(filename)[1].lower()
        if ext in ext2io:
            io_class = ext2io[ext][0]
            file_format = ext2io[ext][1]
            kwargs['file_format'] = file_format
        else:
            raise ValueError('unknown mesh format! (%s)' % ext)

    if mode == 'w' and 'w' not in supported_formats[file_format][2]:
        output('writable mesh formats:')
        output_mesh_formats('w')
        msg = 'write support not implemented for output mesh format "%s",' \
              ' see above!' % file_format
        raise ValueError(msg)

    return io_table[io_class](filename, **kwargs)


insert_static_method(MeshIO, any_from_filename)
del any_from_filename

<filename>sfepy/discrete/fem/meshio.py
from __future__ import print_function
from __future__ import absolute_import
import sys
from copy import copy
import logging
import numpy as nm

from sfepy.base.base import (complex_types, dict_from_keys_init,
                             assert_, is_derived_class, ordered_iteritems,
                             insert_static_method, output, get_default,
                             get_default_attr, Struct, basestr)
from sfepy.base.ioutils import (skip_read_line, look_ahead_line, read_token,
                                read_array, pt, enc, dec,
                                edit_filename,
                                read_from_hdf5, write_to_hdf5,
                                HDF5ContextManager, get_or_create_hdf5_group)
import os.path as op
import six
from six.moves import range
import meshio as meshiolib
try:
    from meshio import CellBlock as meshio_Cells  # for meshio >= 4.0.3
except:
    from meshio import Cells as meshio_Cells  # for 4.0.3 > meshio > 4.0.0

_supported_formats = {
    # format name: IO class, suffix, modes[, variants]
    # modes: r = read, w = write, c = test cell groups, v = test vertex groups
    'abaqus': ('meshio', None, 'cv'),
    'exodus': ('meshio', None, 'v'),
    'gmsh': ('meshio', None, 'cv', ['gmsh4-binary', 'gmsh4-ascii',
                                    'gmsh2-binary', 'gmsh2-ascii']),
    'medit': ('meshio', None, 'cv'),
    'nastran': ('meshio', None, 'cv'),
    'vtk': ('meshio', None, 'cv', ['vtk-binnary', 'vtk-ascii']),
    'vtu': ('meshio', None, 'cv'),
    'med': ('meshio', None, 'cv'),
    'xdmf': ('meshio', None, 'cv'),
    'tetgen': ('meshio', None, ''),
    'ansys': ('ansys_cdb', '.cdb', 'r'),
    'hdf5': ('hdf5', '.h5', 'rwcv'),
    'hdf5-xdmf': ('hdf5-xdmf', '.h5x', 'rwcv'),
    'xyz': ('xyz', '.xyz', 'rw'),
    'comsol': ('comsol', '.txt', 'r'),
    'hmascii': ('hmascii', '.hmascii', 'r'),
    'gambit': ('gambit', '.neu', 'r'),
    'mesh3d': ('mesh3d', '.mesh3d', 'r'),
}

def update_supported_formats(formats):
    from meshio._helpers import reader_map, _writer_map,\
        extension_to_filetype

    f2e = {}
    for k, v in extension_to_filetype.items():
        f2e.setdefault(v, []).append(k)

    out = {}
    for format, info in formats.items():
        io, ext, _flag = info[:3]
        variants = info[3] if len(info) >= 4 else []
        for f in [format] + variants:
            if io == 'meshio':
                flag = _flag[:]
                if ext is None:
                    ext = f2e[format]
                if f in _writer_map:
                    flag = 'w' + flag
                if format in reader_map:
                    flag = 'r' + flag
                if not f == format:
                    flag = '*' + flag  # label format variants
            else:
                flag = _flag

            if not isinstance(ext, list):
                ext = [ext]

            out[f] = (io, ext, flag)

    return out

supported_formats = update_supported_formats(_supported_formats)
del _supported_formats


def output_mesh_formats(mode='r'):
    for key, vals in ordered_iteritems(supported_formats):
        if mode in vals[2]:
            output('%s (%s)' % (key,
                vals[1] if len(vals[1]) > 1 else vals[1][0]))


def split_conns_mat_ids(conns_in):
    """
    Split connectivities (columns except the last ones in `conns_in`) from cell
    groups (the last columns of `conns_in`).
    """
    conns, mat_ids = [], []
    for conn in conns_in:
        conn = nm.asarray(conn, dtype=nm.int32)
        conns.append(conn[:, :-1])
        mat_ids.append(conn[:, -1])

    return conns, mat_ids


def convert_complex_output(out_in):
    """
    Convert complex values in the output dictionary `out_in` to pairs of
    real and imaginary parts.
    """
    out = {}
    for key, val in six.iteritems(out_in):

        if val.data.dtype in  complex_types:
            rval = copy(val)
            rval.data = val.data.real
            out['real.%s' % key] = rval

            ival = copy(val)
            ival.data = val.data.imag
            out['imag.%s' % key] = ival

        else:
            out[key] = val

    return out


class MeshIO(Struct):
    """
    The abstract class for importing and exporting meshes.

    Read the docstring of the Mesh() class. Basically all you need to do is to
    implement the read() method::

        def read(self, mesh, **kwargs):
            nodes = ...
            ngroups = ...
            conns = ...
            mat_ids = ...
            descs = ...
            mesh._set_io_data(nodes, ngroups, conns, mat_ids, descs)
            return mesh

    See the Mesh class' docstring how the nodes, ngroups, conns, mat_ids and
    descs should look like. You just need to read them from your specific
    format from disk.

    To write a mesh to disk, just implement the write() method and use the
    information from the mesh instance (e.g. nodes, conns, mat_ids and descs)
    to construct your specific format.

    Optionally, subclasses can implement read_data() to read also computation
    results. This concerns mainly the subclasses with implemented write()
    supporting the 'out' kwarg.

    The default implementation od read_last_step() just returns 0. It should be
    reimplemented in subclasses capable of storing several steps.
    """
    format = None
    call_msg = 'called an abstract MeshIO instance!'

    def __init__(self, filename, **kwargs):
        Struct.__init__(self, filename=filename, **kwargs)
        self.set_float_format()

    def get_filename_trunk(self):
        if isinstance(self.filename, basestr):
            trunk = op.splitext(self.filename)[0]
        else:
            trunk = 'from_descriptor'

        return trunk

    def read_last_step(self):
        """The default implementation: just return 0 as the last step."""
        return 0

    def read_times(self, filename=None):
        """
        Read true time step data from individual time steps.

        Returns
        -------
        steps : array
            The time steps.
        times : array
            The times of the time steps.
        nts : array
            The normalized times of the time steps, in [0, 1].

        Notes
        -----
        The default implementation returns empty arrays.
        """
        aux = nm.array([0.0], dtype=nm.float64)
        return aux.astype(nm.int32), aux, aux

    def read(self, mesh, omit_facets=False, **kwargs):
        raise ValueError(MeshIO.call_msg)

    def write(self, filename, mesh, **kwargs):
        raise ValueError(MeshIO.call_msg)

    def read_data(self, step, filename=None, cache=None):
        raise ValueError(MeshIO.call_msg)

    def set_float_format(self, format=None):
        self.float_format = get_default(format, '%e')

    def get_vector_format(self, dim):
        return ' '.join([self.float_format] * dim)

class UserMeshIO(MeshIO):
    """
    Special MeshIO subclass that enables reading and writing a mesh using a
    user-supplied function.
    """
    format = 'function'

    def __init__(self, filename, **kwargs):
        assert_(hasattr(filename, '__call__'))
        self.function = filename

        MeshIO.__init__(self, filename='function:%s' % self.function.__name__,
                        **kwargs)

    def get_filename_trunk(self):
        return self.filename

    def read(self, mesh, *args, **kwargs):
        aux = self.function(mesh, mode='read')
        if aux is not None:
            mesh = aux

        self.filename = mesh.name

        return mesh

    def write(self, filename, mesh, *args, **kwargs):
        self.function(mesh, mode='write')


def _suppress_meshio_warnings(f):
    def __suppress_meshio_warnings(*args, **kwargs):
        logger = logging.getLogger()
        level = logger.level
        logger.setLevel(logging.ERROR)
        out = f(*args, **kwargs)
        logger.setLevel(level)

        return out

    return __suppress_meshio_warnings

def _decorate_all(module, decorator):
    import types
    for name in dir(module):
        obj = getattr(module, name)
        if isinstance(obj, types.FunctionType):
            setattr(module, name, decorator(obj))

_decorate_all(meshiolib, _suppress_meshio_warnings)
del _decorate_all, _suppress_meshio_warnings


class MeshioLibIO(MeshIO):
    format = 'meshio'

    cell_types = {
        ('hexahedron', 3): '3_8',
        ('tetra', 3): '3_4',
        ('triangle', 3): '2_3',
        ('triangle', 2): '2_3',
        ('quad', 3): '2_4',
        ('quad', 2): '2_4',
        ('line', 3): '3_2',
        ('line', 2): '2_2',
        ('line', 1): '1_2',
    }

    def __init__(self, filename, file_format=None, **kwargs):
        MeshIO.__init__(self, filename=filename, **kwargs)
        from meshio._helpers import _filetype_from_path
        import pathlib

        if file_format is None:
            file_format = _filetype_from_path(pathlib.Path(filename))

        self.file_format = file_format

    def read_bounding_box(self, ret_dim=False):
        m = meshiolib.read(self.filename, file_format=self.file_format)

        bbox = nm.vstack([nm.amin(m.points, 0), nm.amax(m.points, 0)])

        if ret_dim:
            return bbox, m.points.shape[1]
        else:
            return bbox

    @staticmethod
    def _get_dimension(points):
        dim = nm.sum(nm.max(points, axis=0)
                     - nm.min(points, axis=0) > 1e-15)
        return dim

    def read_dimension(self):
        m = meshiolib.read(self.filename, file_format=self.file_format)
        dim = self._get_dimension(m.points)

        return dim

    def read(self, mesh, omit_facets=False, **kwargs):
        m = meshiolib.read(self.filename, file_format=self.file_format)
        dim = self._get_dimension(m.points)

        ngkey = None
        for k in m.point_data.keys():
            if k == 'node_groups' or k.endswith(':ref'):
                ngkey = k
                break

        if ngkey is not None:
            ngroups = nm.asarray(m.point_data[ngkey]).flatten()
        elif hasattr(m, 'point_sets') and len(m.point_sets) > 0:
            ngroups = nm.zeros((len(m.points),), dtype=nm.int32)
            keys = list(m.point_sets.keys())
            keys.sort()
            try:
                ngrps = [int(ii) for ii in keys]
            except:
                ngrps = nm.arange(len(keys)) + 1

            for ik, k in enumerate(keys):
                ngroups[m.point_sets[k]] = ngrps[ik]
        else:
            ngroups = None

        cells, cgroups, cell_types = [], [], []

        # meshio.__version__ > 3.3.2
        cgkey = None
        for k in list(m.cell_data.keys()):
            if k == 'mat_id' or k.endswith(':ref'):
                cgkey = k
                break

        if cgkey is not None:
            cgdata = m.cell_data[cgkey]
        elif len(m.cell_sets) > 0:
            cgdata = []
            keys = list(m.cell_sets.keys())
            keys.sort()
            try:
                cgrps = [int(ii) for ii in keys]
            except:
                cgrps = nm.arange(len(keys)) + 1

            for ic, c in enumerate(m.cells):
                cgdata0 = nm.zeros((len(c.data),), dtype=nm.int32)
                for ik, k in enumerate(keys):
                    cgdata0[m.cell_sets[k][ic]] = cgrps[ik]
                cgdata.append(cgdata0)
        else:
            cgdata = None

        for ic, c in enumerate(m.cells):
            cells.append(c.data)
            cell_types.append(self.cell_types[(c.type, dim)])

            if cgdata is not None:
                cgroups.append(nm.asarray(cgdata[ic]).flatten())
            else:
                cgroups.append(nm.ones((len(c.data),), dtype=nm.int32))

        mesh._set_io_data(m.points[:,:dim], ngroups,
                          cells, cgroups, cell_types)

        output('number of vertices: %d' % m.points.shape[0])
        output('number of cells:')
        for ii, k in enumerate(cell_types):
            output('  %s: %d' % (k, cells[ii].shape[0]))

        return mesh

    def write(self, filename, mesh, out=None, ts=None, **kwargs):
        inv_cell_types = {v: k for k, v in self.cell_types.items()}

        coors, ngroups, conns, _, descs = mesh._get_io_data()

        out = {} if out is None else out

        point_data = {k: v.data for k, v in out.items() if v.mode == 'vertex'}
        cell_data_keys = [k for k, v in out.items() if v.mode == 'cell']

        if self.file_format in ['vtk', 'vtu']:
            ngkey = 'node_groups'
            cgkey = 'mat_id'
        else:
            ngkey = '%s:ref' % self.file_format
            cgkey = '%s:ref' % self.file_format

        point_data[ngkey] = ngroups
        point_sets = {str(k): nm.where(ngroups == k)[0]
            for k in nm.unique(ngroups)}

        cmesh = mesh.cmesh
        cell_groups = cmesh.cell_groups
        cgrps = nm.unique(cell_groups)

        # meshio.__version__ > 3.3.2
        cells = []
        cgroups = [ ]
        cell_data = {k: [] for k in cell_data_keys}
        cell_sets = {str(k): [] for k in cgrps}
        for ii, desc in enumerate(descs):
            cells.append(meshio_Cells(type=inv_cell_types[desc][0],
                                      data=conns[ii]))
            cidxs = nm.where(cmesh.cell_types == cmesh.key_to_index[desc])
            cidxs = cidxs[0].astype(nm.uint32)

            cgroups.append(cell_groups[cidxs])
            for k in cell_data_keys:
                cell_data[k].append(out[k].data[cidxs, 0, :, 0])

            for k in cgrps:
                idxs = nm.where(cell_groups[cidxs] == k)[0]
                cell_sets[str(k)].append(cidxs[idxs])

        cell_data[cgkey] = cgroups

        meshiolib.write_points_cells(filename, coors, cells,
                                     point_data=point_data,
                                     point_sets=point_sets,
                                     cell_data=cell_data,
                                     cell_sets=cell_sets,
                                     file_format=self.file_format)

    def read_data(self, step, filename=None, cache=None):
        m = meshiolib.read(self.filename, file_format=self.file_format)
        dim = self._get_dimension(m.points)

        def _fix_shape(data):
            if data.ndim == 2:
                data = data[:, :dim]

            elif data.ndim == 3:
                data = data[:, None, ...]

            return data

        out = {}
        for key, data in m.point_data.items():
            aux = _fix_shape(data).astype(nm.float64)
            out[key] = Struct(name=key, mode='vertex', data=aux)

        for key, data in m.cell_data.items():
            aux = _fix_shape(data[0]).astype(nm.float64)
            out[key] = Struct(name=key, mode='cell', data=aux)

        return out

class ComsolMeshIO(MeshIO):
    format = 'comsol'

    def _read_commented_int(self):
        return int(skip_read_line(self.fd).split('#')[0])

    def _skip_comment(self):
        read_token(self.fd)
        self.fd.readline()

    def read(self, mesh, **kwargs):

        self.fd = fd = open(self.filename, 'r')
        mode = 'header'

        coors = conns = None
        while 1:
            if mode == 'header':
                line = skip_read_line(fd)

                n_tags = self._read_commented_int()
                for ii in range(n_tags):
                    skip_read_line(fd)
                n_types = self._read_commented_int()
                for ii in range(n_types):
                    skip_read_line(fd)

                skip_read_line(fd)
                assert_(skip_read_line(fd).split()[1] == 'Mesh')
                skip_read_line(fd)
                dim = self._read_commented_int()
                assert_((dim == 2) or (dim == 3))
                n_nod = self._read_commented_int()
                i0 = self._read_commented_int()
                mode = 'points'

            elif mode == 'points':
                self._skip_comment()
                coors = read_array(fd, n_nod, dim, nm.float64)
                mode = 'cells'

            elif mode == 'cells':

                n_types = self._read_commented_int()
                conns = []
                descs = []
                mat_ids = []
                for it in range(n_types):
                    t_name = skip_read_line(fd).split()[1]
                    n_ep = self._read_commented_int()
                    n_el = self._read_commented_int()

                    self._skip_comment()
                    aux = read_array(fd, n_el, n_ep, nm.int32)
                    if t_name == 'tri':
                        conns.append(aux)
                        descs.append('2_3')
                        is_conn = True
                    elif t_name == 'quad':
                        # Rearrange element node order to match SfePy.
                        aux = aux[:,(0,1,3,2)]
                        conns.append(aux)
                        descs.append('2_4')
                        is_conn = True
                    elif t_name == 'hex':
                        # Rearrange element node order to match SfePy.
                        aux = aux[:,(0,1,3,2,4,5,7,6)]
                        conns.append(aux)
                        descs.append('3_8')
                        is_conn = True
                    elif t_name == 'tet':
                        conns.append(aux)
                        descs.append('3_4')
                        is_conn = True
                    else:
                        is_conn = False

                    # Skip parameters.
                    n_pv = self._read_commented_int()
                    n_par = self._read_commented_int()
                    for ii in range(n_par):
                        skip_read_line(fd)

                    n_domain = self._read_commented_int()
                    assert_(n_domain == n_el)
                    if is_conn:
                        self._skip_comment()
                        mat_id = read_array(fd, n_domain, 1, nm.int32)
                        mat_ids.append(mat_id.squeeze())
                    else:
                        for ii in range(n_domain):
                            skip_read_line(fd)

                    # Skip up/down pairs.
                    n_ud = self._read_commented_int()
                    for ii in range(n_ud):
                        skip_read_line(fd)
                break

        fd.close()
        self.fd = None

        mesh._set_io_data(coors, None, conns, mat_ids, descs)

        return mesh

    def write(self, filename, mesh, out=None, **kwargs):

        def write_elements(fd, ig, conn, mat_ids, type_name,
                           npe, format, norder, nm_params):
            fd.write("# Type #%d\n\n" % ig)
            fd.write("%s # type name\n\n\n" % type_name)
            fd.write("%d # number of nodes per element\n" % npe)
            fd.write("%d # number of elements\n" % conn.shape[0])
            fd.write("# Elements\n")
            for ii in range(conn.shape[0]):
                nn = conn[ii] # Zero based
                fd.write(format % tuple(nn[norder]))
            fd.write("\n%d # number of parameter values per element\n"
                     % nm_params)
            # Top level always 0?
            fd.write("0 # number of parameters\n")
            fd.write("# Parameters\n\n")
            fd.write("%d # number of domains\n"
                     % sum([mi.shape[0] for mi in mat_ids]))
            fd.write("# Domains\n")
            for mi in mat_ids:
                # Domains in comsol have to be > 0
                if (mi <= 0).any():
                    mi += mi.min() + 1
                for dom in mi:
                    fd.write("%d\n" % abs(dom))
            fd.write("\n0 # number of up/down pairs\n")
            fd.write("# Up/down\n")

        fd = open(filename, 'w')

        coors, ngroups, conns, mat_ids, desc = mesh._get_io_data()

        n_nod, dim = coors.shape

        # Header
        fd.write("# Created by SfePy\n\n\n")
        fd.write("# Major & minor version\n")
        fd.write("0 1\n")
        fd.write("1 # number of tags\n")
        fd.write("# Tags\n")
        fd.write("2 m1\n")
        fd.write("1 # number of types\n")
        fd.write("# Types\n")
        fd.write("3 obj\n\n")

        # Record
        fd.write("# --------- Object 0 ----------\n\n")
        fd.write("0 0 1\n") # version unused serializable
        fd.write("4 Mesh # class\n")
        fd.write("1 # version\n")
        fd.write("%d # sdim\n" % dim)
        fd.write("%d # number of mesh points\n" % n_nod)
        fd.write("0 # lowest mesh point index\n\n") # Always zero in SfePy

        fd.write("# Mesh point coordinates\n")

        format = self.get_vector_format(dim) + '\n'
        for ii in range(n_nod):
            nn = tuple(coors[ii])
            fd.write(format % tuple(nn))

        fd.write("\n%d # number of element types\n\n\n" % len(conns))

        for ig, conn in enumerate(conns):
            if (desc[ig] == "2_4"):
                write_elements(fd, ig, conn, mat_ids,
                               "4 quad", 4, "%d %d %d %d\n", [0, 1, 3, 2], 8)

            elif (desc[ig] == "2_3"):
                # TODO: Verify number of parameters for tri element
                write_elements(fd, ig, conn, mat_ids,
                               "3 tri", 3, "%d %d %d\n", [0, 1, 2], 4)

            elif (desc[ig] == "3_4"):
                # TODO: Verify number of parameters for tet element
                write_elements(fd, ig, conn, mat_ids,
                               "3 tet", 4, "%d %d %d %d\n", [0, 1, 2, 3], 16)

            elif (desc[ig] == "3_8"):
                write_elements(fd, ig, conn, mat_ids,
                               "3 hex", 8, "%d %d %d %d %d %d %d %d\n",
                               [0, 1, 3, 2, 4, 5, 7, 6], 24)

            else:
                raise ValueError('unknown element type! (%s)' % desc[ig])

        fd.close()

        if out is not None:
            for key, val in six.iteritems(out):
                raise NotImplementedError


class HDF5MeshIO(MeshIO):
    format = "hdf5"

    import string
    _all = ''.join(map(chr, list(range(256))))
    _letters = string.ascii_letters + string.digits + '_'
    _rubbish = ''.join([ch for ch in set(_all) - set(_letters)])
    if sys.version_info[0] >= 3:
        _tr = str.maketrans(_rubbish, '_' * len(_rubbish))
    else:
        _tr = string.maketrans(_rubbish, '_' * len(_rubbish))

    @staticmethod
    def read_mesh_from_hdf5(filename, group=None, mesh=None):
        """
        Read the mesh from a HDF5 file.

        filename: str or tables.File
            The HDF5 file to read the mesh from.
        group: tables.group.Group or str, optional
            The HDF5 file group to read the mesh from.
            If None, the root group is used.
        mesh: sfepy.dicrete.fem.Mesh or None
            If None, the new mesh is created and returned, otherwise
            content of this argument is replaced by the read mesh.

        Returns
        -------
        sfepy.dicrete.fem.Mesh
            readed mesh
        """
        with HDF5ContextManager(filename, mode='r') as fd:
            if group is None:
                group = fd.root
            elif not isinstance(group, pt.group.Group):
                group = fd.get_node(group)

            set_shape_info = mesh is None
            if mesh is None:
                from .mesh import Mesh
                mesh = Mesh('mesh')

            mesh.name = dec(group.name.read())
            coors = group.coors.read()
            ngroups = group.ngroups.read()

            n_gr = group.n_gr.read()

            conns = []
            descs = []
            mat_ids = []
            for ig in range(n_gr):
                gr_name = 'group%d' % ig
                conn_group = group._f_get_child(gr_name)
                conns.append(conn_group.conn.read())
                mat_ids.append(conn_group.mat_id.read())
                descs.append(dec(conn_group.desc.read()))

            nodal_bcs = {}
            try:
                node_sets_groups = group.node_sets

            except:
                pass

            else:
                for group in node_sets_groups:
                    key = dec(group.key.read())
                    nods = group.nods.read()
                    nodal_bcs[key] = nods

            mesh._set_io_data(coors, ngroups, conns, mat_ids, descs,
                              nodal_bcs=nodal_bcs)

            if set_shape_info:
                mesh._set_shape_info()
        return mesh

    @staticmethod
    def write_mesh_to_hdf5(filename, group, mesh, force_3d=False):
        """
        Write mesh to a hdf5 file.

        filename: str or tables.File
            The HDF5 file to write the mesh to.
        group: tables.group.Group or None or str
            The HDF5 file group to write the mesh to.
            If None, the root group is used.
            The group can be given as a path from root, e.g. /path/to/mesh
        mesh: sfepy.dicrete.fem.Mesh
            The mesh to write.
        """
        with HDF5ContextManager(filename, mode='w') as fd:
            if group is None:
                group = fd.root
            elif not isinstance(group, pt.group.Group):
                group = get_or_create_hdf5_group(fd, group)

            coors, ngroups, conns, mat_ids, descs = mesh._get_io_data()
            if force_3d and coors.shape[1] == 2:
                coors = nm.hstack([coors, nm.zeros((len(coors), 1))])
            fd.create_array(group, 'name', enc(mesh.name), 'name')
            fd.create_array(group, 'coors', coors, 'coors')
            fd.create_array(group, 'ngroups', ngroups, 'ngroups')
            fd.create_array(group, 'n_gr', len(conns), 'n_gr')
            for ig, conn in enumerate(conns):
                conn_group = fd.create_group(group, 'group%d' % ig,
                                             'connectivity group')
                fd.create_array(conn_group, 'conn', conn, 'connectivity')
                fd.create_array(conn_group, 'mat_id', mat_ids[ig],
                                'material id')
                fd.create_array(conn_group, 'desc', enc(descs[ig]),
                                'element Type')

            node_sets_groups = fd.create_group(group, 'node_sets',
                                               'node sets groups')
            ii = 0
            for key, nods in six.iteritems(mesh.nodal_bcs):
                group = fd.create_group(node_sets_groups, 'group%d' % ii,
                                        'node sets group')
                fd.create_array(group, 'key', enc(key), 'key')
                fd.create_array(group, 'nods', nods, 'nods')
                ii += 1

    def read_dimension(self, ret_fd=False):
        fd = pt.open_file(self.filename, mode="r")

        dim = fd.root.mesh.coors.shape[1]

        if ret_fd:
            return dim, fd

        else:
            fd.close()
            return dim

    def read_bounding_box(self, ret_fd=False, ret_dim=False):
        fd = pt.open_file(self.filename, mode="r")

        mesh_group = fd.root.mesh

        coors = mesh_group.coors.read()
        bbox = nm.vstack((nm.amin(coors, 0),
                          nm.amax(coors, 0)))

        if ret_dim:
            dim = coors.shape[1]

            if ret_fd:
                return bbox, dim, fd

            else:
                fd.close()
                return bbox, dim

        else:
            if ret_fd:
                return bbox, fd

            else:
                fd.close()
                return bbox

    def read(self, mesh=None, **kwargs):
        return self.read_mesh_from_hdf5(self.filename, '/mesh', mesh=mesh)

    @staticmethod
    def write_xdmf_file(filename, **kwargs):
        def get_path(node):
            _, fname = op.split(filename)
            return '%s:%s' % (fname, node._v_pathname)

        def get_data_dim(shape):
            if len(shape) == 4:
               return shape[2]
            if len(shape) == 2:
               return shape[1]
            else:
               return 1

        def data_item(data):
            dtype = data.dtype
            if nm.issubdtype(dtype, nm.integer):
                data_type = 'Int'
            elif nm.issubdtype(dtype, nm.floating):
                data_type = 'Float'
            else:
                raise ValueError('wrong data type! (%s)' % dtype)

            dim = get_data_dim(data.shape)
            sh = (data.shape[0], dim)
            ditem = et.Element('DataItem',
                               attrib={'DataType': data_type,
                                       'Dimensions': '%d %d' % sh,
                                       'Format': 'HDF',
                                       'Precision': str(dtype.itemsize)})
            ditem.text = get_path(data)

            return ditem

        def attr_item(data, center=None, name=None):
            if isinstance(data, pt.Group):
                if center is None and 'mode' in data:
                    mode = data.mode.read().decode('ascii')
                    if mode == 'custom':
                        return None
                    center = {'vertex': 'Node',
                              'cell': 'Cell'}[mode]
                if name is None and 'dname' in data:
                    name = data.dname.read().decode('ascii').lstrip('_')
                data = data.data

            dim = get_data_dim(data.shape)
            atype = {1: 'Scalar', 3: 'Vector', 6: 'Tensor6', 9: 'Tensor'}[dim]

            aitem = et.Element('Attribute',
                               attrib={'AttributeType': atype,
                                       'Center': center,
                                       'Name': name})
            aitem.append(data_item(data))

            return aitem

        import xml.etree.ElementTree as et
        from xml.dom import minidom

        topology_table = {
            '2_2': 'Line',
            '3_2': 'Line',
            '2_3': 'Triangle',
            '2_4': 'Quadrilateral',
            '3_4': 'Tetrahedron',
            '3_8': 'Hexahedron',
        }

        et_root = et.Element('Xdmf', attrib={'Version': '3.0'})
        if 'extra_data' in kwargs:
            for k, v in kwargs['extra_data'].items():
                d = et.SubElement(et_root, 'Information', attrib={'Name': k})
                d.text = str(v)

        et_domain = et.SubElement(et_root, 'Domain')

        with HDF5ContextManager(filename, mode='r') as fd:
            root = fd.root
            mesh = root.mesh
            name = mesh.name.read().decode('ascii')
            et_mesh = []
            geom = et.Element('Geometry', attrib={'GeometryType': 'XYZ'})
            geom.append(data_item(mesh.coors))
            et_mesh.append(geom)
            et_mesh.append(attr_item(mesh.ngroups, 'Node', 'node_groups'))

            n_gr = mesh.n_gr.read()
            for ig in range(n_gr):
                gr_name = 'group%d' % ig
                conn_group = mesh._f_get_child(gr_name)
                nc, nnd = conn_group.conn.shape
                ttype = topology_table[dec(conn_group.desc.read())]
                et_conn = et.Element('Topology',
                                     attrib={'NumberOfElements': str(nc),
                                             'TopologyType': ttype})
                et_conn.append(data_item(conn_group.conn))
                et_mesh.append(et_conn)
                et_mesh.append(attr_item(conn_group.mat_id, 'Cell', 'mat_id'))

            steps = [k for k in root if k._v_name.startswith('step')]
            et_ts = et.SubElement(et_domain, 'Grid',
                                    attrib={'Name': 'TimeSeries',
                                            'GridType': 'Collection',
                                            'CollectionType': 'Temporal'})

            for step in steps:
                istep = int(step._v_name[4:])
                et_grid = et.SubElement(et_ts, 'Grid',
                                        attrib={'Name': 'grid%d' % istep,
                                                'GridType': 'Uniform'})
                et.SubElement(et_grid, 'Time', attrib={'Value': '%d' % istep})
                et_grid.extend(et_mesh)

                for val in filter(lambda x: x._v_name.startswith('__'), step):
                    aitem = attr_item(val)
                    if aitem is not None:
                        et_grid.append(aitem)

        out = minidom.parseString(et.tostring(et_root)).toprettyxml(indent="  ")
        xdmf_filename = op.splitext(filename)[0] + '.xdmf'
        with open(xdmf_filename, 'w') as f:
            f.write(out[(out.find('\n') + 1):])

    def write(self, filename, mesh, out=None, ts=None, cache=None,
              xdmf=False, **kwargs):
        def expand_data_3d(data):
            expand_tab = {
                2: (3, [0, 1]),
                3: (6, [0, 1, 3]),
                4: (9, [0, 1, 3, 4]),
            }

            n, dim = data.shape
            if dim in expand_tab:
                dim3, order = expand_tab[dim]
                out = nm.zeros((n, dim3), dtype=data.dtype)
                out[:, order] = data
                return out
            else:
                return data

        from time import asctime

        if pt is None:
            raise ValueError('pytables not imported!')

        step = get_default_attr(ts, 'step', 0)
        if (step == 0) or not op.exists(filename):
            # A new file.
            with pt.open_file(filename, mode="w",
                              title="SfePy output file") as fd:
                mesh_group = fd.create_group('/', 'mesh', 'mesh')
                self.write_mesh_to_hdf5(fd, mesh_group, mesh, force_3d=xdmf)

                if ts is not None:
                    ts_group = fd.create_group('/', 'ts', 'time stepper')
                    fd.create_array(ts_group, 't0', ts.t0, 'initial time')
                    fd.create_array(ts_group, 't1', ts.t1, 'final time' )
                    fd.create_array(ts_group, 'dt', ts.dt, 'time step')
                    fd.create_array(ts_group, 'n_step', ts.n_step, 'n_step')

                tstat_group = fd.create_group('/', 'tstat',
                                              'global time statistics')
                fd.create_array(tstat_group, 'created', enc(asctime()),
                                'file creation time')
                fd.create_array(tstat_group, 'finished', enc('.' * 24),
                                'file closing time')

                fd.create_array(fd.root, 'last_step',
                                nm.array([step], dtype=nm.int32),
                                'last saved step')

        if out is not None:
            if ts is None:
                step, time, nt  = 0, 0.0, 0.0
            else:
                step, time, nt = ts.step, ts.time, ts.nt

            # Existing file.
            fd = pt.open_file(filename, mode="r+")

            step_group_name = 'step%d' % step
            if step_group_name in fd.root:
                raise ValueError('step %d is already saved in "%s" file!'
                                 ' Possible help: remove the old file or'
                                 ' start saving from the initial time.'
                                 % (step, filename))
            step_group = fd.create_group('/', step_group_name, 'time step data')

            ts_group = fd.create_group(step_group, 'ts', 'time stepper')
            fd.create_array(ts_group, 'step', step, 'step')
            fd.create_array(ts_group, 't', time, 'time')
            fd.create_array(ts_group, 'nt', nt, 'normalized time')

            name_dict = {}
            for key, val in six.iteritems(out):
                if xdmf and mesh.coors.shape[1] == 2:
                    data = expand_data_3d(val.data)
                else:
                    data = val.data

                group_name = '__' + key.translate(self._tr)
                data_group = fd.create_group(step_group, group_name,
                                             '%s data' % key)
                fd.create_array(data_group, 'dname', enc(key), 'data name')
                fd.create_array(data_group, 'mode', enc(val.mode), 'mode')
                name = val.get('name', 'output_data')
                fd.create_array(data_group, 'name', enc(name), 'object name')
                if val.mode == 'custom':
                    write_to_hdf5(fd, data_group, 'data', data,
                                  cache=cache,
                                  unpack_markers=getattr(val, 'unpack_markers',
                                                         False))
                    continue

                shape = val.get('shape', data.shape)
                dofs = val.get('dofs', None)
                if dofs is None:
                    dofs = [''] * nm.squeeze(shape)[-1]
                var_name = val.get('var_name', '')

                fd.create_array(data_group, 'data', data, 'data')
                fd.create_array(data_group, 'dofs', [enc(ic) for ic in dofs],
                                'dofs')
                fd.create_array(data_group, 'shape', shape, 'shape')
                fd.create_array(data_group, 'var_name',
                                enc(var_name), 'object parent name')
                if val.mode == 'full':
                    fd.create_array(data_group, 'field_name',
                                    enc(val.field_name), 'field name')

                name_dict[key] = group_name

            step_group._v_attrs.name_dict = name_dict
            fd.root.last_step[0] = step

            fd.remove_node(fd.root.tstat.finished)
            fd.create_array(fd.root.tstat, 'finished', enc(asctime()),
                            'file closing time')
            fd.close()

        if xdmf:
            self.write_xdmf_file(filename, **kwargs)

    def read_last_step(self, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")
        last_step = fd.root.last_step[0]
        fd.close()
        return last_step

    def read_time_stepper(self, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        try:
            ts_group = fd.root.ts
            out =  (ts_group.t0.read(), ts_group.t1.read(),
                    ts_group.dt.read(), ts_group.n_step.read())

        except:
            raise ValueError('no time stepper found!')

        finally:
            fd.close()

        return out

    def _get_step_group_names(self, fd):
        return sorted([name for name in fd.root._v_groups.keys()
                       if name.startswith('step')],
                      key=lambda name: int(name[4:]))

    def read_times(self, filename=None):
        """
        Read true time step data from individual time steps.

        Returns
        -------
        steps : array
            The time steps.
        times : array
            The times of the time steps.
        nts : array
            The normalized times of the time steps, in [0, 1].
        """
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode='r')

        steps = []
        times = []
        nts = []
        for gr_name in self._get_step_group_names(fd):
            ts_group = fd.get_node(fd.root, gr_name + '/ts')

            steps.append(ts_group.step.read())
            times.append(ts_group.t.read())
            nts.append(ts_group.nt.read())
        fd.close()

        steps = nm.asarray(steps, dtype=nm.int32)
        times = nm.asarray(times, dtype=nm.float64)
        nts = nm.asarray(nts, dtype=nm.float64)

        return steps, times, nts

    def _get_step_group(self, step, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        if step is None:
            step = int(self._get_step_group_names(fd)[0][4:])

        gr_name = 'step%d' % step
        try:
            step_group = fd.get_node(fd.root, gr_name)
        except:
            output('step %d data not found - premature end of file?' % step)
            fd.close()
            return None, None

        return fd, step_group

    def read_data(self, step, filename=None, cache=None):
        fd, step_group = self._get_step_group(step, filename=filename)
        if fd is None: return None

        out = {}
        for data_group in step_group:
            try:
                key = dec(data_group.dname.read())

            except pt.exceptions.NoSuchNodeError:
                continue

            mode = dec(data_group.mode.read())
            if mode == 'custom':
               out[key] = read_from_hdf5(fd, data_group.data, cache=cache)
               continue

            name = dec(data_group.name.read())
            data = data_group.data.read()
            dofs = tuple([dec(ic) for ic in data_group.dofs.read()])
            try:
                shape = tuple(int(ii) for ii in data_group.shape.read())

            except pt.exceptions.NoSuchNodeError:
                shape = data.shape

            if mode == 'full':
                field_name = dec(data_group.field_name.read())

            else:
                field_name = None

            out[key] = Struct(name=name, mode=mode, data=data,
                              dofs=dofs, shape=shape, field_name=field_name)

            if out[key].dofs == (-1,):
                out[key].dofs = None

        fd.close()

        return out

    def read_data_header(self, dname, step=None, filename=None):
        fd, step_group = self._get_step_group(step, filename=filename)
        if fd is None: return None

        groups = step_group._v_groups
        for name, data_group in six.iteritems(groups):
            try:
                key = dec(data_group.dname.read())

            except pt.exceptions.NoSuchNodeError:
                continue

            if key == dname:
                mode = dec(data_group.mode.read())
                fd.close()
                return mode, name

        fd.close()
        raise KeyError('non-existent data: %s' % dname)

    def read_time_history(self, node_name, indx, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        th = dict_from_keys_init(indx, list)
        for gr_name in self._get_step_group_names(fd):
            step_group = fd.get_node(fd.root, gr_name)
            data = step_group._f_get_child(node_name).data

            for ii in indx:
                th[ii].append(nm.array(data[ii]))

        fd.close()

        for key, val in six.iteritems(th):
            aux = nm.array(val)
            if aux.ndim == 4: # cell data.
                aux = aux[:,0,:,0]
            th[key] = aux

        return th

    def read_variables_time_history(self, var_names, ts, filename=None):
        filename = get_default(filename, self.filename)
        fd = pt.open_file(filename, mode="r")

        assert_((fd.root.last_step[0] + 1) == ts.n_step)

        ths = dict_from_keys_init(var_names, list)

        arr = nm.asarray
        for step in range(ts.n_step):
            gr_name = 'step%d' % step
            step_group = fd.get_node(fd.root, gr_name)
            name_dict = step_group._v_attrs.name_dict
            for var_name in var_names:
                data = step_group._f_get_child(name_dict[var_name]).data
                ths[var_name].append(arr(data.read()))

        fd.close()

        return ths


class HDF5XdmfMeshIO(HDF5MeshIO):
    format="hdf5-xdmf"

    def write(self, filename, mesh, out=None, ts=None, cache=None, **kwargs):
        HDF5MeshIO.write(self, filename, mesh, out=out, ts=ts, cache=cache,
                         xdmf=True, **kwargs)


class Mesh3DMeshIO(MeshIO):
    format = "mesh3d"

    def read(self, mesh, **kwargs):
        f = open(self.filename)
        # read the whole file:
        vertices = self._read_section(f, integer=False)
        tetras = self._read_section(f)
        hexes = self._read_section(f)
        prisms = self._read_section(f)
        tris = self._read_section(f)
        quads = self._read_section(f)

        # substract 1 from all elements, because we count from 0:
        conns = []
        mat_ids = []
        descs = []
        if len(tetras) > 0:
            conns.append(tetras - 1)
            mat_ids.append([0]*len(tetras))
            descs.append("3_4")
        if len(hexes) > 0:
            conns.append(hexes - 1)
            mat_ids.append([0]*len(hexes))
            descs.append("3_8")
        mesh._set_io_data(vertices, None, conns, mat_ids, descs)
        return mesh

    def read_dimension(self):
        return 3

    def _read_line(self, f):
        """
        Reads one non empty line (if it's a comment, it skips it).
        """
        l = f.readline().strip()
        while l == "" or l[0] == "#": # comment or an empty line
            l = f.readline().strip()
        return l

    def _read_section(self, f, integer=True):
        """
        Reads one section from the mesh3d file.

        integer ... if True, all numbers are passed to int(), otherwise to
            float(), before returning

        Some examples how a section can look like:

        2
        1 2 5 4 7 8 11 10
        2 3 6 5 8 9 12 11

        or

        5
        1 2 3 4     1
        1 2 6 5     1
        2 3 7 6     1
        3 4 8 7     1
        4 1 5 8     1

        or

        0

        """
        if integer:
            dtype=int
        else:
            dtype=float
        l = self._read_line(f)
        N = int(l)
        rows = []
        for i in range(N):
            l = self._read_line(f)
            row = nm.fromstring(l, sep=" ", dtype=dtype)
            rows.append(row)
        return nm.array(rows)

def mesh_from_groups(mesh, ids, coors, ngroups,
                     tris, mat_tris, quads, mat_quads,
                     tetras, mat_tetras, hexas, mat_hexas, remap=None):
    ids = nm.asarray(ids, dtype=nm.int32)
    coors = nm.asarray(coors, dtype=nm.float64)

    if remap is None:
        n_nod = coors.shape[0]
        remap = nm.zeros((ids.max()+1,), dtype=nm.int32)
        remap[ids] = nm.arange(n_nod, dtype=nm.int32)

    tris = remap[nm.array(tris, dtype=nm.int32)]
    quads = remap[nm.array(quads, dtype=nm.int32)]
    tetras = remap[nm.array(tetras, dtype=nm.int32)]
    hexas = remap[nm.array(hexas, dtype=nm.int32)]

    conns = [tris, quads, tetras, hexas]
    mat_ids = [nm.array(ar, dtype=nm.int32)
               for ar in [mat_tris, mat_quads, mat_tetras, mat_hexas]]
    descs = ['2_3', '2_4', '3_4', '3_8']

    # Remove empty groups.
    conns, mat_ids, descs = zip(*[(conns[ig], mat_ids[ig], descs[ig])
                                  for ig in range(4)
                                  if conns[ig].shape[0] > 0])

    mesh._set_io_data(coors, ngroups, conns, mat_ids, descs)
    return mesh


class HypermeshAsciiMeshIO(MeshIO):
    format = 'hmascii'

    def read(self, mesh, **kwargs):
        fd = open(self.filename, 'r')

        ids = []
        coors = []
        tetras = []
        mat_tetras = []
        hexas = []
        mat_hexas = []
        quads = []
        mat_quads = []
        trias = []
        mat_trias = []
        mat_id = 0

        for line in fd:
            if line and (line[0] == '*'):
                if line[1:10] == 'component':
                    line = line.strip()[11:-1].split(',')
                    mat_id = int(line[0])
                if line[1:5] == 'node':
                    line = line.strip()[6:-1].split(',')
                    ids.append(int(line[0]))
                    coors.append([float(coor) for coor in line[1:4]])

                elif line[1:7] == 'tetra4':
                    line = line.strip()[8:-1].split(',')
                    mat_tetras.append(mat_id)
                    tetras.append([int(ic) for ic in line[2:6]])

                elif line[1:6] == 'hexa8':
                    line = line.strip()[7:-1].split(',')
                    mat_hexas.append(mat_id)
                    hexas.append([int(ic) for ic in line[2:10]])

                elif line[1:6] == 'quad4':
                    line = line.strip()[7:-1].split(',')
                    mat_quads.append(mat_id)
                    quads.append([int(ic) for ic in line[2:6]])

                elif line[1:6] == 'tria3':
                    line = line.strip()[7:-1].split(',')
                    mat_trias.append(mat_id)
                    trias.append([int(ic) for ic in line[2:5]])
        fd.close()

        mesh = mesh_from_groups(mesh, ids, coors, None,
                                trias, mat_trias, quads, mat_quads,
                                tetras, mat_tetras, hexas, mat_hexas)

        return mesh

    def read_dimension(self):
        return 3

    def write(self, filename, mesh, out=None, **kwargs):
        raise NotImplementedError


class NEUMeshIO(MeshIO):
    format = 'gambit'

    def read_dimension(self, ret_fd=False):

        fd = open(self.filename, 'r')

        row = fd.readline().split()
        while 1:
            if not row: break
            if len(row) == 0: continue

            if (row[0] == 'NUMNP'):
                row = fd.readline().split()
                n_nod, n_el, dim = row[0], row[1], int(row[4])
                break

        if ret_fd:
            return dim, fd
        else:
            fd.close()
            return dim

    def read(self, mesh, **kwargs):

        el = {'3_8' : [], '3_4' : [], '2_4' : [], '2_3' : []}
        nod = []

        conns_in = []
        descs = []

        group_ids = []
        group_n_els = []
        groups = []
        nodal_bcs = {}

        fd = open(self.filename, 'r')

        row = fd.readline()
        while 1:
            if not row: break
            row = row.split()

            if len(row) == 0:
                row = fd.readline()
                continue

            if (row[0] == 'NUMNP'):
                row = fd.readline().split()
                n_nod, n_el, dim = int(row[0]), int(row[1]), int(row[4])

            elif (row[0] == 'NODAL'):
                row = fd.readline().split()
                while not(row[0] == 'ENDOFSECTION'):
                    nod.append(row[1:])
                    row = fd.readline().split()

            elif (row[0] == 'ELEMENTS/CELLS'):
                row = fd.readline().split()
                while not(row[0] == 'ENDOFSECTION'):
                    elid = [row[0]]
                    gtype = int(row[1])
                    if gtype == 6:
                        el['3_4'].append(row[3:]+elid)
                    elif gtype == 4:
                        rr = row[3:]
                        if (len(rr) < 8):
                            rr.extend(fd.readline().split())
                        el['3_8'].append(rr+elid)
                    elif gtype == 3:
                        el['2_3'].append(row[3:]+elid)
                    elif gtype == 2:
                        el['2_4'].append(row[3:]+elid)
                    row = fd.readline().split()

            elif (row[0] == 'GROUP:'):
                group_ids.append(row[1])
                g_n_el = int(row[3])
                group_n_els.append(g_n_el)
                name = fd.readline().strip()

                els = []
                row = fd.readline().split()
                row = fd.readline().split()
                while not(row[0] == 'ENDOFSECTION'):
                    els.extend(row)
                    row = fd.readline().split()
                if g_n_el != len(els):
                    msg = 'wrong number of group elements! (%d == %d)'\
                          % (n_el, len(els))
                    raise ValueError(msg)
                groups.append(els)

            elif (row[0] == 'BOUNDARY'):
                row = fd.readline().split()
                key = row[0]
                num = int(row[2])
                inod = read_array(fd, num, None, nm.int32) - 1
                nodal_bcs[key] = inod.squeeze()

                row = fd.readline().split()
                assert_(row[0] == 'ENDOFSECTION')

            row = fd.readline()

        fd.close()

        if int(n_el) != sum(group_n_els):
            print('wrong total number of group elements! (%d == %d)'\
                  % (int(n_el), len(group_n_els)))

        mat_ids = nm.zeros(n_el, dtype=nm.int32)
        for ii, els in enumerate(groups):
            els = nm.array(els, dtype=nm.int32)
            mat_ids[els - 1] = group_ids[ii]

        for elem in el.keys():
            if len(el[elem]) > 0:
                els = nm.array(el[elem], dtype=nm.int32)
                els[:, :-1] -= 1
                els[:, -1] = mat_ids[els[:, -1]-1]

                if elem == '3_8':
                    els = els[:, [0, 1, 3, 2, 4, 5, 7, 6, 8]]

                conns_in.append(els)
                descs.append(elem)

        nod = nm.array(nod, nm.float64)

        conns, mat_ids = split_conns_mat_ids(conns_in)
        mesh._set_io_data(nod, None, conns, mat_ids, descs, nodal_bcs=nodal_bcs)

        return mesh

    def write(self, filename, mesh, out=None, **kwargs):
        raise NotImplementedError

class ANSYSCDBMeshIO(MeshIO):
    format = 'ansys_cdb'

    @staticmethod
    def guess(filename):
        fd = open(filename, 'r')

        for ii in range(1000):
            row = fd.readline()
            if not row: break
            if len(row) == 0: continue

            row = row.split(',')
            kw = row[0].lower()

            if (kw == 'nblock'):
                ok = True
                break

        else:
            ok = False

        fd.close()

        return ok

    @staticmethod
    def make_format(format, nchar=1000):
        idx = []
        dtype = []
        start = 0

        for iform in format:
            ret = iform.partition('i')
            if not ret[1]:
                ret = iform.partition('e')
            if not ret[1]:
                raise ValueError
            aux = ret[2].partition('.')
            step = int(aux[0])
            for j in range(int(ret[0])):
                if (start + step) > nchar:
                    break
                idx.append((start, start+step))
                start += step
                dtype.append(ret[1])

        return idx, dtype

    def write(self, filename, mesh, out=None, **kwargs):
        raise NotImplementedError

    def read_bounding_box(self):
        raise NotImplementedError

    def read_dimension(self, ret_fd=False):
        return 3

    def read(self, mesh, **kwargs):
        ids = []
        coors = []
        tetras = []
        hexas = []
        qtetras = []
        qhexas = []
        nodal_bcs = {}

        fd = open(self.filename, 'r')

        while True:
            row = fd.readline()
            if not row: break
            if len(row) == 0: continue

            row = row.split(',')
            kw = row[0].lower()

            if (kw == 'nblock'):
                # Solid keyword -> 3, otherwise 1 is the starting coors index.
                ic = 3 if len(row) == 3 else 1
                fmt = fd.readline()
                fmt = fmt.strip()[1:-1].split(',')
                row = look_ahead_line(fd)
                nchar = len(row)
                idx, dtype = self.make_format(fmt, nchar)
                ii0, ii1 = idx[0]
                while True:
                    row = fd.readline()
                    if ((row[0] == '!') or (row[:2] == '-1')
                        or len(row) != nchar):
                        break

                    line = [float(row[i0:i1]) for i0, i1 in idx[ic:]]

                    ids.append(int(row[ii0:ii1]))
                    coors.append(line)

            elif (kw == 'eblock'):
                if (len(row) <= 2) or row[2].strip().lower() != 'solid':
                    continue

                fmt = fd.readline()
                fmt = [fmt.strip()[1:-1]]
                row = look_ahead_line(fd)
                nchar = len(row)
                idx, dtype = self.make_format(fmt, nchar)

                imi0, imi1 = idx[0] # Material id.
                inn0, inn1 = idx[8] # Number of nodes in line.
                ien0, ien1 = idx[10] # Element number.
                ic0 = 11
                while True:
                    row = fd.readline()
                    if ((row[0] == '!') or (row[:2] == '-1')
                        or (len(row) != nchar)):
                        break

                    line = [int(row[imi0:imi1])]
                    n_nod = int(row[inn0:inn1])

                    line.extend(int(row[i0:i1])
                                for i0, i1 in idx[ic0 : ic0 + n_nod])
                    if n_nod == 4:
                        tetras.append(line)

                    elif n_nod == 8:
                        hexas.append(line)

                    elif n_nod == 10:
                        row = fd.readline()
                        line.extend(int(row[i0:i1])
                                    for i0, i1 in idx[:2])
                        qtetras.append(line)

                    elif n_nod == 20:
                        row = fd.readline()
                        line.extend(int(row[i0:i1])
                                    for i0, i1 in idx[:12])
                        qhexas.append(line)

                    else:
                        raise ValueError('unsupported element type! (%d nodes)'
                                         % n_nod)

            elif kw == 'cmblock':
                if row[2].lower() != 'node': # Only node sets support.
                    continue

                n_nod = int(row[3].split('!')[0])
                fd.readline() # Format line not needed.

                nods = read_array(fd, n_nod, 1, nm.int32)
                nodal_bcs[row[1].strip()] = nods.ravel()

        fd.close()

        coors = nm.array(coors, dtype=nm.float64)

        tetras = nm.array(tetras, dtype=nm.int32)
        if len(tetras):
            mat_ids_tetras = tetras[:, 0]
            tetras = tetras[:, 1:]

        else:
            tetras.shape = (0, 4)
            mat_ids_tetras = nm.array([])

        hexas = nm.array(hexas, dtype=nm.int32)
        if len(hexas):
            mat_ids_hexas = hexas[:, 0]
            hexas = hexas[:, 1:]

        else:
            hexas.shape = (0, 8)
            mat_ids_hexas = nm.array([])

        if len(qtetras):
            qtetras = nm.array(qtetras, dtype=nm.int32)
            tetras.shape = (max(0, tetras.shape[0]), 4)
            tetras = nm.r_[tetras, qtetras[:, 1:5]]
            mat_ids_tetras = nm.r_[mat_ids_tetras, qtetras[:, 0]]

        if len(qhexas):
            qhexas = nm.array(qhexas, dtype=nm.int32)
            hexas.shape = (max(0, hexas.shape[0]), 8)
            hexas = nm.r_[hexas, qhexas[:, 1:9]]
            mat_ids_hexas = nm.r_[mat_ids_hexas, qhexas[:, 0]]

        if len(qtetras) or len(qhexas):
            ii = nm.union1d(tetras.ravel(), hexas.ravel())
            n_nod = len(ii)

            remap = nm.zeros((ii.max()+1,), dtype=nm.int32)
            remap[ii] = nm.arange(n_nod, dtype=nm.int32)

            ic = nm.searchsorted(ids, ii)
            coors = coors[ic]

        else:
            n_nod = coors.shape[0]
            remap = nm.zeros((nm.array(ids).max() + 1,), dtype=nm.int32)
            remap[ids] = nm.arange(n_nod, dtype=nm.int32)

        # Convert tetras as degenerate hexas to true tetras.
        ii = nm.where((hexas[:, 2] == hexas[:, 3])
                      & (hexas[:, 4] == hexas[:, 5])
                      & (hexas[:, 4] == hexas[:, 6])
                      & (hexas[:, 4] == hexas[:, 7]))[0]

        if len(ii) == len(hexas):
            tetras = nm.r_[tetras, hexas[ii[:, None], [0, 1, 2, 4]]]
            mat_ids_tetras = nm.r_[mat_ids_tetras, mat_ids_hexas[ii]]

            hexas = nm.delete(hexas, ii, axis=0)
            mat_ids_hexas = nm.delete(mat_ids_hexas, ii)

        else:
            output('WARNING: mesh "%s" has both tetrahedra and hexahedra!'
                   % mesh.name)

        ngroups = nm.zeros(len(coors), dtype=nm.int32)

        mesh = mesh_from_groups(mesh, ids, coors, ngroups,
                                [], [], [], [],
                                tetras, mat_ids_tetras,
                                hexas, mat_ids_hexas, remap=remap)

        mesh.nodal_bcs = {}
        for key, nods in six.iteritems(nodal_bcs):
            nods = nods[nods < len(remap)]
            mesh.nodal_bcs[key] = remap[nods]

        return mesh

class XYZMeshIO(MeshIO):
    """
    Trivial XYZ format working only with coordinates (in a .XYZ file) and the
    connectivity stored in another file with the same base name and .IEN
    suffix.
    """
    format = 'xyz'

    def _read_coors(self):
        coors = nm.loadtxt(self.filename)
        if (coors[:, -1] == 0).all():
            coors = coors[:, :-1].copy()

        return coors

    def read_dimension(self, ret_fd=False):
        coors = self._read_coors()
        dim = coors.shape[1]

        if ret_fd:
            fd = open(self.filename, 'r')
            return dim, fd

        else:
            return dim

    def read_bounding_box(self, ret_fd=False, ret_dim=False):
        coors = self._read_coors()
        bbox = nm.vstack((nm.amin(coors, 0), nm.amax(coors, 0)))

        if ret_fd: fd = open(self.filename, 'r')
        if ret_dim:
            dim = coors.shape[1]
            if ret_fd:
                return bbox, dim, fd
            else:
                return bbox, dim
        else:
            if ret_fd:
                return bbox, fd
            else:
                return bbox

    def read(self, mesh, omit_facets=False, **kwargs):
        coors = self._read_coors()
        n_nod, dim = coors.shape

        conn_ext = '.IEN' if op.splitext(self.filename)[1].isupper() else '.ien'
        conn = nm.loadtxt(edit_filename(self.filename, new_ext=conn_ext),
                          dtype=nm.int32) - 1
        desc = '%d_%d' % (dim, conn.shape[1])

        mesh._set_io_data(coors, nm.zeros(n_nod, dtype=nm.int32),
                          [conn], [nm.zeros(conn.shape[0], dtype=nm.int32)],
                          [desc])
        return mesh

    def write(self, filename, mesh, out=None, **kwargs):
        coors, ngroups, conns, mat_ids, desc = mesh._get_io_data()
        n_nod, dim = coors.shape

        zz = nm.zeros((n_nod, 3-dim))
        nm.savetxt(filename, nm.c_[coors, zz])

        conn_ext = '.IEN' if op.splitext(filename)[1].isupper() else '.ien'
        nm.savetxt(edit_filename(self.filename, new_ext=conn_ext),
                   conns[0] + 1)

        if out is not None:
            raise NotImplementedError


var_dict = list(vars().items())
io_table = {}

for key, var in var_dict:
    try:
        if is_derived_class(var, MeshIO):
            io_table[var.format] = var
    except TypeError:
        pass
del var_dict

def any_from_filename(filename, prefix_dir=None, file_format=None, mode='r'):
    """
    Create a MeshIO instance according to the kind of `filename`.

    Parameters
    ----------
    filename : str, function or MeshIO subclass instance
        The name of the mesh file. It can be also a user-supplied function
        accepting two arguments: `mesh`, `mode`, where `mesh` is a Mesh
        instance and `mode` is one of 'read','write', or a MeshIO subclass
        instance.
    prefix_dir : str
        The directory name to prepend to `filename`.

    Returns
    -------
    io : MeshIO subclass instance
        The MeshIO subclass instance corresponding to the kind of `filename`.
    """
    if not isinstance(filename, basestr):
        if isinstance(filename, MeshIO):
            return filename

        else:
            return UserMeshIO(filename)

    if prefix_dir is not None:
        filename = op.normpath(op.join(prefix_dir, filename))

    kwargs = {}
    if file_format is not None:
        if file_format in supported_formats:
            io_class = supported_formats[file_format][0]
            kwargs['file_format'] = file_format
        else:
            raise ValueError('unknown mesh format! (%s)' % file_format)
    else:
        ext2io = {e: (v[0], k) for k, v in supported_formats.items()
            for e in v[1] if '*' not in v[2]}
        ext = op.splitext(filename)[1].lower()
        if ext in ext2io:
            io_class = ext2io[ext][0]
            file_format = ext2io[ext][1]
            kwargs['file_format'] = file_format
        else:
            raise ValueError('unknown mesh format! (%s)' % ext)

    if mode == 'w' and 'w' not in supported_formats[file_format][2]:
        output('writable mesh formats:')
        output_mesh_formats('w')
        msg = 'write support not implemented for output mesh format "%s",' \
              ' see above!' % file_format
        raise ValueError(msg)

    return io_table[io_class](filename, **kwargs)


insert_static_method(MeshIO, any_from_filename)
del any_from_filename

# for meshio >= 4.0.3 # for 4.0.3 > meshio > 4.0.0 # format name: IO class, suffix, modes[, variants] # modes: r = read, w = write, c = test cell groups, v = test vertex groups # label format variants Split connectivities (columns except the last ones in `conns_in`) from cell     groups (the last columns of `conns_in`). Convert complex values in the output dictionary `out_in` to pairs of     real and imaginary parts. The abstract class for importing and exporting meshes.      Read the docstring of the Mesh() class. Basically all you need to do is to     implement the read() method::          def read(self, mesh, **kwargs):             nodes = ...             ngroups = ...             conns = ...             mat_ids = ...             descs = ...             mesh._set_io_data(nodes, ngroups, conns, mat_ids, descs)             return mesh      See the Mesh class' docstring how the nodes, ngroups, conns, mat_ids and     descs should look like. You just need to read them from your specific     format from disk.      To write a mesh to disk, just implement the write() method and use the     information from the mesh instance (e.g. nodes, conns, mat_ids and descs)     to construct your specific format.      Optionally, subclasses can implement read_data() to read also computation     results. This concerns mainly the subclasses with implemented write()     supporting the 'out' kwarg.      The default implementation od read_last_step() just returns 0. It should be     reimplemented in subclasses capable of storing several steps. The default implementation: just return 0 as the last step. Read true time step data from individual time steps.          Returns         -------         steps : array             The time steps.         times : array             The times of the time steps.         nts : array             The normalized times of the time steps, in [0, 1].          Notes         -----         The default implementation returns empty arrays. Special MeshIO subclass that enables reading and writing a mesh using a     user-supplied function. # meshio.__version__ > 3.3.2 # meshio.__version__ > 3.3.2 # Rearrange element node order to match SfePy. # Rearrange element node order to match SfePy. # Skip parameters. # Skip up/down pairs. #%d\n\n" % ig) # type name\n\n\n" % type_name) # number of nodes per element\n" % npe) # number of elements\n" % conn.shape[0]) # Zero based # number of parameter values per element\n" # Top level always 0? # number of parameters\n") # number of domains\n" # Domains in comsol have to be > 0 # number of up/down pairs\n") # Header # number of tags\n") # number of types\n") # Record # version unused serializable # class\n") # version\n") # sdim\n" % dim) # number of mesh points\n" % n_nod) # lowest mesh point index\n\n") # Always zero in SfePy # number of element types\n\n\n" % len(conns)) # TODO: Verify number of parameters for tri element # TODO: Verify number of parameters for tet element Read the mesh from a HDF5 file.          filename: str or tables.File             The HDF5 file to read the mesh from.         group: tables.group.Group or str, optional             The HDF5 file group to read the mesh from.             If None, the root group is used.         mesh: sfepy.dicrete.fem.Mesh or None             If None, the new mesh is created and returned, otherwise             content of this argument is replaced by the read mesh.          Returns         -------         sfepy.dicrete.fem.Mesh             readed mesh Write mesh to a hdf5 file.          filename: str or tables.File             The HDF5 file to write the mesh to.         group: tables.group.Group or None or str             The HDF5 file group to write the mesh to.             If None, the root group is used.             The group can be given as a path from root, e.g. /path/to/mesh         mesh: sfepy.dicrete.fem.Mesh             The mesh to write. # A new file. # Existing file. Read true time step data from individual time steps.          Returns         -------         steps : array             The time steps.         times : array             The times of the time steps.         nts : array             The normalized times of the time steps, in [0, 1]. # cell data. # read the whole file: # substract 1 from all elements, because we count from 0: Reads one non empty line (if it's a comment, it skips it). # comment or an empty line Reads one section from the mesh3d file.          integer ... if True, all numbers are passed to int(), otherwise to             float(), before returning          Some examples how a section can look like:          2         1 2 5 4 7 8 11 10         2 3 6 5 8 9 12 11          or          5         1 2 3 4     1         1 2 6 5     1         2 3 7 6     1         3 4 8 7     1         4 1 5 8     1          or          0 # Remove empty groups. # Solid keyword -> 3, otherwise 1 is the starting coors index. # Material id. # Number of nodes in line. # Element number. # Only node sets support. # Format line not needed. # Convert tetras as degenerate hexas to true tetras. Trivial XYZ format working only with coordinates (in a .XYZ file) and the     connectivity stored in another file with the same base name and .IEN     suffix. Create a MeshIO instance according to the kind of `filename`.      Parameters     ----------     filename : str, function or MeshIO subclass instance         The name of the mesh file. It can be also a user-supplied function         accepting two arguments: `mesh`, `mode`, where `mesh` is a Mesh         instance and `mode` is one of 'read','write', or a MeshIO subclass         instance.     prefix_dir : str         The directory name to prepend to `filename`.      Returns     -------     io : MeshIO subclass instance         The MeshIO subclass instance corresponding to the kind of `filename`.