Source code for ogs5py.fileclasses.msh.core

# -*- coding: utf-8 -*-
"""Core module for the ogs5py mesh file."""
from copy import deepcopy as dcp

import numpy as np

from ogs5py.fileclasses.base import File
from ogs5py.fileclasses.msh import generator as gen
from ogs5py.fileclasses.msh.checker import check_mesh_dict, check_mesh_list
from ogs5py.fileclasses.msh.msh_io import (
    export_mesh,
    import_mesh,
    load_ogs5msh,
    remove_dim,
    save_ogs5msh,
)
from ogs5py.fileclasses.msh.tools import (
    combine,
    gen_std_elem_id,
    gen_std_mat_id,
    get_centroids,
    get_mesh_center,
    get_node_centroids,
    get_volumes,
    no_of_elements,
    rotate_mesh,
    shift_mesh,
    transform_mesh,
)
from ogs5py.tools.types import ELEM_NAMES, EMPTY_MSH


class MSHsgl(File):
    """
    Class for a single mesh file.

    Parameters
    ----------
    mesh_dict : dict or None, optional
        Dictionary contains one '#FEM_MSH' block of the mesh file with
        the following information (sorted by keys):
            - mesh_data : dictionary
                contains information about (sorted by keys):
                    - AXISYMMETRY: bool (just true, otherwise not present)
                    - CROSS_SECTION: bool (just true, otherwise not present)
                    - PCS_TYPE: string
                    - GEO_TYPE: str
                    - GEO_NAME: string
                    - LAYER: int
            - nodes : ndarray
                Array with all node postions
            - elements : dictionary
                contains array of nodelists for elements sorted by element type
            - material_id : dictionary
                contains material ids for each element sorted by element type
    task_root : str, optional
        Path to the destiny model folder.
        Default: cwd+"ogs5model"
    task_id : str, optional
        Name for the ogs task.
        Default: "model"
    """

    def __init__(self, mesh_dict=None, **OGS_Config):
        super().__init__(**OGS_Config)
        self.file_ext = ".msh"
        self.force_writing = True

        if mesh_dict is None:
            self.__dict = EMPTY_MSH
        else:
            if isinstance(mesh_dict, list):
                raise ValueError(
                    "The given mesh may contains a multi-layer "
                    + "mesh. Try loading within "
                    + "the multimesh class."
                )
            else:
                if check_mesh_dict(mesh_dict):
                    self.__dict = mesh_dict
                else:
                    print("given mesh is not valid")
        self._block = 0

    # Pretend that there is a main keyword in the standard BASE-FORMAT
    @property
    def is_empty(self):
        """State if file is empty."""
        return not bool(self._meshlist[0]["elements"])

    #######################
    ### meshlist
    #######################
    # this is a workaround to make multi-layer and single-layer meshes usable
    @property
    def _meshlist(self):
        """list: mesh as list of dicts."""
        return [self._dict]

    @_meshlist.setter
    def _meshlist(self, value):
        self._dict = value[0]

    @_meshlist.deleter
    def _meshlist(self):
        self._dict = EMPTY_MSH

    @property
    def _dict(self):
        return self.__dict

    @_dict.setter
    def _dict(self, value):
        self.__dict = value

    #######################
    ### AXISYMMETRY
    #######################
    @property
    def AXISYMMETRY(self):
        """bool: AXISYMMETRY attribute."""
        if "AXISYMMETRY" in self._dict["mesh_data"]:
            return self._dict["mesh_data"]["AXISYMMETRY"]
        return False

    @AXISYMMETRY.setter
    def AXISYMMETRY(self, value):
        if value:
            self._dict["mesh_data"]["AXISYMMETRY"] = True
        else:
            del self.AXISYMMETRY

    @AXISYMMETRY.deleter
    def AXISYMMETRY(self):
        if "AXISYMMETRY" in self._dict["mesh_data"]:
            del self._dict["mesh_data"]["AXISYMMETRY"]

    #######################
    ### CROSS_SECTION
    #######################
    @property
    def CROSS_SECTION(self):
        """bool: CROSS_SECTION attribute."""
        if "CROSS_SECTION" in self._dict["mesh_data"]:
            return self._dict["mesh_data"]["CROSS_SECTION"]
        return False

    @CROSS_SECTION.setter
    def CROSS_SECTION(self, value):
        if value:
            self._dict["mesh_data"]["CROSS_SECTION"] = True
        else:
            del self.CROSS_SECTION

    @CROSS_SECTION.deleter
    def CROSS_SECTION(self):
        if "CROSS_SECTION" in self._dict["mesh_data"]:
            del self._dict["mesh_data"]["CROSS_SECTION"]

    #######################
    ### PCS_TYPE
    #######################
    @property
    def PCS_TYPE(self):
        """str: PCS_TYPE."""
        if "PCS_TYPE" in self._dict["mesh_data"]:
            return self._dict["mesh_data"]["PCS_TYPE"]
        return None

    @PCS_TYPE.setter
    def PCS_TYPE(self, value):
        if value is not None:
            value = str(value)
        if value:
            self._dict["mesh_data"]["PCS_TYPE"] = value
        else:
            del self.PCS_TYPE

    @PCS_TYPE.deleter
    def PCS_TYPE(self):
        if "PCS_TYPE" in self._dict["mesh_data"]:
            del self._dict["mesh_data"]["PCS_TYPE"]

    #######################
    ### GEO_NAME
    #######################
    @property
    def GEO_NAME(self):
        """str: GEO_NAME."""
        if "GEO_NAME" in self._dict["mesh_data"]:
            return self._dict["mesh_data"]["GEO_NAME"]
        return None

    @GEO_NAME.setter
    def GEO_NAME(self, value):
        if value is not None:
            value = str(value)
        if value:
            self._dict["mesh_data"]["GEO_NAME"] = value
        else:
            del self.GEO_NAME

    @GEO_NAME.deleter
    def GEO_NAME(self):
        if "GEO_NAME" in self._dict["mesh_data"]:
            del self._dict["mesh_data"]["GEO_NAME"]
            del self.GEO_TYPE

    #######################
    ### GEO_TYPE
    #######################
    @property
    def GEO_TYPE(self):
        """str: GEO_TYPE."""
        if "GEO_TYPE" in self._dict["mesh_data"]:
            return self._dict["mesh_data"]["GEO_TYPE"]
        return None

    @GEO_TYPE.setter
    def GEO_TYPE(self, value):
        if value is not None:
            value = str(value)
        if value:
            self._dict["mesh_data"]["GEO_TYPE"] = value
        else:
            del self.GEO_TYPE

    @GEO_TYPE.deleter
    def GEO_TYPE(self):
        if "GEO_TYPE" in self._dict["mesh_data"]:
            del self._dict["mesh_data"]["GEO_TYPE"]

    #######################
    ### LAYER
    #######################
    @property
    def LAYER(self):
        """int: LAYER."""
        if "LAYER" in self._dict["mesh_data"]:
            return self._dict["mesh_data"]["LAYER"]
        return None

    @LAYER.setter
    def LAYER(self, value):
        if value is not None:
            self._dict["mesh_data"]["LAYER"] = int(value)
        else:
            del self.LAYER

    @LAYER.deleter
    def LAYER(self):
        if "LAYER" in self._dict["mesh_data"]:
            del self._dict["mesh_data"]["LAYER"]

    #######################
    ### NODES
    #######################
    @property
    def NODES(self):
        """ndarray: (n,3) NODES of the mesh by its xyz-coordinates."""
        return self._dict["nodes"]

    @NODES.setter
    def NODES(self, value):
        self._dict["nodes"] = np.array(value)

    @NODES.deleter
    def NODES(self):
        self._dict["nodes"] = None
        del self.ELEMENTS
        del self.MATERIAL_ID
        del self.ELEMENT_ID

    #######################
    ### NODE_NO
    #######################
    @property
    def NODE_NO(self):
        """int: number of NODES."""
        if self.NODES is None:
            return 0
        return self.NODES.shape[0]

    #######################
    ### ELEMENTS
    #######################
    @property
    def ELEMENTS(self):
        """
        Get and set the ELEMENTS of the mesh.

        Notes
        -----
        Type : dict of ndarrays
            The elements are a dictionary sorted by their element-type

                "line" : ndarray of shape (n_line,2)
                    1D element with 2 nodes
                "tri" : ndarray of shape (n_tri,3)
                    2D element with 3 nodes
                "quad" : ndarray of shape (n_quad,4)
                    2D element with 4 nodes
                "tet" : ndarray of shape (n_tet,4)
                    3D element with 4 nodes
                "pyra" : ndarray of shape (n_pyra,5)
                    3D element with 5 nodes
                "pris" : ndarray of shape (n_pris,6)
                    3D element with 6 nodes
                "hex" : ndarray of shape (n_hex,8)
                    3D element with 8 nodes
        """
        return self._dict["elements"]

    @ELEMENTS.setter
    def ELEMENTS(self, value):
        self._dict["elements"] = value

    @ELEMENTS.deleter
    def ELEMENTS(self):
        self._dict["elements"] = {}
        del self.MATERIAL_ID
        del self.ELEMENT_ID

    #######################
    ### ELEMENT_TYPES
    #######################
    @property
    def ELEMENT_TYPES(self):
        """:class:`set`: ELEMENT types of the mesh."""
        out = []
        for elem in ELEM_NAMES:
            if elem not in self.ELEMENTS:
                continue
            out.append(elem)
        # return set(self.ELEMENTS)
        return out

    #######################
    ### ELEMENT_ID
    #######################
    @property
    def ELEMENT_ID(self):
        """
        Get and set the ELEMENT_IDs of the mesh.

        Standard element id order is given by:

            "line" "tri" "quad" "tet" "pyra" "pris" "hex"

        Notes
        -----
        Type : dict of ndarrays
            The element IDs are a dictionary containing ints
            sorted by their element-type

                "line" : ndarray of shape (n_line,)
                    1D element with 2 nodes
                "tri" : ndarray of shape (n_tri,)
                    2D element with 3 nodes
                "quad" : ndarray of shape (n_quad,)
                    2D element with 4 nodes
                "tet" : ndarray of shape (n_tet,)
                    3D element with 4 nodes
                "pyra" : ndarray of shape (n_pyra,)
                    3D element with 5 nodes
                "pris" : ndarray of shape (n_pris,)
                    3D element with 6 nodes
                "hex" : ndarray of shape (n_hex,)
                    3D element with 8 nodes
        """
        return self._dict["element_id"]

    @ELEMENT_ID.setter
    def ELEMENT_ID(self, value):
        self._dict["element_id"] = value

    @ELEMENT_ID.deleter
    def ELEMENT_ID(self):
        self._dict["element_id"] = gen_std_elem_id(self.ELEMENTS)

    #######################
    ### ELEMENT_NO
    #######################
    @property
    def ELEMENT_NO(self):
        """int: number of ELEMENTS."""
        return no_of_elements(self._dict)

    #######################
    ### MATERIAL_ID
    #######################
    @property
    def MATERIAL_ID(self):
        """
        Get and set the MATERIAL_IDs of the mesh.

        Notes
        -----
        Type : dict of ndarrays
            The material IDs are a dictionary containing ints
            sorted by their element-type

                "line" : ndarray of shape (n_line,)
                    1D element with 2 nodes
                "tri" : ndarray of shape (n_tri,)
                    2D element with 3 nodes
                "quad" : ndarray of shape (n_quad,)
                    2D element with 4 nodes
                "tet" : ndarray of shape (n_tet,)
                    3D element with 4 nodes
                "pyra" : ndarray of shape (n_pyra,)
                    3D element with 5 nodes
                "pris" : ndarray of shape (n_pris,)
                    3D element with 6 nodes
                "hex" : ndarray of shape (n_hex,)
                    3D element with 8 nodes
        """
        return self._dict["material_id"]

    @MATERIAL_ID.setter
    def MATERIAL_ID(self, value):
        if isinstance(value, int):
            self._dict["material_id"] = gen_std_mat_id(self.ELEMENTS, value)
        else:
            self._dict["material_id"] = value

    @MATERIAL_ID.deleter
    def MATERIAL_ID(self):
        self._dict["material_id"] = gen_std_mat_id(self.ELEMENTS)

    @property
    def MATERIAL_ID_flat(self):
        """
        Get flat version of the MATERIAL_IDs of the mesh.

        See "mesh.MATERIAL_ID"
        This flattend MATERIAL_IDs are a stacked version of MATERIAL_ID, to get
        one continous array. They are stacked in order of the ELEMENT_IDs.
        Standard stack order is given by:

            "line" "tri" "quad" "tet" "pyra" "pris" "hex"

        Notes
        -----
        Type : ndarray
            The centroids are a list containing xyz-coordiantes
        """
        out = np.empty(self.ELEMENT_NO, dtype=int)
        for elem in ELEM_NAMES:
            if elem not in self.ELEMENTS:
                continue
            out[self.ELEMENT_ID[elem]] = self.MATERIAL_ID[elem]
        return out

    #######################
    ### centroids
    #######################
    @property
    def centroids(self):
        """
        Get the centroids of the mesh.

        Notes
        -----
        Type : dict of ndarrays
            The centroids are a dictionary containing xyz-coordiantes
            sorted by their element-type

                "line" : ndarray of shape (n_line,3)
                    1D element with 2 nodes
                "tri" : ndarray of shape (n_tri,3)
                    2D element with 3 nodes
                "quad" : ndarray of shape (n_quad,3)
                    2D element with 4 nodes
                "tet" : ndarray of shape (n_tet,3)
                    3D element with 4 nodes
                "pyra" : ndarray of shape (n_pyra,3)
                    3D element with 5 nodes
                "pris" : ndarray of shape (n_pris,3)
                    3D element with 6 nodes
                "hex" : ndarray of shape (n_hex,3)
                    3D element with 8 nodes
        """
        return get_centroids(self._dict)

    @property
    def centroids_flat(self):
        """
        Get flat version of the centroids of the mesh.

        See the "mesh.get_centroids" method.
        This flattend centroids are a stacked version of centroids, to get
        one continous array. They are stacked in order of the element ids.
        Standard stack order is given by:

            "line" "tri" "quad" "tet" "pyra" "pris" "hex"

        Notes
        -----
        Type : ndarray
            The centroids are a list containing xyz-coordiantes
        """
        # just call centroids once
        tmp = dcp(self.centroids)
        out = np.empty((self.ELEMENT_NO, 3), dtype=float)
        for elem in ELEM_NAMES:
            if elem not in self.ELEMENTS:
                continue
            out[self.ELEMENT_ID[elem]] = tmp[elem]
        return out

    @property
    def node_centroids(self):
        """
        Get the node centroids of the mesh.

        Notes
        -----
        Type : dict of ndarrays
            The centroids are a dictionary containing xyz-coordiantes
            sorted by their element-type

                "line" : ndarray of shape (n_line,3)
                    1D element with 2 nodes
                "tri" : ndarray of shape (n_tri,3)
                    2D element with 3 nodes
                "quad" : ndarray of shape (n_quad,3)
                    2D element with 4 nodes
                "tet" : ndarray of shape (n_tet,3)
                    3D element with 4 nodes
                "pyra" : ndarray of shape (n_pyra,3)
                    3D element with 5 nodes
                "pris" : ndarray of shape (n_pris,3)
                    3D element with 6 nodes
                "hex" : ndarray of shape (n_hex,3)
                    3D element with 8 nodes
        """
        return get_node_centroids(self._dict)

    @property
    def node_centroids_flat(self):
        """
        Get flat version of the node centroids of the mesh.

        See the "mesh.get_centroids" method.
        This flattend centroids are a stacked version of centroids, to get
        one continous array. They are stacked in order of the element ids.
        Standard stack order is given by:

            "line" "tri" "quad" "tet" "pyra" "pris" "hex"

        Notes
        -----
        Type : ndarray
            The centroids are a list containing xyz-coordiantes
        """
        # just call centroids once
        tmp = dcp(self.node_centroids)
        out = np.empty((self.ELEMENT_NO, 3), dtype=float)
        for elem in ELEM_NAMES:
            if elem not in self.ELEMENTS:
                continue
            out[self.ELEMENT_ID[elem]] = tmp[elem]
        return out

    #######################
    ### volumes
    #######################
    @property
    def volumes(self):
        """
        Get the volumes of the mesh-elements.

        Notes
        -----
        Type : dict of ndarrays
            The volumes are a dictionary containing the n-dimension volumes
            sorted by their element-type

                "line" : ndarray of shape (n_line,3)
                    1D element with 2 nodes
                "tri" : ndarray of shape (n_tri,3)
                    2D element with 3 nodes
                "quad" : ndarray of shape (n_quad,3)
                    2D element with 4 nodes
                "tet" : ndarray of shape (n_tet,3)
                    3D element with 4 nodes
                "pyra" : ndarray of shape (n_pyra,3)
                    3D element with 5 nodes
                "pris" : ndarray of shape (n_pris,3)
                    3D element with 6 nodes
                "hex" : ndarray of shape (n_hex,3)
                    3D element with 8 nodes
        """
        return get_volumes(self._dict)

    @property
    def volumes_flat(self):
        """
        Get flat version of the volumes of the mesh-elements.

        This flattend volumes are a stacked version of centroids, to get
        one continous array. They are stacked in order of the element ids.
        Standard stack order is given by:

            "line" "tri" "quad" "tet" "pyra" "pris" "hex"

        Notes
        -----
        Type : ndarray
            The volumes are a list containing the n-dimensional element volume
        """
        # just call volumes once
        tmp = dcp(self.volumes)
        out = np.empty(self.ELEMENT_NO, dtype=float)
        for elem in ELEM_NAMES:
            if elem not in self.ELEMENTS:
                continue
            out[self.ELEMENT_ID[elem]] = tmp[elem]
        return out

    #######################
    ### center
    #######################
    @property
    def center(self):
        """Get the mesh center."""
        return get_mesh_center(self._dict)

    #######################
    ### Class methods
    #######################
    def reset(self):
        """Delete every content."""
        self._block = 0
        self._meshlist = [EMPTY_MSH]

    def load(self, filepath, **kwargs):
        r"""
        Load an OGS5 mesh from file.

        kwargs will be forwarded to "tools.load_ogs5msh"

        Parameters
        ----------
        filepath : string
            path to the '\*.msh' OGS5 mesh file to load
        verbose : bool, optional
            Print information of the reading process. Default: True
        ignore_unknown : bool, optional
            Unknown data in the file will be ignored. Default: False
        max_node_no : int, optional
            If you know the maximal node number per elements in the mesh file,
            you can optimise the reading a bit. By default the algorithm will
            assume hexahedrons as 'largest' elements in the mesh. Default: 8
        encoding : str or None, optional
            encoding of the given file. If ``None`` is given, the system
            standard is used. Default: ``None``

        Notes
        -----
        The $AREA keyword within the Nodes definition is NOT supported
        and will violate the read data if present.
        """
        #        kwargs["verbose"] = True
        tmp = load_ogs5msh(filepath, **kwargs)
        if isinstance(tmp, list) and not isinstance(self, MSH):
            print(
                filepath
                + " may contains a multi-layer mesh."
                + "Try loading within the multimesh class."
            )
        elif not isinstance(tmp, list):
            tmp = [tmp]
        if "verbose" in kwargs:
            verbose = kwargs["verbose"]
        else:
            verbose = False
        if check_mesh_list(tmp, verbose=verbose):
            self._block = 0
            self._meshlist = tmp
        else:
            raise ValueError("MSH: " + filepath + ": given mesh is not valid")

    def read_file(self, path, encoding=None, verbose=False):
        r"""
        Load an OGS5 mesh from file.

        Parameters
        ----------
        path : str
            path to the '\*.msh' OGS5 mesh file to load
        encoding : str or None, optional
            encoding of the given file. If ``None`` is given, the system
            standard is used. Default: ``None``
        verbose : bool, optional
            Print information of the reading process. Default: True
        """
        self.load(
            path, verbose=verbose, ignore_unknown=True, encoding=encoding
        )

    def set_dict(self, mesh_dict):
        """
        Set an mesh as returned by tools methods or generators.

        Mesh will be checked for validity.

        Parameters
        ----------
        mesh_dict : dict or None, optional
            Contains one '#FEM_MSH' block of an OGS5 mesh file
            with the following information (sorted by keys):

                mesh_data : dict
                    dictionary containing information about

                    - AXISYMMETRY (bool)
                    - CROSS_SECTION (bool)
                    - PCS_TYPE (str)
                    - GEO_TYPE (str)
                    - GEO_NAME (str)
                    - LAYER (int)

                nodes : ndarray
                    Array with all node postions
                elements : dict
                    contains nodelists for elements sorted by element types
                material_id : dict
                    contains material ids for each element sorted by types
                element_id : dict
                    contains element ids for each element sorted by types
        """
        if check_mesh_dict(mesh_dict):
            self._dict = mesh_dict
        else:
            print("given mesh is not valid")

    def save(self, path, **kwargs):
        r"""
        Save the mesh to an OGS5 mesh file.

        kwargs will be forwarded to "tools.save_ogs5msh"

        Parameters
        ----------
        path : string
            path to the '\*.msh' OGS5 mesh file to save
        verbose : bool, optional
            Print information of the writing process. Default: True
        """
        # no top-comment allowed in MSH file
        if "verbose" in kwargs:
            verbose = kwargs["verbose"]
        else:
            kwargs["verbose"] = verbose = False
        if self.check(verbose=verbose):
            save_ogs5msh(
                path,
                self._meshlist,
                top_com=None,
                bot_com=self.bot_com,
                **kwargs,
            )
        else:
            print("the mesh could not be saved since it is not valid")

    def import_mesh(self, filepath, **kwargs):
        """
        Import an external unstructured mesh from diffrent file-formats.

        kwargs will be forwarded to "tools.import_mesh"

        Parameters
        ----------
        filepath : string or meshio.Mesh instance
            mesh to import
        file_format : str, optional
            Here you can specify the fileformat. If 'None' it will be
            determined by file extension. Default: None
        ignore_unknown : bool, optional
            Unknown data in the file will be ignored. Default: False
        import_dim : iterable of int, optional
            State which elements should be imported by dimensionality.
            Can be used to sort out unneeded elements for example from gmsh.
            Default: (1, 2, 3)

        Notes
        -----
        This routine calls the 'read' function from the meshio package
        and converts the output (see here: https://github.com/nschloe/meshio)
        If there is any "vertex" (0D element) in the element data,
        it will be removed.
        """
        tmp = import_mesh(filepath, **kwargs)
        if check_mesh_dict(tmp):
            self._meshlist = [tmp]
        else:
            raise ValueError("MSH:" + filepath + ": given mesh is not valid")

    def export_mesh(self, filepath, verbose=False, **kwargs):
        """
        Export the mesh to an unstructured mesh in diffrent file-formats.

        kwargs will be forwarded to "tools.export_mesh"

        Parameters
        ----------
        filepath : string
            path to the file to export
        file_format : str, optional
            Here you can specify the fileformat. If 'None' it will be
            determined by file extension. Default: None
        export_material_id : bool, optional
            Here you can specify if the material_id should be exported.
            Default: True
        export_element_id : bool, optional
            Here you can specify if the element_id should be exported.
            Default: True
        cell_data_by_id : ndarray or dict, optional
            Here you can specify additional element data sorted by their IDs.
            It can be a dictionary with data-name as key and
            the ndarray as value.
            Default: None
        point_data : ndarray or dict, optional
            Here you can specify additional point data sorted by their IDs.
            It can be a dictionary with data-name as key and
            the ndarray as value.
            Default: None
        field_data : ndarray or dict, optional
            Here you can specify additional field data of the mesh.
            It can be a dictionary with data-name as key and
            the ndarray as value.
            Default: None

        Notes
        -----
        This routine calls the 'write' function from the meshio package
        and converts the input (see here: https://github.com/nschloe/meshio)
        """
        if self.check(verbose=verbose):
            export_mesh(filepath, self._dict, **kwargs)
        else:
            print("the mesh could not be exported since it is not valid")

    def combine_mesh(self, ext_mesh, **kwargs):
        """
        Combine this mesh with an external mesh.

        The node list will be
        updated to eliminate duplicates.
        Element intersections are not checked.
        kwargs will be forwarded to "tools.combine"

        Parameters
        ----------
        ext_mesh: mesh, dict or file
            This is the mesh that should be added to the existing one.
        decimals : int, optional
            Number of decimal places to round the nodes to (default: 3).
            This will not round the output, it is just for comparison of the
            node vectors.
        fast : bool, optional
            If fast is True, the vector comparison is executed by a
            decimal comparison. If fast is False, all pairwise distances
            are calculated. Default: False
        """
        if isinstance(ext_mesh, MSH):
            tmp_mesh = ext_mesh()
        elif isinstance(ext_mesh, dict):
            tmp_mesh = ext_mesh
        else:
            try:
                tmp_mesh = load_ogs5msh(ext_mesh)
            except Exception:
                try:
                    tmp_mesh = import_mesh(ext_mesh)
                except Exception:
                    print("Could not interpret the mesh that should be added")
                    return

        if check_mesh_dict(tmp_mesh, verbose=False):
            self._dict = combine(self._dict, tmp_mesh, **kwargs)
        else:
            print("given mesh to add is not valid")

    def check(self, verbose=True):
        """
        Check if the mesh is valid.

        Checked in the sence, that the contained data is consistent.
        Checks for correct element definitions or Node duplicates
        are not carried out.

        Parameters
        ----------
        verbose : bool, optional
            Print information for the executed checks. Default: True

        Returns
        -------
        result : bool
            Validity of the given mesh.
        """
        return check_mesh_list(self._meshlist, verbose=verbose)

    def swap_axis(self, axis1="y", axis2="z"):
        """
        Swap axis of the coordinate system.

        Parameters
        ----------
        axis1 : :class:`str` or :class:`int`, optional
            First selected Axis.
            Either in ["x", "y", "z"] or in [0, 1, 2]. Default: "y"
        axis2 : :class:`str` or :class:`int`, optional
            Second selected Axis.
            Either in ["x", "y", "z"] or in [0, 1, 2]. Default: "z"
        """
        axis = ["x", "y", "z"]
        if axis1 in range(3):
            axis1 = axis[axis1]
        if axis2 in range(3):
            axis2 = axis[axis2]
        if axis1 not in axis or axis2 not in axis:
            raise ValueError(
                "MSH.swap_axis: axis need to be 'x', 'y' or 'z': "
                + str(axis1)
                + ", "
                + str(axis2)
            )
        if axis1 == axis2:
            raise ValueError(
                "MSH.swap_axis: please select distict axis: "
                + str(axis1)
                + " = "
                + str(axis2)
            )
        ax1 = axis.index(axis1)
        ax2 = axis.index(axis2)
        self.NODES[:, [ax1, ax2]] = self.NODES[:, [ax2, ax1]]

    def rotate(
        self,
        angle,
        rotation_axis=(0.0, 0.0, 1.0),
        rotation_point=(0.0, 0.0, 0.0),
    ):
        """
        Rotate a given mesh around a given rotation axis with a given angle.

        Parameters
        ----------
        angle : float
            rotation angle given in radial length
        rotation_axis : array_like, optional
            Array containing the vector for ratation axis. Default: (0,0,1)
        rotation_point : array_like, optional
            Vector of the ratation base point. Default:(0,0,0)
        """
        rotate_mesh(self._dict, angle, rotation_axis, rotation_point)

    def shift(self, vector):
        """
        Shift a given mesh with a given vector.

        Parameters
        ----------
        vector : ndarray
            array containing the shifting vector
        """
        shift_mesh(self._dict, vector)

    def transform(self, xyz_func, **kwargs):
        """
        Transform a given mesh with a given function "xyz_func".

        kwargs will be forwarded to "xyz_func".

        Parameters
        ----------
        xyz_func : function
            the function transforming the points:
            ``x_new, y_new, z_new = f(x_old, y_old, z_old, **kwargs)``
        """
        transform_mesh(self._dict, xyz_func, **kwargs)

    def remove_dim(self, remove):
        """
        Remove elements by given dimensions from a mesh.

        Parameters
        ----------
        remove : iterable of int or single int
            State which elements should be removed by dimensionality (1, 2, 3).
        """
        remove_dim(self._dict, remove=remove)

    def generate(self, generator="rectangular", **kwargs):
        """
        Use a mesh-generator from the generator module.

        See: :any:`ogs5py.fileclasses.msh.generator`

        Parameters
        ----------
        generator : str
            set the generator from the generator module
        **kwargs
            kwargs will be forwarded to the generator in use

        Notes
        -----
        .. currentmodule:: ogs5py.fileclasses.msh.generator

        The following generators are available:

        .. autosummary::
           rectangular
           radial
           grid_adapter2D
           grid_adapter3D
           block_adapter3D
           gmsh
        """
        self._dict = getattr(gen, generator)(**kwargs)

    def show(
        self,
        show_cell_data=None,
        show_material_id=False,
        show_element_id=False,
        log_scale=False,
    ):
        """
        Display the mesh colored by its material ID.

        Parameters
        ----------
        show_cell_data : ndarray or dict, optional
            Here you can specify additional
            element/cell data sorted by their IDs.
            It can be a dictionary with data-name as key
            and the ndarray as value. Default: None
        show_material_id : bool, optional
            Here you can specify if the material_id should be shown.
            Default: False
        show_element_id : bool, optional
            Here you can specify if the element_id should be shown.
            Default: False
        log_scale : bool, optional
            State if the cell_data should be shown in log scale.
            Default: False

        Notes
        -----
        This routine needs "mayavi" to display the mesh.
        (see here: https://github.com/enthought/mayavi)
        """
        from ogs5py.fileclasses.msh.viewer import show_mesh

        show_mesh(
            self._dict,
            show_cell_data,
            show_material_id,
            show_element_id,
            log_scale,
        )

    def set_material_id(
        self, material_id=0, element_id=None, element_mask=None
    ):
        """
        Set material IDs by the corresponding element IDs.

        Parameters
        ----------
        material_id : :class:`int` or ndarray, optional
            The new material IDs. Either one value or an array.
            Default: 0
        element_id : ndarray or None, optional
            The corresponding element IDs, where to set the material IDs.
            If None, all elements are assumed and the material IDs are added
            by their index.
            Default: None
        element_mask : ndarray or None, optional
            Instead of the element IDs, one can specify a mask to select the
            element IDs.
            Default: None
        """
        if np.ndim(np.squeeze(material_id)) == 0:
            material_id = int(np.squeeze(material_id))
        if (
            isinstance(material_id, int)
            and element_id is None
            and element_mask is None
        ):
            self.MATERIAL_ID = material_id
        else:
            if element_id is None and element_mask is None:
                element_id = np.arange(self.ELEMENT_NO, dtype=int)
            else:
                if element_mask is None:
                    element_id = np.array(element_id, ndmin=1, dtype=int)
                else:
                    element_mask = np.array(element_mask, ndmin=1, dtype=bool)
                    assert len(element_mask) == self.ELEMENT_NO, "Wrong mask."
                    element_id = np.nonzero(element_mask)[0]
                assert (
                    np.max(element_id) < self.ELEMENT_NO
                ), "Element ID out of range."
                assert np.max(element_id) >= 0, "Element ID out of range."
                assert len(element_id) == len(
                    np.unique(element_id)
                ), "Element IDs not unique."
            if not isinstance(material_id, int):
                material_id = np.array(material_id, ndmin=1, dtype=int)
                assert len(element_id) == len(
                    material_id
                ), "Different number of material and element IDs."
            mat_id = self.MATERIAL_ID_flat
            mat_id[element_id] = material_id
            for elem in self.ELEMENT_TYPES:
                self._dict["material_id"][elem] = mat_id[self.ELEMENT_ID[elem]]

    #######################
    ### Special methods
    #######################
    def __call__(self):
        """Get a copy of the underlying dictionary."""
        return dcp(self._dict)

    def __repr__(self):
        """Representation."""
        out = "#FEM_MSH\n"
        if self.AXISYMMETRY:
            out += " $AXISYMMETRY\n"
        if self.CROSS_SECTION:
            out += " $CROSS_SECTION\n"
        if self.PCS_TYPE is not None:
            out += " $PCS_TYPE\n"
            out += "  " + self.PCS_TYPE + "\n"
        if self.GEO_NAME is not None:
            out += " $GEO_NAME\n"
            out += "  " + self.GEO_NAME + "\n"
        if self.GEO_TYPE is not None:
            out += " $GEO_TYPE\n"
            out += "  " + self.GEO_TYPE + " " + self.GEO_NAME + "\n"
        if self.LAYER is not None:
            out += " $LAYER\n"
            out += "  " + str(self.LAYER) + "\n"
        out += " $NODES\n"
        if self.NODES is None:
            out += "  None\n"
        else:
            out += "  " + str(self.NODES.shape[0]) + "\n"
            out += "   ...\n"
        out += " $ELEMENTS\n"
        if self.ELEMENTS is None:
            out += "  None\n"
        else:
            elem_no = 0
            for elem in self.ELEMENTS:
                elem_no += self.ELEMENTS[elem].shape[0]
            out += "  " + str(elem_no) + "\n"
            out += "   ...\n"
        out += "#STOP\n"
        return out


