mesh.py

from __future__ import annotations
import warnings
from abc import ABC, abstractmethod
from collections.abc import Iterable, Sequence, Mapping
from functools import wraps
from math import pi, sqrt, atan2
from numbers import Integral, Real
from pathlib import Path
from typing import Protocol

import h5py
import lxml.etree as ET
import numpy as np
from pathlib import Path

import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.utility_funcs import change_directory
from .bounding_box import BoundingBox
from ._xml import get_elem_list, get_text
from .mixin import IDManagerMixin
from .surface import _BOUNDARY_TYPES
from .utility_funcs import input_path


class MeshMaterialVolumes(Mapping):
    """Results from a material volume in mesh calculation.

    This class provides multiple ways of accessing information about material
    volumes in individual mesh elements. First, the class behaves like a
    dictionary that maps material IDs to an array of volumes equal in size to
    the number of mesh elements. Second, the class provides a :meth:`by_element`
    method that gives all the material volumes for a specific mesh element.

    .. versionadded:: 0.15.1

    Parameters
    ----------
    materials : numpy.ndarray
        Array of shape (elements, max_materials) storing material IDs
    volumes : numpy.ndarray
        Array of shape (elements, max_materials) storing material volumes
    bboxes : numpy.ndarray, optional
        Array of shape (elements, max_materials, 6) storing axis-aligned
        bounding boxes for each (element, material) combination with ordering
        (xmin, ymin, zmin, xmax, ymax, zmax). Bounding boxes enclose the
        ray-estimator prisms used to compute volumes.

    See Also
    --------
    openmc.MeshBase.material_volumes

    Examples
    --------
    If you want to get the volume of a specific material in every mesh element,
    index the object with the material ID:

    >>> volumes = mesh.material_volumes(...)
    >>> volumes
    {1: <32121 nonzero volumes>
     2: <338186 nonzero volumes>
     3: <49120 nonzero volumes>}

    If you want the volume of all materials in a specific mesh element, use the
    :meth:`by_element` method:

    >>> volumes = mesh.material_volumes(...)
    >>> volumes.by_element(42)
    [(2, 31.87963824195591), (1, 6.129949130817542)]

    """
    def __init__(
        self,
        materials: np.ndarray,
        volumes: np.ndarray,
        bboxes: np.ndarray | None = None
    ):
        self._materials = materials
        self._volumes = volumes
        self._bboxes = bboxes

        if self._bboxes is not None:
            if self._bboxes.shape[:2] != self._materials.shape:
                raise ValueError(
                    'bboxes must have shape (elements, max_materials, 6) '
                    'matching materials/volumes.'
                )
            if self._bboxes.shape[2] != 6:
                raise ValueError(
                    'bboxes must have shape (elements, max_materials, 6).'
                )

    @property
    def has_bounding_boxes(self) -> bool:
        return self._bboxes is not None

    @property
    def num_elements(self) -> int:
        return self._volumes.shape[0]

    def __iter__(self):
        for mat in np.unique(self._materials):
            if mat > 0:
                yield mat

    def __len__(self) -> int:
        return (np.unique(self._materials) > 0).sum()

    def __repr__(self) -> str:
        ids, counts = np.unique(self._materials, return_counts=True)
        return '{' + '\n '.join(
            f'{id}: <{count} nonzero volumes>' for id, count in zip(ids, counts) if id > 0) + '}'

    def __getitem__(self, material_id: int) -> np.ndarray:
        volumes = np.zeros(self.num_elements)
        for i in range(self._volumes.shape[1]):
            indices = (self._materials[:, i] == material_id)
            volumes[indices] = self._volumes[indices, i]
        return volumes

    def by_element(
        self,
        index_elem: int,
        include_bboxes: bool = False
    ) -> list[tuple[int | None, float] | tuple[int | None, float, BoundingBox | None]]:
        """Get a list of volumes for each material within a specific element.

        Parameters
        ----------
        index_elem : int
            Mesh element index

        Returns
        -------
        list of tuple
            If ``include_bboxes`` is False (default), returns tuples of
            (material ID, volume). If ``include_bboxes`` is True, returns
            tuples of (material ID, volume, bounding box).

        """
        table_size = self._volumes.shape[1]
        if include_bboxes and self._bboxes is None:
            raise ValueError('Bounding boxes were not computed for this object.')

        results = []
        for i in range(table_size):
            m = self._materials[index_elem, i]
            if m == -2:
                continue
            mat_id = m if m > -1 else None
            vol = self._volumes[index_elem, i]

            if include_bboxes:
                vals = self._bboxes[index_elem, i]
                bbox = BoundingBox(vals[0:3], vals[3:6])
                results.append((mat_id, vol, bbox))
            else:
                results.append((mat_id, vol))

        return results

    def save(self, filename: PathLike):
        """Save material volumes to a .npz file.

        Parameters
        ----------
        filename : path-like
            Filename where data will be saved
        """
        kwargs = {'materials': self._materials, 'volumes': self._volumes}
        if self._bboxes is not None:
            kwargs['bboxes'] = self._bboxes
        np.savez_compressed(filename, **kwargs)

    @classmethod
    def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes:
        """Generate material volumes from a .npz file

        Parameters
        ----------
        filename : path-like
            File where data will be read from

        """
        filedata = np.load(filename)
        bboxes = filedata['bboxes'] if 'bboxes' in filedata.files else None
        return cls(filedata['materials'], filedata['volumes'], bboxes)


class MeshBase(IDManagerMixin, ABC):
    """A mesh that partitions geometry for tallying purposes.

    Parameters
    ----------
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    lower_left : Iterable of float
        The lower-left coordinates
    upper_right : Iterable of float
        The upper-right coordinates
    bounding_box : openmc.BoundingBox
        Axis-aligned bounding box of the mesh as defined by the upper-right and
        lower-left coordinates.
    indices : Iterable of tuple
        An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...]
    n_elements : int
        Number of elements in the mesh
    """

    next_id = 1
    used_ids = set()

    def __init__(self, mesh_id: int | None = None, name: str = ''):
        # Initialize Mesh class attributes
        self.id = mesh_id
        self.name = name

    @property
    def name(self):
        return self._name

    @name.setter
    def name(self, name: str):
        if name is not None:
            cv.check_type(f'name for mesh ID="{self._id}"', name, str)
            self._name = name
        else:
            self._name = ''
            
    @property
    @abstractmethod
    def lower_left(self):
        pass
        
    @property
    @abstractmethod
    def upper_right(self):
        pass

    @property
    def bounding_box(self) -> openmc.BoundingBox:
        return openmc.BoundingBox(self.lower_left, self.upper_right)

    @property
    @abstractmethod
    def indices(self):
        pass
        
    @property
    @abstractmethod
    def n_elements(self):
        pass

    def __repr__(self):
        string = type(self).__name__ + '\n'
        string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
        string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
        return string

    def _volume_dim_check(self):
        if self.n_dimension != 3 or \
           any([d == 0 for d in self.dimension]):
            raise RuntimeError(f'Mesh {self.id} is not 3D. '
                               'Volumes cannot be provided.')

    @classmethod
    def from_hdf5(cls, group: h5py.Group):
        """Create mesh from HDF5 group

        Parameters
        ----------
        group : h5py.Group
            Group in HDF5 file

        Returns
        -------
        openmc.MeshBase
            Instance of a MeshBase subclass

        """
        mesh_type = 'regular' if 'type' not in group.keys() else group['type'][()].decode()
        mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
        mesh_name = '' if not 'name' in group else group['name'][()].decode()

        if mesh_type == 'regular':
            return RegularMesh.from_hdf5(group, mesh_id, mesh_name)
        elif mesh_type == 'rectilinear':
            return RectilinearMesh.from_hdf5(group, mesh_id, mesh_name)
        elif mesh_type == 'cylindrical':
            return CylindricalMesh.from_hdf5(group, mesh_id, mesh_name)
        elif mesh_type == 'spherical':
            return SphericalMesh.from_hdf5(group, mesh_id, mesh_name)
        elif mesh_type == 'unstructured':
            return UnstructuredMesh.from_hdf5(group, mesh_id, mesh_name)
        else:
            raise ValueError('Unrecognized mesh type: "' + mesh_type + '"')

    def to_xml_element(self):
        """Return XML representation of the mesh

        Returns
        -------
        element : lxml.etree._Element
            XML element containing mesh data

        """
        elem = ET.Element("mesh")

        elem.set("id", str(self._id))
        if self.name:
            elem.set("name", self.name)

        return elem

    @classmethod
    def from_xml_element(cls, elem: ET.Element):
        """Generates a mesh from an XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element

        Returns
        -------
        openmc.MeshBase
            an openmc mesh object

        """
        mesh_type = get_text(elem, 'type')

        if mesh_type == 'regular' or mesh_type is None:
            mesh = RegularMesh.from_xml_element(elem)
        elif mesh_type == 'rectilinear':
            mesh = RectilinearMesh.from_xml_element(elem)
        elif mesh_type == 'cylindrical':
            mesh = CylindricalMesh.from_xml_element(elem)
        elif mesh_type == 'spherical':
            mesh = SphericalMesh.from_xml_element(elem)
        elif mesh_type == 'unstructured':
            mesh = UnstructuredMesh.from_xml_element(elem)
        else:
            raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.')

        mesh.name = get_text(elem, 'name', default='')
        return mesh

    def get_homogenized_materials(
            self,
            model: openmc.Model,
            n_samples: int | tuple[int, int, int] = 10_000,
            include_void: bool = True,
            material_volumes: MeshMaterialVolumes | None = None,
            **kwargs
    ) -> list[openmc.Material]:
        """Generate homogenized materials over each element in a mesh.

        .. versionadded:: 0.15.0

        Parameters
        ----------
        model : openmc.Model
            Model containing materials to be homogenized and the associated
            geometry.
        n_samples : int or 2-tuple of int
            Total number of rays to sample. The number of rays in each direction
            is determined by the aspect ratio of the mesh bounding box. When
            specified as a 3-tuple, it is interpreted as the number of rays in
            the x, y, and z dimensions.
        include_void : bool, optional
            Whether homogenization should include voids.
        material_volumes : MeshMaterialVolumes, optional
            Previously computed mesh material volumes to use for homogenization.
            If not provided, they will be computed by calling
            :meth:`material_volumes`.
        **kwargs
            Keyword-arguments passed to :meth:`material_volumes`.

        Returns
        -------
        list of openmc.Material
            Homogenized material in each mesh element

        """
        if material_volumes is None:
            vols = self.material_volumes(model, n_samples, **kwargs)
        else:
            vols = material_volumes
        mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]

        # Get dictionary of all materials
        materials = model._get_all_materials()

        # Create homogenized material for each element
        homogenized_materials = []
        for mat_volume_list in mat_volume_by_element:
            material_ids, volumes = [list(x) for x in zip(*mat_volume_list)]
            total_volume = sum(volumes)

            # Check for void material and remove
            try:
                index_void = material_ids.index(None)
            except ValueError:
                pass
            else:
                material_ids.pop(index_void)
                volumes.pop(index_void)

            # If void should be excluded, adjust total volume
            if not include_void:
                total_volume = sum(volumes)

            # Compute volume fractions
            volume_fracs = np.array(volumes) / total_volume

            # Get list of materials and mix 'em up!
            mats = [materials[uid] for uid in material_ids]
            homogenized_mat = openmc.Material.mix_materials(
                mats, volume_fracs, 'vo'
            )
            homogenized_mat.volume = total_volume
            homogenized_materials.append(homogenized_mat)

        return homogenized_materials

    def material_volumes(
            self,
            model: openmc.Model,
            n_samples: int | tuple[int, int, int] = 10_000,
            max_materials: int = 4,
            bounding_boxes: bool = False,
            **kwargs
    ) -> MeshMaterialVolumes:
        """Determine volume of materials in each mesh element.

        This method works by raytracing repeatedly through the mesh to count the
        estimated volume of each material in all mesh elements. Three sets of
        rays are used: one set parallel to the x-axis, one parallel to the
        y-axis, and one parallel to the z-axis.

        .. versionadded:: 0.15.1

        Parameters
        ----------
        model : openmc.Model
            Model containing materials.
        n_samples : int or 3-tuple of int
            Total number of rays to sample. The number of rays in each direction
            is determined by the aspect ratio of the mesh bounding box. When
            specified as a 3-tuple, it is interpreted as the number of rays in
            the x, y, and z dimensions.
        max_materials : int, optional
            Estimated maximum number of materials in any given mesh element.
        bounding_boxes : bool, optional
            Whether to compute an axis-aligned bounding box for each
            (mesh element, material) combination. When enabled, the bounding
            box encloses the ray-estimator prisms used for the volume
            estimation.
        **kwargs : dict
            Keyword arguments passed to :func:`openmc.lib.init`

        Returns
        -------
        Dictionary-like object that maps material IDs to an array of volumes
        equal in size to the number of mesh elements.

        """
        import openmc.lib

        # In order to get mesh into model, we temporarily replace the
        # tallies with a single mesh tally using the current mesh
        original_tallies = list(model.tallies)
        new_tally = openmc.Tally()
        new_tally.filters = [openmc.MeshFilter(self)]
        new_tally.scores = ['flux']
        model.tallies = [new_tally]

        # Set default arguments
        kwargs.setdefault('output', True)
        if 'args' in kwargs:
            kwargs['args'] = ['-c'] + kwargs['args']
        kwargs.setdefault('args', ['-c'])

        with openmc.lib.TemporarySession(model, **kwargs):
            # Get mesh from single tally
            mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh

            # Compute material volumes
            volumes = mesh.material_volumes(
                n_samples, max_materials, output=kwargs['output'],
                bounding_boxes=bounding_boxes)

        # Restore original tallies
        model.tallies = original_tallies
        return volumes