[docs]class MSH(MSHsgl): """ Class for a multi layer mesh file that contains multiple '#FEM_MSH' Blocks. Parameters ---------- mesh_list : list of dict or None, optional each dictionary contains one '#FEM_MSH' block of the mesh file with with the following information (sorted by keys): mesh_data : dict dictionary containing information about - AXISYMMETRY (bool) - CROSS_SECTION (bool) - PCS_TYPE (str) - GEO_TYPE (str) - GEO_NAME (str) - LAYER (int) nodes : ndarray Array with all node postions elements : dict contains nodelists for elements sorted by element types material_id : dict contains material ids for each element sorted by element types element_id : dict contains element ids for each element sorted by element types task_root : str, optional Path to the destiny model folder. Default: cwd+"ogs5model" task_id : str, optional Name for the ogs task. Default: "model" """ def __init__(self, mesh_list=None, **OGS_Config): super().__init__(None, **OGS_Config) if mesh_list is None: self.__meshlist = [EMPTY_MSH] else: self.__meshlist = mesh_list # meshlist (override the property of MSHsgl) # @property @MSHsgl._meshlist.getter def _meshlist(self): return self.__meshlist @_meshlist.setter def _meshlist(self, value): self.__meshlist = value if len(self.__meshlist) >= self._block: self._block = 0 @_meshlist.deleter def _meshlist(self): self.__meshlist = [EMPTY_MSH] self._block = 0 # override the _dict property @MSHsgl._dict.getter def _dict(self): return self.__meshlist[self.block] @_dict.setter def _dict(self, value): self.__meshlist[self.block] = value @property def block_no(self): """:class:`int`: The number of blocks in the file.""" return len(self.__meshlist) # select the block to be edited @property def block(self): """:class:`int`: The actual block to access in the file.""" return self._block @block.setter def block(self, value): value = int(value) if 0 <= value < len(self.__meshlist): self._block = value if -len(self.__meshlist) <= value < 0: self._block = len(self.__meshlist) - value if value == len(self.__meshlist): self._block = value self.__meshlist.append(EMPTY_MSH) @block.deleter def block(self): self._block = 0 def __repr__(self): """Representation.""" out = "" old_block = self.block for i in range(len(self._meshlist)): self.block = i out += "#FEM_MSH\n" if self.AXISYMMETRY: out += " $AXISYMMETRY\n" if self.CROSS_SECTION: out += " $CROSS_SECTION\n" if self.PCS_TYPE is not None: out += " $PCS_TYPE\n" out += " " + self.PCS_TYPE + "\n" if self.GEO_NAME is not None: out += " $GEO_NAME\n" out += " " + self.GEO_NAME + "\n" if self.GEO_TYPE is not None: out += " $GEO_TYPE\n" out += " " + self.GEO_TYPE + " " + self.GEO_NAME + "\n" if self.LAYER is not None: out += " $LAYER\n" out += " " + str(self.LAYER) + "\n" out += " $NODES\n" if self.NODES is None: out += " None\n" else: out += " " + str(self.NODES.shape[0]) + "\n" out += " ...\n" out += " $ELEMENTS\n" if self.ELEMENTS is None: out += " None\n" else: out += " " + str(self.ELEMENT_NO) + "\n" out += " ...\n" if self._meshlist: out += "#STOP\n" self.block = old_block return out