class StructuredMesh(MeshBase):
    """A base class for structured mesh functionality

    Parameters
    ----------
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    @property
    @abstractmethod
    def dimension(self):
        pass

    @property
    @abstractmethod
    def n_dimension(self):
        pass

    @property
    @abstractmethod
    def _grids(self):
        pass

    @property
    def vertices(self):
        """Return coordinates of mesh vertices in Cartesian coordinates. Also
           see :meth:`CylindricalMesh.vertices_cylindrical` and
           :meth:`SphericalMesh.vertices_spherical` for coordinates in other coordinate
           systems.

        Returns
        -------
        vertices : numpy.ndarray
            Returns a numpy.ndarray representing the coordinates of the mesh
            vertices with a shape equal to (dim1 + 1, ..., dimn + 1, ndim). X, Y, Z values
            can be unpacked with xx, yy, zz = np.rollaxis(mesh.vertices, -1).

        """
        return self._generate_vertices(*self._grids)

    @staticmethod
    def _generate_vertices(i_grid, j_grid, k_grid):
        """Returns an array with shape (i_grid.size, j_grid.size, k_grid.size, 3)
           containing the corner vertices of mesh elements.
        """
        return np.stack(np.meshgrid(i_grid, j_grid, k_grid, indexing='ij'), axis=-1)

    @staticmethod
    def _generate_edge_midpoints(grids):
        """Generates the midpoints of mesh element edges for each dimension of the mesh.

        Parameters
        ----------
        grids : numpy.ndarray
            The vertex grids along each dimension of the mesh.

        Returns
        -------
        midpoint_grids : list of numpy.ndarray
            The edge midpoints for the i, j, and k midpoints of each element in
            i, j, k ordering. The shapes of the resulting grids are
            [(ni-1, nj, nk, 3), (ni, nj-1, nk, 3), (ni, nj, nk-1, 3)]
        """
        # generate a set of edge midpoints for each dimension
        midpoint_grids = []
        # generate the element edge midpoints in order s.t.
        # the epxected element ordering is preserved with respect to the corner vertices

        # each grid is comprised of the mid points for one dimension and the
        # corner vertices of the other two
        for dims in ((0, 1, 2), (1, 0, 2), (2, 0, 1)):
            # compute the midpoints along the last dimension
            midpoints = grids[dims[0]][:-1] + 0.5 * np.diff(grids[dims[0]])

            coords = (midpoints, grids[dims[1]], grids[dims[2]])

            i_grid, j_grid, k_grid = [coords[dims.index(i)] for i in range(3)]

            # re-use the generate vertices method to create the full mesh grid
            # transpose to get (i, j, k) ordering of the gridpoints
            midpoint_grid = StructuredMesh._generate_vertices(i_grid, j_grid, k_grid)
            midpoint_grids.append(midpoint_grid)

        return midpoint_grids

    @property
    def midpoint_vertices(self):
        """Create vertices that lie on the midpoint of element edges
        """
        # generate edge midpoints needed for curvilinear element definition
        midpoint_vertices = self._generate_edge_midpoints(self._grids)

        # convert each of the midpoint grids to cartesian coordinates
        for vertices in midpoint_vertices:
            self._convert_to_cartesian(vertices, self.origin)

        return midpoint_vertices

    @property
    def centroids(self):
        """Return coordinates of mesh element centroids.

        Returns
        -------
        centroids : numpy.ndarray
            Returns a numpy.ndarray representing the mesh element centroid
            coordinates with a shape equal to (dim1, ..., dimn, ndim). X,
            Y, Z values can be unpacked with xx, yy, zz =
            np.rollaxis(mesh.centroids, -1).
        """
        ndim = self.n_dimension
        # this line ensures that the vertices aren't adjusted by the origin or
        # converted to the Cartesian system for cylindrical and spherical meshes
        vertices = StructuredMesh.vertices.fget(self)
        s0 = (slice(0, -1),)*ndim + (slice(None),)
        s1 = (slice(1, None),)*ndim + (slice(None),)
        return (vertices[s0] + vertices[s1]) / 2
    
    @property
    def n_elements(self):
        return np.prod(self.dimension)

    @property
    def num_mesh_cells(self):
        warnings.warn(
            "The 'num_mesh_cells' attribute is deprecated and will be removed in a future version. "
            "Use 'n_elements' instead.",
            FutureWarning, stacklevel=2
        )
        return self.n_elements

    def write_data_to_vtk(self,
                          filename: PathLike,
                          datasets: dict | None = None,
                          volume_normalization: bool = True,
                          curvilinear: bool = False):
        """Creates a VTK object of the mesh

        Parameters
        ----------
        filename : str
            Name of the VTK file to write.
        datasets : dict
            Dictionary whose keys are the data labels and values are the data
            sets. 1D datasets are expected to be extracted directly from
            statepoint data without reordering/reshaping. Multidimensional
            datasets are expected to have the same dimensions as the mesh itself
            with structured indexing in "C" ordering. See the "expand_dims" flag
            of :meth:`~openmc.Tally.get_reshaped_data` on reshaping tally data when using
            :class:`~openmc.MeshFilter`'s.
        volume_normalization : bool, optional
            Whether or not to normalize the data by the volume of the mesh
            elements.
        curvilinear : bool
            Whether or not to write curvilinear elements. Only applies to
            ``SphericalMesh`` and ``CylindricalMesh``.

        Raises
        ------
        ValueError
            When the size of a dataset doesn't match the number of mesh cells

        Returns
        -------
        vtk.StructuredGrid or vtk.UnstructuredGrid
            a VTK grid object representing the mesh

        Examples
        --------
        1D data from a tally with only a mesh filter and heating score:

            # pass the tally mean property of shape (N, 1, 1) directly to this
            # method; dimensions of size 1 will automatically removed
            >>> heating = tally.mean
            >>> mesh.write_data_to_vtk({'heating': heating})

        Multidimensional data from a tally with only a mesh

           # retrieve a data array with the mesh filter expanded into three
           # dimensions, ijk; additional dimensions of size one will
           # automatically be removed
           >>> heating = tally.get_reshaped_data(expand_dims=True)
           >>> mesh.write_data_to_vtk({'heating': heating})
        """
        import vtk
        from vtk.util import numpy_support as nps

        # write linear elements using a structured grid
        if not curvilinear or isinstance(self, (RegularMesh, RectilinearMesh)):
            vtk_grid = self._create_vtk_structured_grid()
            writer = vtk.vtkStructuredGridWriter()
        # write curvilinear elements using an unstructured grid
        else:
            vtk_grid = self._create_vtk_unstructured_grid()
            writer = vtk.vtkUnstructuredGridWriter()

        if datasets is not None:
            # maintain a list of the datasets as added to the VTK arrays to
            # ensure they persist in memory until the file is written
            datasets_out = []
            for label, dataset in datasets.items():
                dataset = self._reshape_vtk_dataset(dataset)
                self._check_vtk_dataset(label, dataset)
                # If the array data is 3D, assume is in C ordering and transpose
                # before flattening to match the ordering expected by the VTK
                # array based on the way mesh indices are ordered in the Python
                # API
                # TODO: update to "C" ordering throughout
                if dataset.ndim == 3:
                    dataset = dataset.T.ravel()
                datasets_out.append(dataset)

                if volume_normalization:
                    dataset /= self.volumes.T.ravel()

                dataset_array = vtk.vtkDoubleArray()
                dataset_array.SetName(label)
                dataset_array.SetArray(nps.numpy_to_vtk(dataset), dataset.size, True)
                vtk_grid.GetCellData().AddArray(dataset_array)

        writer.SetFileName(str(filename))
        writer.SetInputData(vtk_grid)
        writer.Write()

        return vtk_grid

    def _create_vtk_structured_grid(self):
        """Create a structured grid

        Returns
        -------
        vtk.vtkStructuredGrid
            a VTK structured grid object representing the mesh
        """
        import vtk
        from vtk.util import numpy_support as nps

        vtkPts = vtk.vtkPoints()
        vtkPts.SetData(nps.numpy_to_vtk(np.swapaxes(self.vertices, 0, 2).reshape(-1, 3), deep=True))
        vtk_grid = vtk.vtkStructuredGrid()
        vtk_grid.SetPoints(vtkPts)
        vtk_grid.SetDimensions(*[dim + 1 for dim in self.dimension])

        return vtk_grid

    def _create_vtk_unstructured_grid(self):
        """Create an unstructured grid of curvilinear elements
           representing the mesh

        Returns
        -------
        vtk.vtkUnstructuredGrid
            a VTK unstructured grid object representing the mesh
        """
        import vtk
        from vtk.util import numpy_support as nps

        corner_vertices = np.swapaxes(self.vertices, 0, 2).reshape(-1, 3)

        vtkPts = vtk.vtkPoints()
        vtk_grid = vtk.vtkUnstructuredGrid()
        vtk_grid.SetPoints(vtkPts)
        # add corner vertices to the point set for the unstructured grid
        # only insert unique points, we'll get their IDs in the point set to
        # define element connectivity later
        vtkPts.SetData(nps.numpy_to_vtk(np.unique(corner_vertices, axis=0), deep=True))

        # create a locator to assist with duplicate points
        locator = vtk.vtkPointLocator()
        locator.SetDataSet(vtk_grid)
        locator.AutomaticOn() # autmoatically adds points to locator
        locator.InitPointInsertion(vtkPts, vtkPts.GetBounds())
        locator.BuildLocator()

        # this function is used to add new points to the unstructured
        # grid. It will return an existing point ID if the point is alread present
        def _insert_point(pnt):
            result = locator.IsInsertedPoint(pnt)
            if result == -1:
                point_id = vtkPts.InsertNextPoint(pnt)
                locator.InsertPoint(point_id, pnt)
                return point_id
            else:
                return result

        # Add all points to the unstructured grid, maintaining a flat list of IDs as we go ###

        # flat array storing point IDs for a given vertex
        # in the grid
        point_ids = []

        # add element corner vertices to array
        for pnt in corner_vertices:
            point_ids.append(_insert_point(pnt))

        # get edge midpoints and add them to the
        # list of point IDs
        midpoint_vertices = self.midpoint_vertices
        for edge_grid in midpoint_vertices:
            for pnt in np.swapaxes(edge_grid, 0, 2).reshape(-1, 3):
                point_ids.append(_insert_point(pnt))

        # determine how many elements in each dimension
        # and how many points in each grid
        n_elem = np.asarray(self.dimension)
        n_pnts = n_elem + 1

        # create hexes and set points for corner
        # vertices
        for i, j, k in self.indices:
            # handle indices indexed from one
            i -= 1
            j -= 1
            k -= 1

            # create a new vtk hex
            hex = vtk.vtkQuadraticHexahedron()

            # set connectivity the hex corners
            for n, (di, dj, dk) in enumerate(_HEX_VERTEX_CONN):
                # compute flat index into the point ID list based on i, j, k
                # of the vertex
                flat_idx = np.ravel_multi_index((i+di, j+dj, k+dk), n_pnts, order='F')
                # set corner vertices
                hex.GetPointIds().SetId(n, point_ids[flat_idx])

            # set connectivity of the hex midpoints
            n_midpoint_vertices = [v.size // 3 for v in midpoint_vertices]
            for n, (dim, (di, dj, dk)) in enumerate(_HEX_MIDPOINT_CONN):
                # initial offset for corner vertices and midpoint dimension
                flat_idx = corner_vertices.shape[0] + sum(n_midpoint_vertices[:dim])
                # generate a flat index into the table of point IDs
                midpoint_shape = midpoint_vertices[dim].shape[:-1]
                flat_idx += np.ravel_multi_index((i+di, j+dj, k+dk),
                                                 midpoint_shape,
                                                 order='F')
                # set hex midpoint connectivity
                hex.GetPointIds().SetId(_N_HEX_VERTICES + n, point_ids[flat_idx])

            # add the hex to the grid
            vtk_grid.InsertNextCell(hex.GetCellType(), hex.GetPointIds())

        return vtk_grid

    @staticmethod
    def _reshape_vtk_dataset(dataset):
        """Reshape a dataset to be compatible with VTK output

        This method performs the following operations on a dataset:
        1. Convert to numpy array if not already
        2. Remove any trailing dimensions of size 1
        3. Squeeze out any extra dimensions of size 1 beyond the first 3

        Parameters
        ----------
        dataset : array-like
            The dataset to reshape

        Returns
        -------
        numpy.ndarray
            The reshaped dataset
        """
        reshaped_data = np.asarray(dataset)

        # detect flat array with extra dims
        if all(d == 1 for d in reshaped_data.shape[1:]):
            reshaped_data = reshaped_data.squeeze()

        # remove any higher dimensions with size 1
        if reshaped_data.ndim > 3 and all(d == 1 for d in reshaped_data.shape[3:]):
            reshaped_data = reshaped_data.reshape(reshaped_data.shape[:3])

        if np.shares_memory(reshaped_data, dataset):
            return np.copy(reshaped_data)
        else:
            return reshaped_data

    def _check_vtk_dataset(self, label: str, dataset: np.ndarray):
        """Perform some basic checks that a dataset is valid for this Mesh

        Parameters
        ----------
        label : str
            The label for the dataset being checked
        dataset : numpy.ndarray
            The dataset array to check against this mesh's dimensions

        """
        cv.check_type('data label', label, str)

        if dataset.size != self.n_elements:
            raise ValueError(
                f"The size of the dataset '{label}' ({dataset.size}) should be"
                f" equal to the number of mesh cells ({self.n_elements})"
            )

        # accept a flat array as-is, assuming it is in the correct order
        if dataset.ndim == 1:
            return

        if dataset.shape != self.dimension:
            raise ValueError(
                f'Cannot apply multidimensional dataset "{label}" with '
                f"shape {dataset.shape} to mesh {self.id} "
                f"with dimensions {self.dimension}"
            )


class HasBoundingBox(Protocol):
    """Object that has a ``bounding_box`` attribute."""
    bounding_box: openmc.BoundingBox


class RegularMesh(StructuredMesh):
    """A regular Cartesian mesh in one, two, or three dimensions

    Parameters
    ----------
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    dimension : Iterable of int
        The number of mesh cells in each direction (x, y, z).
    n_dimension : int
        Number of mesh dimensions.
    lower_left : Iterable of float
        The lower-left corner of the structured mesh. If only two coordinate
        are given, it is assumed that the mesh is an x-y mesh.
    upper_right : Iterable of float
        The upper-right corner of the structured mesh. If only two coordinate
        are given, it is assumed that the mesh is an x-y mesh.
    bounding_box : openmc.BoundingBox
        Axis-aligned bounding box of the mesh as defined by the upper-right and
        lower-left coordinates.
    width : Iterable of float
        The width of mesh cells in each direction.
    indices : Iterable of tuple
        An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
        (2, 1, 1), ...]

    """

    def __init__(self, mesh_id: int | None = None, name: str = ''):
        super().__init__(mesh_id, name)

        self._dimension = None
        self._lower_left = None
        self._upper_right = None
        self._width = None

    @property
    def dimension(self):
        return tuple(self._dimension)

    @dimension.setter
    def dimension(self, dimension: Iterable[int]):
        cv.check_type('mesh dimension', dimension, Iterable, Integral)
        cv.check_length('mesh dimension', dimension, 1, 3)
        self._dimension = dimension

    @property
    def n_dimension(self):
        if self._dimension is not None:
            return len(self._dimension)
        else:
            return None

    @property
    def lower_left(self):
        return self._lower_left

    @lower_left.setter
    def lower_left(self, lower_left: Iterable[Real]):
        cv.check_type('mesh lower_left', lower_left, Iterable, Real)
        cv.check_length('mesh lower_left', lower_left, 1, 3)
        self._lower_left = lower_left

        if self.upper_right is not None and any(np.isclose(self.upper_right, lower_left)):
            raise ValueError("Mesh cannot have zero thickness in any dimension")

    @property
    def upper_right(self):
        if self._upper_right is not None:
            return self._upper_right
        elif self._width is not None:
            if self._lower_left is not None and self._dimension is not None:
                ls = self._lower_left
                ws = self._width
                dims = self._dimension
                return [l + w * d for l, w, d in zip(ls, ws, dims)]

    @upper_right.setter
    def upper_right(self, upper_right: Iterable[Real]):
        cv.check_type('mesh upper_right', upper_right, Iterable, Real)
        cv.check_length('mesh upper_right', upper_right, 1, 3)
        self._upper_right = upper_right

        if self._width is not None:
            self._width = None
            warnings.warn("Unsetting width attribute.")

        if self.lower_left is not None and any(np.isclose(self.lower_left, upper_right)):
            raise ValueError("Mesh cannot have zero thickness in any dimension")

    @property
    def width(self):
        if self._width is not None:
            return self._width
        elif self._upper_right is not None:
            if self._lower_left is not None and self._dimension is not None:
                us = self._upper_right
                ls = self._lower_left
                dims =  self._dimension
                return [(u - l) / d for u, l, d in zip(us, ls, dims)]

    @width.setter
    def width(self, width: Iterable[Real]):
        cv.check_type('mesh width', width, Iterable, Real)
        cv.check_length('mesh width', width, 1, 3)
        self._width = width

        if self._upper_right is not None:
            self._upper_right = None
            warnings.warn("Unsetting upper_right attribute.")

    @property
    def volumes(self):
        """Return Volumes for every mesh cell

        Returns
        -------
        volumes : numpy.ndarray
            Volumes

        """
        self._volume_dim_check()
        return np.full(self.dimension, np.prod(self.width))

    @property
    def total_volume(self):
        return np.prod(self.dimension) * np.prod(self.width)

    @property
    def indices(self):
        ndim = len(self._dimension)
        if ndim == 3:
            nx, ny, nz = self.dimension
            return ((x, y, z)
                    for z in range(1, nz + 1)
                    for y in range(1, ny + 1)
                    for x in range(1, nx + 1))
        elif ndim == 2:
            nx, ny = self.dimension
            return ((x, y)
                    for y in range(1, ny + 1)
                    for x in range(1, nx + 1))
        else:
            nx, = self.dimension
            return ((x,) for x in range(1, nx + 1))

    @property
    def _grids(self):
        ndim = len(self._dimension)
        if ndim == 3:
            x0, y0, z0 = self.lower_left
            x1, y1, z1 = self.upper_right
            nx, ny, nz = self.dimension
            xarr = np.linspace(x0, x1, nx + 1)
            yarr = np.linspace(y0, y1, ny + 1)
            zarr = np.linspace(z0, z1, nz + 1)
            return (xarr, yarr, zarr)
        elif ndim == 2:
            x0, y0 = self.lower_left
            x1, y1 = self.upper_right
            nx, ny = self.dimension
            xarr = np.linspace(x0, x1, nx + 1)
            yarr = np.linspace(y0, y1, ny + 1)
            return (xarr, yarr)
        else:
            nx, = self.dimension
            x0, = self.lower_left
            x1, = self.upper_right
            return (np.linspace(x0, x1, nx + 1),)

    def __repr__(self):
        string = super().__repr__()
        string += '{0: <16}{1}{2}\n'.format('\tDimensions', '=\t', self.n_dimension)
        string += '{0: <16}{1}{2}\n'.format('\tVoxels', '=\t', self._dimension)
        string += '{0: <16}{1}{2}\n'.format('\tLower left', '=\t', self._lower_left)
        string += '{0: <16}{1}{2}\n'.format('\tUpper Right', '=\t', self.upper_right)
        string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self.width)
        return string

    @classmethod
    def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
        # Read and assign mesh properties
        mesh = cls(mesh_id=mesh_id, name=name)
        mesh.dimension = group['dimension'][()]
        mesh.lower_left = group['lower_left'][()]
        if 'width' in group:
            mesh.width = group['width'][()]
        elif 'upper_right' in group:
            mesh.upper_right = group['upper_right'][()]
        else:
            raise IOError('Invalid mesh: must have one of "upper_right" or "width"')

        return mesh

    @classmethod
    def from_rect_lattice(
        cls,
        lattice: 'openmc.RectLattice',
        division: int = 1,
        mesh_id: int | None = None,
        name: str = ''
    ):
        """Create mesh from an existing rectangular lattice

        Parameters
        ----------
        lattice : openmc.RectLattice
            Rectangular lattice used as a template for this mesh
        division : int
            Number of mesh cells per lattice cell.
            If not specified, there will be 1 mesh cell per lattice cell.
        mesh_id : int
            Unique identifier for the mesh
        name : str
            Name of the mesh

        Returns
        -------
        openmc.RegularMesh
            RegularMesh instance

        """
        cv.check_type('rectangular lattice', lattice, openmc.RectLattice)

        shape = np.array(lattice.shape)
        width = lattice.pitch*shape

        mesh = cls(mesh_id=mesh_id, name=name)
        mesh.lower_left = lattice.lower_left
        mesh.upper_right = lattice.lower_left + width
        mesh.dimension = shape*division

        return mesh

    @classmethod
    def from_domain(
        cls,
        domain: HasBoundingBox | BoundingBox,
        dimension: Sequence[int] | int = 1000,
        mesh_id: int | None = None,
        name: str = ''
    ):
        """Create RegularMesh from a domain using its bounding box.

        Parameters
        ----------
        domain : HasBoundingBox | openmc.BoundingBox
            The object passed in will be used as a template for this mesh. The
            bounding box of the property of the object passed will be used to
            set the lower_left and upper_right and of the mesh instance.
            Alternatively, a :class:`openmc.BoundingBox` can be passed
            directly.
        dimension : Iterable of int | int
            The number of mesh cells in total or number of mesh cells in each
            direction (x, y, z). If a single integer is provided, the domain
            will will be divided into that many mesh cells with roughly equal
            lengths in each direction (cubes).
        mesh_id : int
            Unique identifier for the mesh
        name : str
            Name of the mesh

        Returns
        -------
        openmc.RegularMesh
            RegularMesh instance

        """
        if isinstance(domain, BoundingBox):
            bb = domain
        elif hasattr(domain, 'bounding_box'):
            bb = domain.bounding_box
        else:
            raise TypeError("Domain must be a BoundingBox or have a "
                            "bounding_box property")

        mesh = cls(mesh_id=mesh_id, name=name)
        mesh.lower_left = bb[0]
        mesh.upper_right = bb[1]
        if isinstance(dimension, int):
            cv.check_greater_than("dimension", dimension, 1, equality=True)
            # If a single integer is provided, divide the domain into that many
            # mesh cells with roughly equal lengths in each direction
            ideal_cube_volume = bb.volume / dimension
            ideal_cube_size = ideal_cube_volume ** (1 / 3)
            dimension = [
                max(1, int(round(side / ideal_cube_size)))
                for side in bb.width
            ]
        mesh.dimension = dimension

        return mesh

    def to_xml_element(self):
        """Return XML representation of the mesh

        Returns
        -------
        element : lxml.etree._Element
            XML element containing mesh data

        """
        element = super().to_xml_element()

        if self._dimension is not None:
            subelement = ET.SubElement(element, "dimension")
            subelement.text = ' '.join(map(str, self._dimension))

        subelement = ET.SubElement(element, "lower_left")
        subelement.text = ' '.join(map(str, self._lower_left))

        if self._upper_right is not None:
            subelement = ET.SubElement(element, "upper_right")
            subelement.text = ' '.join(map(str, self._upper_right))
        if self._width is not None:
            subelement = ET.SubElement(element, "width")
            subelement.text = ' '.join(map(str, self._width))

        return element

    @classmethod
    def from_xml_element(cls, elem: ET.Element):
        """Generate mesh from an XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element

        Returns
        -------
        openmc.Mesh
            Mesh generated from XML element

        """
        mesh_id = int(get_text(elem, 'id'))
        mesh = cls(mesh_id=mesh_id)

        dimension = get_elem_list(elem, "dimension", int)
        if dimension is not None:
            mesh.dimension = dimension

        lower_left = get_elem_list(elem, "lower_left", float)
        if lower_left is not None:
            mesh.lower_left = lower_left

        upper_right = get_elem_list(elem, "upper_right", float)
        if upper_right is not None:
            mesh.upper_right = upper_right

        width = get_elem_list(elem, "width", float)
        if width is not None:
            mesh.width = width

        return mesh

    def build_cells(self, bc: str | None = None):
        """Generates a lattice of universes with the same dimensionality
        as the mesh object.  The individual cells/universes produced
        will not have material definitions applied and so downstream code
        will have to apply that information.

        Parameters
        ----------
        bc : iterable of {'reflective', 'periodic', 'transmission', 'vacuum', or 'white'}
            Boundary conditions for each of the four faces of a rectangle
            (if applying to a 2D mesh) or six faces of a parallelepiped
            (if applying to a 3D mesh) provided in the following order:
            [x min, x max, y min, y max, z min, z max].  2-D cells do not
            contain the z min and z max entries. Defaults to 'reflective' for
            all faces.

        Returns
        -------
        root_cell : openmc.Cell
            The cell containing the lattice representing the mesh geometry;
            this cell is a single parallelepiped with boundaries matching
            the outermost mesh boundary with the boundary conditions from bc
            applied.
        cells : iterable of openmc.Cell
            The list of cells within each lattice position mimicking the mesh
            geometry.

        """
        if bc is None:
            bc = ['reflective'] * 6
        if len(bc) not in (4, 6):
            raise ValueError('Boundary condition must be of length 4 or 6')
        for entry in bc:
            cv.check_value('bc', entry, _BOUNDARY_TYPES)

        n_dim = self.n_dimension

        # Build the cell which will contain the lattice
        xplanes = [openmc.XPlane(self.lower_left[0], boundary_type=bc[0]),
                   openmc.XPlane(self.upper_right[0], boundary_type=bc[1])]
        if n_dim == 1:
            yplanes = [openmc.YPlane(-1e10, boundary_type='reflective'),
                       openmc.YPlane(1e10, boundary_type='reflective')]
        else:
            yplanes = [openmc.YPlane(self.lower_left[1], boundary_type=bc[2]),
                       openmc.YPlane(self.upper_right[1], boundary_type=bc[3])]

        if n_dim <= 2:
            # Would prefer to have the z ranges be the max supported float, but
            # these values are apparently different between python and Fortran.
            # Choosing a safe and sane default.
            # Values of +/-1e10 are used here as there seems to be an
            # inconsistency between what numpy uses as the max float and what
            # Fortran expects for a real(8), so this avoids code complication
            # and achieves the same goal.
            zplanes = [openmc.ZPlane(-1e10, boundary_type='reflective'),
                       openmc.ZPlane(1e10, boundary_type='reflective')]
        else:
            zplanes = [openmc.ZPlane(self.lower_left[2], boundary_type=bc[4]),
                       openmc.ZPlane(self.upper_right[2], boundary_type=bc[5])]
        root_cell = openmc.Cell()
        root_cell.region = ((+xplanes[0] & -xplanes[1]) &
                            (+yplanes[0] & -yplanes[1]) &
                            (+zplanes[0] & -zplanes[1]))

        # Build the universes which will be used for each of the (i,j,k)
        # locations within the mesh.
        # We will concurrently build cells to assign to these universes
        cells = []
        universes = []
        for _ in self.indices:
            cells.append(openmc.Cell())
            universes.append(openmc.Universe())
            universes[-1].add_cell(cells[-1])

        lattice = openmc.RectLattice()
        lattice.lower_left = self.lower_left

        # Assign the universe and rotate to match the indexing expected for
        # the lattice
        if n_dim == 1:
            universe_array = np.array([universes])
        elif n_dim == 2:
            universe_array = np.empty(self.dimension[::-1],
                                      dtype=openmc.Universe)
            i = 0
            for y in range(self.dimension[1] - 1, -1, -1):
                for x in range(self.dimension[0]):
                    universe_array[y][x] = universes[i]
                    i += 1
        else:
            universe_array = np.empty(self.dimension[::-1],
                                      dtype=openmc.Universe)
            i = 0
            for z in range(self.dimension[2]):
                for y in range(self.dimension[1] - 1, -1, -1):
                    for x in range(self.dimension[0]):
                        universe_array[z][y][x] = universes[i]
                        i += 1
        lattice.universes = universe_array

        if self.width is not None:
            lattice.pitch = self.width
        else:
            dx = ((self.upper_right[0] - self.lower_left[0]) /
                  self.dimension[0])

            if n_dim == 1:
                lattice.pitch = [dx]
            elif n_dim == 2:
                dy = ((self.upper_right[1] - self.lower_left[1]) /
                      self.dimension[1])
                lattice.pitch = [dx, dy]
            else:
                dy = ((self.upper_right[1] - self.lower_left[1]) /
                      self.dimension[1])
                dz = ((self.upper_right[2] - self.lower_left[2]) /
                      self.dimension[2])
                lattice.pitch = [dx, dy, dz]

        # Fill Cell with the Lattice
        root_cell.fill = lattice

        return root_cell, cells


def Mesh(*args, **kwargs):
    warnings.warn("Mesh has been renamed RegularMesh. Future versions of "
                  "OpenMC will not accept the name Mesh.")
    return RegularMesh(*args, **kwargs)


class RectilinearMesh(StructuredMesh):
    """A 3D rectilinear Cartesian mesh

    Parameters
    ----------
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    dimension : Iterable of int
        The number of mesh cells in each direction (x, y, z).
    n_dimension : int
        Number of mesh dimensions (always 3 for a RectilinearMesh).
    x_grid : numpy.ndarray
        1-D array of mesh boundary points along the x-axis.
    y_grid : numpy.ndarray
        1-D array of mesh boundary points along the y-axis.
    z_grid : numpy.ndarray
        1-D array of mesh boundary points along the z-axis.
    indices : Iterable of tuple
        An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
        (2, 1, 1), ...]
    bounding_box : openmc.BoundingBox
        Axis-aligned bounding box of the mesh as defined by the upper-right and
        lower-left coordinates.

    """

    def __init__(self, mesh_id: int = None, name: str = ''):
        super().__init__(mesh_id, name)

        self._x_grid = None
        self._y_grid = None
        self._z_grid = None

    @property
    def dimension(self):
        return (len(self.x_grid) - 1,
                len(self.y_grid) - 1,
                len(self.z_grid) - 1)

    @property
    def n_dimension(self):
        return 3

    @property
    def x_grid(self):
        return self._x_grid

    @x_grid.setter
    def x_grid(self, grid):
        cv.check_type('mesh x_grid', grid, Iterable, Real)
        self._x_grid = np.asarray(grid, dtype=float)

    @property
    def y_grid(self):
        return self._y_grid

    @y_grid.setter
    def y_grid(self, grid):
        cv.check_type('mesh y_grid', grid, Iterable, Real)
        self._y_grid = np.asarray(grid, dtype=float)

    @property
    def z_grid(self):
        return self._z_grid

    @z_grid.setter
    def z_grid(self, grid):
        cv.check_type('mesh z_grid', grid, Iterable, Real)
        self._z_grid = np.asarray(grid, dtype=float)

    @property
    def _grids(self):
        return (self.x_grid, self.y_grid, self.z_grid)

    @property
    def lower_left(self):
        return np.array([self.x_grid[0], self.y_grid[0], self.z_grid[0]])

    @property
    def upper_right(self):
        return np.array([self.x_grid[-1], self.y_grid[-1], self.z_grid[-1]])

    @property
    def volumes(self):
        """Return Volumes for every mesh cell

        Returns
        -------
        volumes : numpy.ndarray
            Volumes

        """
        self._volume_dim_check()
        V_x = np.diff(self.x_grid)
        V_y = np.diff(self.y_grid)
        V_z = np.diff(self.z_grid)

        return np.multiply.outer(np.outer(V_x, V_y), V_z)

    @property
    def total_volume(self):
        return np.sum(self.volumes)

    @property
    def indices(self):
        nx = len(self.x_grid) - 1
        ny = len(self.y_grid) - 1
        nz = len(self.z_grid) - 1
        return ((x, y, z)
                for z in range(1, nz + 1)
                for y in range(1, ny + 1)
                for x in range(1, nx + 1))

    def __repr__(self):
        fmt = '{0: <16}{1}{2}\n'
        string = super().__repr__()
        string += fmt.format('\tDimensions', '=\t', self.n_dimension)
        x_grid_str = str(self._x_grid) if self._x_grid is None else len(self._x_grid)
        string += fmt.format('\tN X pnts:', '=\t', x_grid_str)
        if self._x_grid is not None:
            string += fmt.format('\tX Min:', '=\t', self._x_grid[0])
            string += fmt.format('\tX Max:', '=\t', self._x_grid[-1])
        y_grid_str = str(self._y_grid) if self._y_grid is None else len(self._y_grid)
        string += fmt.format('\tN Y pnts:', '=\t', y_grid_str)
        if self._y_grid is not None:
            string += fmt.format('\tY Min:', '=\t', self._y_grid[0])
            string += fmt.format('\tY Max:', '=\t', self._y_grid[-1])
        z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
        string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
        if self._z_grid is not None:
            string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
            string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
        return string

    @classmethod
    def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
        # Read and assign mesh properties
        mesh = cls(mesh_id=mesh_id, name=name)
        mesh.x_grid = group['x_grid'][()]
        mesh.y_grid = group['y_grid'][()]
        mesh.z_grid = group['z_grid'][()]

        return mesh

    @classmethod
    def from_xml_element(cls, elem: ET.Element):
        """Generate a rectilinear mesh from an XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element

        Returns
        -------
        openmc.RectilinearMesh
            Rectilinear mesh object

        """
        mesh_id = int(get_text(elem, 'id'))
        mesh = cls(mesh_id=mesh_id)
        mesh.x_grid = get_elem_list(elem, "x_grid", float)
        mesh.y_grid = get_elem_list(elem, "y_grid", float)
        mesh.z_grid = get_elem_list(elem, "z_grid", float)

        return mesh

    def to_xml_element(self):
        """Return XML representation of the mesh

        Returns
        -------
        element : lxml.etree._Element
            XML element containing mesh data

        """

        element = super().to_xml_element()
        element.set("type", "rectilinear")

        subelement = ET.SubElement(element, "x_grid")
        subelement.text = ' '.join(map(str, self.x_grid))

        subelement = ET.SubElement(element, "y_grid")
        subelement.text = ' '.join(map(str, self.y_grid))

        subelement = ET.SubElement(element, "z_grid")
        subelement.text = ' '.join(map(str, self.z_grid))

        return element


class CylindricalMesh(StructuredMesh):
    """A 3D cylindrical mesh

    Parameters
    ----------
    r_grid : numpy.ndarray
        1-D array of mesh boundary points along the r-axis
        Requirement is r >= 0.
    z_grid : numpy.ndarray
        1-D array of mesh boundary points along the z-axis relative to the
        origin.
    phi_grid : numpy.ndarray
        1-D array of mesh boundary points along the phi-axis in radians.
        The default value is [0, 2π], i.e. the full phi range.
    origin : numpy.ndarray
        1-D array of length 3 the (x,y,z) origin of the mesh in
        cartesian coordinates
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    dimension : Iterable of int
        The number of mesh cells in each direction (r_grid, phi_grid, z_grid).
    n_dimension : int
        Number of mesh dimensions (always 3 for a CylindricalMesh).
    r_grid : numpy.ndarray
        1-D array of mesh boundary points along the r-axis.
        Requirement is r >= 0.
    phi_grid : numpy.ndarray
        1-D array of mesh boundary points along the phi-axis in radians.
        The default value is [0, 2π], i.e. the full phi range.
    z_grid : numpy.ndarray
        1-D array of mesh boundary points along the z-axis relative to the
        origin.
    origin : numpy.ndarray
        1-D array of length 3 the (x,y,z) origin of the mesh in
        cartesian coordinates
    indices : Iterable of tuple
        An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
        (2, 1, 1), ...]
    lower_left : Iterable of float
        The lower-left corner of the structured mesh. If only two coordinate
        are given, it is assumed that the mesh is an x-y mesh.
    upper_right : Iterable of float
        The upper-right corner of the structured mesh. If only two coordinate
        are given, it is assumed that the mesh is an x-y mesh.
    bounding_box : openmc.BoundingBox
        Axis-aligned bounding box of the mesh as defined by the upper-right and
        lower-left coordinates.

    """

    def __init__(
        self,
        r_grid: Sequence[float],
        z_grid: Sequence[float],
        phi_grid: Sequence[float] = (0, 2*pi),
        origin: Sequence[float] = (0., 0., 0.),
        mesh_id: int | None = None,
        name: str = '',
    ):
        super().__init__(mesh_id, name)

        self.r_grid = r_grid
        self.phi_grid = phi_grid
        self.z_grid = z_grid
        self.origin = origin

    @property
    def dimension(self):
        return (len(self.r_grid) - 1,
                len(self.phi_grid) - 1,
                len(self.z_grid) - 1)

    @property
    def n_dimension(self):
        return 3

    @property
    def origin(self):
        return self._origin

    @origin.setter
    def origin(self, coords):
        cv.check_type('mesh origin', coords, Iterable, Real)
        cv.check_length("mesh origin", coords, 3)
        self._origin = np.asarray(coords)

    @property
    def r_grid(self):
        return self._r_grid

    @r_grid.setter
    def r_grid(self, grid):
        cv.check_type('mesh r_grid', grid, Iterable, Real)
        cv.check_length('mesh r_grid', grid, 2)
        cv.check_increasing('mesh r_grid', grid)
        self._r_grid = np.asarray(grid, dtype=float)

    @property
    def phi_grid(self):
        return self._phi_grid

    @phi_grid.setter
    def phi_grid(self, grid):
        cv.check_type('mesh phi_grid', grid, Iterable, Real)
        cv.check_length('mesh phi_grid', grid, 2)
        cv.check_increasing('mesh phi_grid', grid)
        grid = np.asarray(grid, dtype=float)
        if np.any((grid < 0.0) | (grid > 2*pi)):
            raise ValueError("phi_grid values must be in [0, 2π].")
        self._phi_grid = grid

    @property
    def z_grid(self):
        return self._z_grid

    @z_grid.setter
    def z_grid(self, grid):
        cv.check_type('mesh z_grid', grid, Iterable, Real)
        cv.check_length('mesh z_grid', grid, 2)
        cv.check_increasing('mesh z_grid', grid)
        self._z_grid = np.asarray(grid, dtype=float)

    @property
    def _grids(self):
        return (self.r_grid, self.phi_grid, self.z_grid)

    @property
    def indices(self):
        nr, np, nz = self.dimension
        np = len(self.phi_grid) - 1
        nz = len(self.z_grid) - 1
        return ((r, p, z)
                for z in range(1, nz + 1)
                for p in range(1, np + 1)
                for r in range(1, nr + 1))

    @property
    def lower_left(self):
        return np.array((
            self.origin[0] - self.r_grid[-1],
            self.origin[1] - self.r_grid[-1],
            self.origin[2] + self.z_grid[0]
        ))

    @property
    def upper_right(self):
        return np.array((
            self.origin[0] + self.r_grid[-1],
            self.origin[1] + self.r_grid[-1],
            self.origin[2] + self.z_grid[-1]
        ))

    def __repr__(self):
        fmt = '{0: <16}{1}{2}\n'
        string = super().__repr__()
        string += fmt.format('\tDimensions', '=\t', self.n_dimension)
        string += fmt.format('\tOrigin', '=\t', self.origin)
        r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
        string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
        if self._r_grid is not None:
            string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
            string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
        phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
        string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
        if self._phi_grid is not None:
            string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
            string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
        z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
        string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
        if self._z_grid is not None:
            string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
            string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
        return string

    def get_indices_at_coords(
            self,
            coords: Sequence[float]
        ) -> tuple[int, int, int]:
        """Finds the index of the mesh voxel at the specified x,y,z coordinates.

        .. versionadded:: 0.15.0

        Parameters
        ----------
        coords : Sequence[float]
            The x, y, z axis coordinates

        Returns
        -------
        tuple[int, int, int]
            The r, phi, z indices

        """
        r_value_from_origin = sqrt((coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)

        if r_value_from_origin < self.r_grid[0] or r_value_from_origin > self.r_grid[-1]:
            raise ValueError(
                f'The specified x, y ({coords[0]}, {coords[1]}) combine to give an r value of '
                f'{r_value_from_origin} from the origin of {self.origin}.which '
                f'is outside the origin absolute r grid values {self.r_grid}.'
            )

        r_index = np.searchsorted(self.r_grid, r_value_from_origin) - 1

        z_grid_values = np.array(self.z_grid) + self.origin[2]

        if coords[2] < z_grid_values[0] or coords[2] > z_grid_values[-1]:
            raise ValueError(
                f'The specified z value ({coords[2]}) from the z origin of '
                f'{self.origin[-1]} is outside of the absolute z grid range {z_grid_values}.'
            )

        z_index = np.argmax(z_grid_values > coords[2]) - 1

        delta_x = coords[0] - self.origin[0]
        delta_y = coords[1] - self.origin[1]
        # atan2 returns values in -pi to +pi range
        phi_value = atan2(delta_y, delta_x)
        if delta_x < 0 and delta_y < 0:
            # returned phi_value anticlockwise and negative
            phi_value += 2 * pi
        if delta_x > 0 and delta_y < 0:
            # returned phi_value anticlockwise and negative
            phi_value += 2 * pi

        phi_grid_values = np.array(self.phi_grid)

        if phi_value < phi_grid_values[0] or phi_value > phi_grid_values[-1]:
            raise ValueError(
                f'The phi value ({phi_value}) resulting from the specified x, y '
                f'values is outside of the absolute  phi grid range {phi_grid_values}.'
            )
        phi_index = np.argmax(phi_grid_values > phi_value) - 1

        return (r_index, phi_index, z_index)

    @classmethod
    def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
        # Read and assign mesh properties
        mesh = cls(
            mesh_id=mesh_id,
            name=name,
            r_grid = group['r_grid'][()],
            phi_grid = group['phi_grid'][()],
            z_grid = group['z_grid'][()],
        )
        if 'origin' in group:
            mesh.origin = group['origin'][()]

        return mesh

    @classmethod
    def from_domain(
        cls,
        domain: HasBoundingBox | BoundingBox,
        dimension: Sequence[int] = (10, 10, 10),
        mesh_id: int | None = None,
        phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
        name: str = '',
        enclose_domain: bool = False
    ):
        """Create CylindricalMesh from a domain using its bounding box.

        Parameters
        ----------
        domain : HasBoundingBox | openmc.BoundingBox
            The object passed in will be used as a template for this mesh. The
            bounding box of the property of the object passed will be used to
            set the r_grid, z_grid ranges. Alternatively, a
            :class:`openmc.BoundingBox` can be passed directly.
        dimension : Iterable of int
            The number of equally spaced mesh cells in each direction (r_grid,
            phi_grid, z_grid)
        mesh_id : int
            Unique identifier for the mesh
        phi_grid_bounds : numpy.ndarray
            Mesh bounds points along the phi-axis in radians. The default value
            is (0, 2π), i.e., the full phi range.
        name : str
            Name of the mesh
        enclose_domain : bool
            If True, the mesh will encompass the bounding box of the domain. If
            False, the mesh will be inscribed within the domain's bounding box.

        Returns
        -------
        openmc.CylindricalMesh
            CylindricalMesh instance

        """
        if isinstance(domain, BoundingBox):
            cached_bb = domain
        elif hasattr(domain, 'bounding_box'):
            cached_bb = domain.bounding_box
        else:
            raise TypeError("Domain must be a BoundingBox or have a "
                            "bounding_box property")

        if enclose_domain:
            outer_radius = 0.5 * np.linalg.norm(cached_bb.width[:2])
        else:
            outer_radius = 0.5 * min(cached_bb.width[:2])

        r_grid = np.linspace(
            0,
            outer_radius,
            num=dimension[0]+1
        )
        phi_grid = np.linspace(
            phi_grid_bounds[0],
            phi_grid_bounds[1],
            num=dimension[1]+1
        )
        z_grid = np.linspace(
            cached_bb[0][2],
            cached_bb[1][2],
            num=dimension[2]+1
        )
        origin = (cached_bb.center[0], cached_bb.center[1], z_grid[0])

        # make z-grid relative to the origin
        z_grid -= origin[2]

        mesh = cls(
            r_grid=r_grid,
            z_grid=z_grid,
            phi_grid=phi_grid,
            mesh_id=mesh_id,
            name=name,
            origin=origin
        )

        return mesh

    def to_xml_element(self):
        """Return XML representation of the mesh

        Returns
        -------
        element : lxml.etree._Element
            XML element containing mesh data

        """

        element = super().to_xml_element()
        element.set("type", "cylindrical")

        subelement = ET.SubElement(element, "r_grid")
        subelement.text = ' '.join(map(str, self.r_grid))

        subelement = ET.SubElement(element, "phi_grid")
        subelement.text = ' '.join(map(str, self.phi_grid))

        subelement = ET.SubElement(element, "z_grid")
        subelement.text = ' '.join(map(str, self.z_grid))

        subelement = ET.SubElement(element, "origin")
        subelement.text = ' '.join(map(str, self.origin))

        return element

    @classmethod
    def from_xml_element(cls, elem: ET.Element):
        """Generate a cylindrical mesh from an XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element

        Returns
        -------
        openmc.CylindricalMesh
            Cylindrical mesh object

        """

        mesh_id = int(get_text(elem, 'id'))
        mesh = cls(
            r_grid = get_elem_list(elem, "r_grid", float),
            phi_grid = get_elem_list(elem, "phi_grid", float),
            z_grid = get_elem_list(elem, "z_grid", float),
            origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
            mesh_id=mesh_id,
        )

        return mesh

    @property
    def volumes(self):
        """Return Volumes for every mesh cell

        Returns
        -------
        volumes : Iterable of float
            Volumes

        """
        self._volume_dim_check()
        V_r = np.diff(np.asarray(self.r_grid)**2 / 2)
        V_p = np.diff(self.phi_grid)
        V_z = np.diff(self.z_grid)

        return np.multiply.outer(np.outer(V_r, V_p), V_z)

    @property
    def vertices(self):
        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
        return self._convert_to_cartesian(self.vertices_cylindrical, self.origin)

    @property
    def vertices_cylindrical(self):
        """Returns vertices of the mesh in cylindrical coordinates.
        """
        return super().vertices

    @property
    def centroids(self):
        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
        return self._convert_to_cartesian(self.centroids_cylindrical, self.origin)

    @property
    def centroids_cylindrical(self):
        """Returns centroids of the mesh in cylindrical coordinates.
        """
        return super().centroids

    @staticmethod
    def _convert_to_cartesian(arr, origin: Sequence[float]):
        """Converts an array with r, phi, z values in the last dimension (shape (..., 3))
        to Cartesian coordinates.
        """
        x = arr[..., 0] * np.cos(arr[..., 1]) + origin[0]
        y = arr[..., 0] * np.sin(arr[..., 1]) + origin[1]
        arr[..., 0] = x
        arr[..., 1] = y
        arr[..., 2] += origin[2]
        return arr


class SphericalMesh(StructuredMesh):
    """A 3D spherical mesh

    Parameters
    ----------
    r_grid : numpy.ndarray
        1-D array of mesh boundary points along the r-axis.
        Requirement is r >= 0.
    phi_grid : numpy.ndarray
        1-D array of mesh boundary points along the phi-axis in radians.
        The default value is [0, 2π], i.e. the full phi range.
    theta_grid : numpy.ndarray
        1-D array of mesh boundary points along the theta-axis in radians.
        The default value is [0, π], i.e. the full theta range.
    origin : numpy.ndarray
        1-D array of length 3 the (x,y,z) origin of the mesh in
        cartesian coordinates
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    dimension : Iterable of int
        The number of mesh cells in each direction (r_grid,
        theta_grid, phi_grid).
    n_dimension : int
        Number of mesh dimensions (always 3 for a SphericalMesh).
    r_grid : numpy.ndarray
        1-D array of mesh boundary points along the r-axis.
        Requirement is r >= 0.
    theta_grid : numpy.ndarray
        1-D array of mesh boundary points along the theta-axis in radians.
        The default value is [0, π], i.e. the full theta range.
    phi_grid : numpy.ndarray
        1-D array of mesh boundary points along the phi-axis in radians.
        The default value is [0, 2π], i.e. the full phi range.
    origin : numpy.ndarray
        1-D array of length 3 the (x,y,z) origin of the mesh in
        cartesian coordinates
    indices : Iterable of tuple
        An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
        (2, 1, 1), ...]
    lower_left : numpy.ndarray
        The lower-left corner of the structured mesh. If only two coordinate
        are given, it is assumed that the mesh is an x-y mesh.
    upper_right : numpy.ndarray
        The upper-right corner of the structured mesh. If only two coordinate
        are given, it is assumed that the mesh is an x-y mesh.
    bounding_box : openmc.BoundingBox
        Axis-aligned bounding box of the mesh as defined by the upper-right and
        lower-left coordinates.

    """

    def __init__(
        self,
        r_grid: Sequence[float],
        phi_grid: Sequence[float] = (0, 2*pi),
        theta_grid: Sequence[float] = (0, pi),
        origin: Sequence[float] = (0., 0., 0.),
        mesh_id: int | None = None,
        name: str = '',
    ):
        super().__init__(mesh_id, name)

        self.r_grid = r_grid
        self.theta_grid = theta_grid
        self.phi_grid = phi_grid
        self.origin = origin

    @property
    def dimension(self):
        return (len(self.r_grid) - 1,
                len(self.theta_grid) - 1,
                len(self.phi_grid) - 1)

    @property
    def n_dimension(self):
        return 3

    @property
    def origin(self):
        return self._origin

    @origin.setter
    def origin(self, coords):
        cv.check_type('mesh origin', coords, Iterable, Real)
        cv.check_length("mesh origin", coords, 3)
        self._origin = np.asarray(coords, dtype=float)

    @property
    def r_grid(self):
        return self._r_grid

    @r_grid.setter
    def r_grid(self, grid):
        cv.check_type('mesh r_grid', grid, Iterable, Real)
        cv.check_length('mesh r_grid', grid, 2)
        cv.check_increasing('mesh r_grid', grid)
        self._r_grid = np.asarray(grid, dtype=float)

    @property
    def theta_grid(self):
        return self._theta_grid

    @theta_grid.setter
    def theta_grid(self, grid):
        cv.check_type('mesh theta_grid', grid, Iterable, Real)
        cv.check_length('mesh theta_grid', grid, 2)
        cv.check_increasing('mesh theta_grid', grid)
        grid = np.asarray(grid, dtype=float)
        if np.any((grid < 0.0) | (grid > pi)):
            raise ValueError("theta_grid values must be in [0, π].")
        self._theta_grid = grid

    @property
    def phi_grid(self):
        return self._phi_grid

    @phi_grid.setter
    def phi_grid(self, grid):
        cv.check_type('mesh phi_grid', grid, Iterable, Real)
        cv.check_length('mesh phi_grid', grid, 2)
        cv.check_increasing('mesh phi_grid', grid)
        grid = np.asarray(grid, dtype=float)
        if np.any((grid < 0.0) | (grid > 2*pi)):
            raise ValueError("phi_grid values must be in [0, 2π].")
        self._phi_grid = grid

    @property
    def _grids(self):
        return (self.r_grid, self.theta_grid, self.phi_grid)

    @property
    def indices(self):
        nr, nt, np = self.dimension
        nt = len(self.theta_grid) - 1
        np = len(self.phi_grid) - 1
        return ((r, t, p)
                for p in range(1, np + 1)
                for t in range(1, nt + 1)
                for r in range(1, nr + 1))

    @property
    def lower_left(self):
        r = self.r_grid[-1]
        return np.array((self.origin[0] - r, self.origin[1] - r, self.origin[2] - r))

    @property
    def upper_right(self):
        r = self.r_grid[-1]
        return np.array((self.origin[0] + r, self.origin[1] + r, self.origin[2] + r))

    def __repr__(self):
        fmt = '{0: <16}{1}{2}\n'
        string = super().__repr__()
        string += fmt.format('\tDimensions', '=\t', self.n_dimension)
        string += fmt.format('\tOrigin', '=\t', self.origin)
        r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
        string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
        if self._r_grid is not None:
            string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
            string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
        theta_grid_str = str(self._theta_grid) if self._theta_grid is None else len(self._theta_grid)
        string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str)
        if self._theta_grid is not None:
            string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0])
            string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1])
        phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
        string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
        if self._phi_grid is not None:
            string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
            string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
        return string

    @classmethod
    def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
        # Read and assign mesh properties
        mesh = cls(
            r_grid = group['r_grid'][()],
            theta_grid = group['theta_grid'][()],
            phi_grid = group['phi_grid'][()],
            mesh_id=mesh_id,
            name=name
        )
        if 'origin' in group:
            mesh.origin = group['origin'][()]

        return mesh

    @classmethod
    def from_domain(
        cls,
        domain: HasBoundingBox | BoundingBox,
        dimension: Sequence[int] = (10, 10, 10),
        mesh_id: int | None = None,
        phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
        theta_grid_bounds: Sequence[float] = (0.0, pi),
        name: str = '',
        enclose_domain: bool = False
    ):
        """Create SphericalMesh from a domain using its bounding box.

        Parameters
        ----------
        domain : HasBoundingBox | openmc.BoundingBox
            The object passed in will be used as a template for this mesh. The
            bounding box of the property of the object passed will be used to
            set the r_grid, phi_grid, and theta_grid ranges. Alternatively, a
            :class:`openmc.BoundingBox` can be passed directly.
        dimension : Iterable of int
            The number of equally spaced mesh cells in each direction (r_grid,
            phi_grid, theta_grid). Spacing is in angular space (radians) for
            phi and theta, and in absolute space for r.
        mesh_id : int
            Unique identifier for the mesh
        phi_grid_bounds : numpy.ndarray
            Mesh bounds points along the phi-axis in radians. The default value
            is (0, 2π), i.e., the full phi range.
        theta_grid_bounds : numpy.ndarray
            Mesh bounds points along the theta-axis in radians. The default value
            is (0, π), i.e., the full theta range.
        name : str
            Name of the mesh
        enclose_domain : bool
            If True, the mesh will encompass the bounding box of the domain. If
            False, the mesh will be inscribed within the domain's bounding box.

        Returns
        -------
        openmc.SphericalMesh
            SphericalMesh instance

        """
        if isinstance(domain, BoundingBox):
            cached_bb = domain
        elif hasattr(domain, 'bounding_box'):
            cached_bb = domain.bounding_box
        else:
            raise TypeError("Domain must be a BoundingBox or have a "
                            "bounding_box property")

        if enclose_domain:
            outer_radius = 0.5 * np.linalg.norm(cached_bb.width)
        else:
            outer_radius = 0.5 * min(cached_bb.width)

        r_grid = np.linspace(0, outer_radius, num=dimension[0] + 1)
        theta_grid = np.linspace(
            theta_grid_bounds[0],
            theta_grid_bounds[1],
            num=dimension[1]+1
        )
        phi_grid = np.linspace(
            phi_grid_bounds[0],
            phi_grid_bounds[1],
            num=dimension[2]+1
        )
        origin = np.array([
            cached_bb.center[0], cached_bb.center[1], cached_bb.center[2]])

        return cls(r_grid=r_grid, phi_grid=phi_grid, theta_grid=theta_grid,
                   origin=origin, mesh_id=mesh_id, name=name)

    def to_xml_element(self):
        """Return XML representation of the mesh

        Returns
        -------
        element : lxml.etree._Element
            XML element containing mesh data

        """

        element = super().to_xml_element()
        element.set("type", "spherical")

        subelement = ET.SubElement(element, "r_grid")
        subelement.text = ' '.join(map(str, self.r_grid))

        subelement = ET.SubElement(element, "theta_grid")
        subelement.text = ' '.join(map(str, self.theta_grid))

        subelement = ET.SubElement(element, "phi_grid")
        subelement.text = ' '.join(map(str, self.phi_grid))

        subelement = ET.SubElement(element, "origin")
        subelement.text = ' '.join(map(str, self.origin))

        return element

    @classmethod
    def from_xml_element(cls, elem: ET.Element):
        """Generate a spherical mesh from an XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element

        Returns
        -------
        openmc.SphericalMesh
            Spherical mesh object

        """
        mesh_id = int(get_text(elem, 'id'))
        mesh = cls(
            mesh_id=mesh_id,
            r_grid = get_elem_list(elem, "r_grid", float),
            theta_grid = get_elem_list(elem, "theta_grid", float),
            phi_grid = get_elem_list(elem, "phi_grid", float),
            origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
        )

        return mesh

    @property
    def volumes(self):
        """Return Volumes for every mesh cell

        Returns
        -------
        volumes : Iterable of float
            Volumes

        """
        self._volume_dim_check()
        V_r = np.diff(np.asarray(self.r_grid)**3 / 3)
        V_t = np.diff(-np.cos(self.theta_grid))
        V_p = np.diff(self.phi_grid)

        return np.multiply.outer(np.outer(V_r, V_t), V_p)

    @property
    def vertices(self):
        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
        return self._convert_to_cartesian(self.vertices_spherical, self.origin)

    @property
    def vertices_spherical(self):
        """Returns vertices of the mesh in cylindrical coordinates.
        """
        return super().vertices

    @property
    def centroids(self):
        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
        return self._convert_to_cartesian(self.centroids_spherical, self.origin)

    @property
    def centroids_spherical(self):
        """Returns centroids of the mesh in cylindrical coordinates.
        """
        return super().centroids


    @staticmethod
    def _convert_to_cartesian(arr, origin: Sequence[float]):
        """Converts an array with r, theta, phi values in the last dimension (shape (..., 3))
        to Cartesian coordinates.
        """
        r_xy = arr[..., 0] * np.sin(arr[..., 1])
        x = r_xy * np.cos(arr[..., 2])
        y = r_xy * np.sin(arr[..., 2])
        z = arr[..., 0] * np.cos(arr[..., 1])
        arr[..., 0] = x + origin[0]
        arr[..., 1] = y + origin[1]
        arr[..., 2] = z + origin[2]
        return arr


def require_statepoint_data(func):
    @wraps(func)
    def wrapper(self: UnstructuredMesh, *args, **kwargs):
        if not self._has_statepoint_data:
            raise AttributeError(f'The "{func.__name__}" property requires '
                                 'information about this mesh to be loaded '
                                 'from a statepoint file.')
        return func(self, *args, **kwargs)
    return wrapper


class UnstructuredMesh(MeshBase):
    """A 3D unstructured mesh

    .. versionadded:: 0.12

    .. versionchanged:: 0.12.2
        Support for libMesh unstructured meshes was added.

    Parameters
    ----------
    filename : path-like
        Location of the unstructured mesh file. Supported files for 'moab'
        library are .h5 and .vtk. Supported files for 'libmesh' library are
        exodus mesh files .exo.
    library : {'moab', 'libmesh'}
        Mesh library used for the unstructured mesh tally
    mesh_id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    length_multiplier: float
        Constant multiplier to apply to mesh coordinates
    options : str, optional
        Special options that control spatial search data structures used. This
        is currently only used to set `parameters
        <https://tinyurl.com/kdtree-params>`_ for MOAB's AdaptiveKDTree. If
        None, OpenMC internally uses a default of "MAX_DEPTH=20;PLANE_SET=2;".

    Attributes
    ----------
    id : int
        Unique identifier for the mesh
    name : str
        Name of the mesh
    filename : str
        Name of the file containing the unstructured mesh
    length_multiplier: float
        Multiplicative factor to apply to mesh coordinates
    library : {'moab', 'libmesh'}
        Mesh library used for the unstructured mesh tally
    options : str
        Special options that control spatial search data structures used. This
        is currently only used to set `parameters
        <https://tinyurl.com/kdtree-params>`_ for MOAB's AdaptiveKDTree. If
        None, OpenMC internally uses a default of "MAX_DEPTH=20;PLANE_SET=2;".
    output : bool
        Indicates whether or not automatic tally output should be generated for
        this mesh
    volumes : Iterable of float
        Volumes of the unstructured mesh elements
    centroids : numpy.ndarray
        Centroids of the mesh elements with array shape (n_elements, 3)

    vertices : numpy.ndarray
        Coordinates of the mesh vertices with array shape (n_elements, 3)

        .. versionadded:: 0.13.1
    connectivity : numpy.ndarray
        Connectivity of the elements with array shape (n_elements, 8)

        .. versionadded:: 0.13.1
    element_types : Iterable of integers
        Mesh element types

        .. versionadded:: 0.13.1
    total_volume : float
        Volume of the unstructured mesh in total
    bounding_box : openmc.BoundingBox
        Axis-aligned bounding box of the mesh as defined by the upper-right and
        lower-left coordinates.

    """

    _UNSUPPORTED_ELEM = -1
    _LINEAR_TET = 0
    _LINEAR_HEX = 1
    _VTK_TETRA = 10

    def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None,
                 name: str = '', length_multiplier: float = 1.0,
                 options: str | None = None):
        super().__init__(mesh_id, name)
        self.filename = filename
        self._volumes = None
        self._n_elements = None
        self._conectivity = None
        self._vertices = None
        self.library = library
        self._output = False
        self.length_multiplier = length_multiplier
        self.options = options
        self._has_statepoint_data = False

    @property
    def filename(self):
        return self._filename

    @filename.setter
    def filename(self, filename):
        cv.check_type('Unstructured Mesh filename', filename, PathLike)
        self._filename = input_path(filename)

    @property
    def library(self):
        return self._library

    @library.setter
    def library(self, lib: str):
        cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh'))
        self._library = lib

    @property
    def options(self) -> str | None:
        return self._options

    @options.setter
    def options(self, options: str | None):
        cv.check_type('options', options, (str, type(None)))
        self._options = options

    @property
    @require_statepoint_data
    def size(self):
        return self._size

    @size.setter
    def size(self, size: int):
        cv.check_type("Unstructured mesh size", size, Integral)
        self._size = size

    @property
    def output(self):
        return self._output

    @output.setter
    def output(self, val: bool):
        cv.check_type("Unstructured mesh output value", val, bool)
        self._output = val

    @property
    @require_statepoint_data
    def volumes(self):
        """Return Volumes for every mesh cell if
        populated by a StatePoint file

        Returns
        -------
        volumes : numpy.ndarray
            Volumes

        """
        return self._volumes

    @volumes.setter
    def volumes(self, volumes: Iterable[Real]):
        cv.check_type("Unstructured mesh volumes", volumes, Iterable, Real)
        self._volumes = volumes

    @property
    @require_statepoint_data
    def total_volume(self):
        return np.sum(self.volumes)

    @property
    @require_statepoint_data
    def vertices(self):
        return self._vertices

    @property
    @require_statepoint_data
    def connectivity(self):
        return self._connectivity

    @property
    @require_statepoint_data
    def element_types(self):
        return self._element_types

    @property
    @require_statepoint_data
    def centroids(self):
        return np.array([self.centroid(i) for i in range(self.n_elements)])

    @property
    @require_statepoint_data
    def n_elements(self):
        if self._n_elements is None:
            raise RuntimeError("No information about this mesh has "
                               "been loaded from a statepoint file.")
        return self._n_elements

    @n_elements.setter
    def n_elements(self, val: int):
        cv.check_type('Number of elements', val, Integral)
        self._n_elements = val

    @property
    def length_multiplier(self):
        return self._length_multiplier

    @length_multiplier.setter
    def length_multiplier(self, length_multiplier):
        cv.check_type("Unstructured mesh length multiplier",
                      length_multiplier,
                      Real)
        self._length_multiplier = length_multiplier

    @property
    def dimension(self):
        return (self.n_elements,)

    @property
    def n_dimension(self):
        return 3

    @property
    @require_statepoint_data
    def indices(self):
        return [(i,) for i in range(self.n_elements)]

    @property
    def has_statepoint_data(self) -> bool:
        return self._has_statepoint_data

    def __repr__(self):
        string = super().__repr__()
        string += '{: <16}=\t{}\n'.format('\tFilename', self.filename)
        string += '{: <16}=\t{}\n'.format('\tMesh Library', self.library)
        if self.length_multiplier != 1.0:
            string += '{: <16}=\t{}\n'.format('\tLength multiplier',
                                              self.length_multiplier)
        if self.options is not None:
            string += '{: <16}=\t{}\n'.format('\tOptions', self.options)
        return string

    @property
    @require_statepoint_data
    def lower_left(self):
        return self.vertices.min(axis=0)

    @property
    @require_statepoint_data
    def upper_right(self):
        return self.vertices.max(axis=0)

    @require_statepoint_data
    def centroid(self, bin: int):
        """Return the vertex averaged centroid of an element

        Parameters
        ----------
        bin : int
            Bin ID for the returned centroid

        Returns
        -------
        numpy.ndarray
            x, y, z values of the element centroid

        """
        conn = self.connectivity[bin]
        # remove invalid connectivity values
        conn = conn[conn >= 0]
        coords = self.vertices[conn]
        return coords.mean(axis=0)

    def write_vtk_mesh(self, **kwargs):
        """Map data to unstructured VTK mesh elements.

        .. deprecated:: 0.13
          Use :func:`UnstructuredMesh.write_data_to_vtk` instead.

        Parameters
        ----------
        filename : str or pathlib.Path
            Name of the VTK file to write.
        datasets : dict
            Dictionary whose keys are the data labels
            and values are the data sets.
        volume_normalization : bool
            Whether or not to normalize the data by the
            volume of the mesh elements
        """
        warnings.warn(
            "The 'UnstructuredMesh.write_vtk_mesh' method has been renamed "
            "to 'write_data_to_vtk' and will be removed in a future version "
            " of OpenMC.",
            FutureWarning,
        )
        self.write_data_to_vtk(**kwargs)

    def write_data_to_vtk(
        self,
        filename: PathLike | None = None,
        datasets: dict | None = None,
        volume_normalization: bool = True,
    ):
        """Map data to unstructured VTK mesh elements.

        If filename is None, then a filename will be generated based on the mesh
        ID, and exported to VTK format.

        Parameters
        ----------
        filename : str or pathlib.Path
            Name of the VTK file to write. If the filename ends in '.vtkhdf'
            then a VTKHDF format file will be written. If the filename ends in
            '.vtu' then a binary VTU format file will be written. If the
            filename ends in '.vtk' then a legacy VTK file will be written.
        datasets : dict
            Dictionary whose keys are the data labels and values are numpy
            appropriately sized arrays of the data
        volume_normalization : bool
            Whether or not to normalize the data by the volume of the mesh
            elements
        """

        if Path(filename).suffix == ".vtkhdf":

            self._write_data_to_vtk_hdf5_format(
                filename=filename,
                datasets=datasets,
                volume_normalization=volume_normalization,
            )

        elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu":

            self._write_data_to_vtk_ascii_format(
                filename=filename,
                datasets=datasets,
                volume_normalization=volume_normalization,
            )

        else:
            raise ValueError(
                "Unsupported file extension, The filename must end with "
                "'.vtkhdf', '.vtu' or '.vtk'"
            )

    def _write_data_to_vtk_ascii_format(
        self,
        filename: PathLike | None = None,
        datasets: dict | None = None,
        volume_normalization: bool = True,
    ):
        from vtkmodules.util import numpy_support
        from vtkmodules import vtkCommonCore
        from vtkmodules import vtkCommonDataModel
        from vtkmodules import vtkIOLegacy
        from vtkmodules import vtkIOXML

        if self.connectivity is None or self.vertices is None:
            raise RuntimeError("This mesh has not been loaded from a statepoint file.")

        if filename is None:
            filename = f"mesh_{self.id}.vtk"

        if Path(filename).suffix == ".vtk":
            writer = vtkIOLegacy.vtkUnstructuredGridWriter()

        elif Path(filename).suffix == ".vtu":
            writer = vtkIOXML.vtkXMLUnstructuredGridWriter()
            writer.SetCompressorTypeToZLib()
            writer.SetDataModeToBinary()

        writer.SetFileName(str(filename))

        grid = vtkCommonDataModel.vtkUnstructuredGrid()

        points = vtkCommonCore.vtkPoints()
        points.SetData(numpy_support.numpy_to_vtk(self.vertices))
        grid.SetPoints(points)

        n_skipped = 0
        for elem_type, conn in zip(self.element_types, self.connectivity):
            if elem_type == self._LINEAR_TET:
                elem = vtkCommonDataModel.vtkTetra()
            elif elem_type == self._LINEAR_HEX:
                elem = vtkCommonDataModel.vtkHexahedron()
            elif elem_type == self._UNSUPPORTED_ELEM:
                n_skipped += 1
                continue
            else:
                raise RuntimeError(f"Invalid element type {elem_type} found")

            for i, c in enumerate(conn):
                if c == -1:
                    break
                elem.GetPointIds().SetId(i, c)

            grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())

        if n_skipped > 0:
            warnings.warn(
                f"{n_skipped} elements were not written because "
                "they are not of type linear tet/hex"
            )

        # check that datasets are the correct size
        datasets_out = []
        if datasets is not None:
            for name, data in datasets.items():
                if data.shape != self.dimension:
                    raise ValueError(
                        f'Cannot apply dataset "{name}" with '
                        f"shape {data.shape} to mesh {self.id} "
                        f"with dimensions {self.dimension}"
                    )

            if volume_normalization:
                for name, data in datasets.items():
                    if np.issubdtype(data.dtype, np.integer):
                        warnings.warn(
                            f'Integer data set "{name}" will '
                            "not be volume-normalized."
                        )
                        continue
                    data /= self.volumes

            # add data to the mesh
            for name, data in datasets.items():
                datasets_out.append(data)
                arr = vtkCommonCore.vtkDoubleArray()
                arr.SetName(name)
                arr.SetNumberOfTuples(data.size)

                for i in range(data.size):
                    arr.SetTuple1(i, data.flat[i])
                grid.GetCellData().AddArray(arr)

        writer.SetInputData(grid)

        writer.Write()

    def _write_data_to_vtk_hdf5_format(
        self,
        filename: PathLike | None = None,
        datasets: dict | None = None,
        volume_normalization: bool = True,
    ):
        def append_dataset(dset, array):
            """Convenience function to append data to an HDF5 dataset"""
            origLen = dset.shape[0]
            dset.resize(origLen + array.shape[0], axis=0)
            dset[origLen:] = array

        if self.library != "moab":
            raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")

        # the self.connectivity contains arrays of length 8 to support hex
        # elements as well, in the case of tetrahedra mesh elements, the
        # last 4 values are -1 and are removed
        trimmed_connectivity = []
        for cell in self.connectivity:
            # Find the index of the first -1 value, if any
            first_negative_index = np.where(cell == -1)[0]
            if first_negative_index.size > 0:
                # Slice the array up to the first -1 value
                trimmed_connectivity.append(cell[: first_negative_index[0]])
            else:
                # No -1 values, append the whole cell
                trimmed_connectivity.append(cell)
        trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()

        # MOAB meshes supports tet elements only so we know it has 4 points per cell
        points_per_cell = 4

        # offsets are the indices of the first point of each cell in the array of points
        offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)

        for name, data in datasets.items():
            if data.shape != self.dimension:
                raise ValueError(
                    f'Cannot apply dataset "{name}" with '
                    f"shape {data.shape} to mesh {self.id} "
                    f"with dimensions {self.dimension}"
                )

        with h5py.File(filename, "w") as f:

            root = f.create_group("VTKHDF")
            vtk_file_format_version = (2, 1)
            root.attrs["Version"] = vtk_file_format_version
            ascii_type = "UnstructuredGrid".encode("ascii")
            root.attrs.create(
                "Type",
                ascii_type,
                dtype=h5py.string_dtype("ascii", len(ascii_type)),
            )

            # create hdf5 file structure
            root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
            root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
            root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
            root.create_dataset(
                "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
            )
            root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
            root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
            root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")

            append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
            append_dataset(root["Points"], self.vertices)
            append_dataset(
                root["NumberOfConnectivityIds"],
                np.array([len(trimmed_connectivity)]),
            )
            append_dataset(root["Connectivity"], trimmed_connectivity)
            append_dataset(root["NumberOfCells"], np.array([self.n_elements]))
            append_dataset(root["Offsets"], offsets)

            append_dataset(
                root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
            )

            cell_data_group = root.create_group("CellData")

            for name, data in datasets.items():

                cell_data_group.create_dataset(
                    name, (0,), maxshape=(None,), dtype="float64", chunks=True
                )

                if volume_normalization:
                    data /= self.volumes
                append_dataset(cell_data_group[name], data)

    @classmethod
    def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
        filename = group["filename"][()].decode()
        library = group["library"][()].decode()
        if "options" in group.attrs:
            options = group.attrs['options'].decode()
        else:
            options = None

        mesh = cls(
            filename=filename,
            library=library,
            mesh_id=mesh_id,
            name=name,
            options=options,
        )
        mesh._has_statepoint_data = True
        vol_data = group["volumes"][()]
        mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
        mesh.n_elements = mesh.volumes.size

        vertices = group["vertices"][()]
        mesh._vertices = vertices.reshape((-1, 3))
        connectivity = group["connectivity"][()]
        mesh._connectivity = connectivity.reshape((-1, 8))
        mesh._element_types = group["element_types"][()]

        if "length_multiplier" in group:
            mesh.length_multiplier = group["length_multiplier"][()]

        return mesh

    def to_xml_element(self):
        """Return XML representation of the mesh

        Returns
        -------
        element : lxml.etree._Element
            XML element containing mesh data

        """

        element = super().to_xml_element()
        element.set("type", "unstructured")

        element.set("library", self._library)
        if self.options is not None:
            element.set("options", self.options)
        subelement = ET.SubElement(element, "filename")
        subelement.text = str(self.filename)

        if self._length_multiplier != 1.0:
            element.set("length_multiplier", str(self.length_multiplier))

        return element

    @classmethod
    def from_xml_element(cls, elem: ET.Element):
        """Generate unstructured mesh object from XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element

        Returns
        -------
        openmc.UnstructuredMesh
            UnstructuredMesh generated from an XML element
        """
        mesh_id = int(get_text(elem, 'id'))
        filename = get_text(elem, 'filename')
        library = get_text(elem, 'library')
        length_multiplier = float(get_text(elem, 'length_multiplier', 1.0))
        options = get_text(elem, "options")

        return cls(filename, library, mesh_id, '', length_multiplier, options)


def _read_meshes(elem):
    """Generate dictionary of meshes from a given XML node

    Parameters
    ----------
    elem : lxml.etree._Element
        XML element

    Returns
    -------
    dict
        A dictionary with mesh IDs as keys and openmc.MeshBase
        instanaces as values
    """
    out = {}
    for mesh_elem in elem.findall('mesh'):
        mesh = MeshBase.from_xml_element(mesh_elem)
        out[mesh.id] = mesh

    return out


# hexahedron element connectivity
# lower-k connectivity offsets
_HEX_VERTEX_CONN = ((0, 0, 0),
                    (1, 0, 0),
                    (1, 1, 0),
                    (0, 1, 0))
# upper-k connectivity offsets
_HEX_VERTEX_CONN += ((0, 0, 1),
                     (1, 0, 1),
                     (1, 1, 1),
                     (0, 1, 1))

_N_HEX_VERTICES = 8

# lower-k connectivity offsets
_HEX_MIDPOINT_CONN = ((0, (0, 0, 0)),
                      (1, (1, 0, 0)),
                      (0, (0, 1, 0)),
                      (1, (0, 0, 0)))
# upper-k connectivity offsets
_HEX_MIDPOINT_CONN += ((0, (0, 0, 1)),
                       (1, (1, 0, 1)),
                       (0, (0, 1, 1)),
                       (1, (0, 0, 1)))
# mid-plane k connectivity
_HEX_MIDPOINT_CONN += ((2, (0, 0, 0)),
                       (2, (1, 0, 0)),
                       (2, (1, 1, 0)),
                       (2, (0, 1, 0)))