import collections
from copy import deepcopy
import meshio
import numpy as np
from ._filter import MeshFilter
from ._properties import (
_connections,
_face_areas,
_face_normals,
_faces,
_materials,
_qualities,
_volumes,
)
__all__ = [
"CellBlock",
"Mesh",
"from_meshio",
"from_pyvista",
]
CellBlock = collections.namedtuple("CellBlock", ["type", "data"])
[docs]
class Mesh(object):
def __init__(
self,
points,
cells,
point_data=None,
cell_data=None,
field_data=None,
point_sets=None,
cell_sets=None,
):
"""
Unstructured toughio mesh.
This class is updated following the latest :mod:`meshio` version and
brings backward compatibility with its previous versions.
Parameters
----------
points : array_like (n_points, 3)
Cooordinates of points.
cells : list of tuples (cell_type, data)
Connectivity of cells.
point_data : dict or None, optional, default None
Point data arrays.
cell_data : dict or None, optional, default None
Cell data arrays.
field_data : dict or None, optional, default None
Field data names.
point_sets : dict or None, optional, default None
Sets of points.
cell_sets : dict or None, optional, default None
Sets of cells.
"""
self.points = points
self.cells = cells
self.point_data = point_data if point_data else {}
self.cell_data = cell_data if cell_data else {}
self.field_data = field_data if field_data else {}
self.point_sets = point_sets if point_sets else {}
self.cell_sets = cell_sets if cell_sets else {}
self.label_length = None
def __repr__(self):
"""Represent a :class:`toughio.Mesh`."""
lines = [
"<toughio mesh object>",
f" Number of points: {len(self.points)}",
]
if len(self.cells) > 0:
lines.append(" Number of cells:")
for tpe, elems in self.cells:
lines.append(f" {tpe}: {len(elems)}")
else:
lines.append(" No cells.")
if self.point_sets:
lines.append(f" Point sets: {', '.join(self.point_sets)}")
if self.point_data:
lines.append(f" Point data: {', '.join(self.point_data)}")
if self.cell_data:
lines.append(f" Cell data: {', '.join(self.cell_data)}")
return "\n".join(lines)
def __getitem__(self, islice):
"""Slice mesh."""
if np.ndim(islice) == 0:
islice = [islice]
if np.ndim(islice) != 1:
raise ValueError()
if isinstance(islice[0], (bool, np.bool_)):
if len(islice) != self.n_cells:
raise ValueError()
elif isinstance(islice[0], (int, np.int8, np.int16, np.int32, np.int64)):
if np.max(islice) >= self.n_cells:
raise ValueError()
tmp = islice
islice = np.zeros(self.n_cells, dtype=bool)
islice[tmp] = True
else:
raise TypeError()
count = 0
cell_idx = []
nodes_to_keep = set()
for _, cell_data in self.cells:
cell_idx.append([])
for j, cell in enumerate(cell_data):
if islice[count]:
cell_idx[-1].append(j)
for node in cell:
nodes_to_keep.add(node)
count += 1
# Prune orphaned nodes
node_map = {}
point_idx = []
count = 0
for i in range(self.n_points):
if i in nodes_to_keep:
node_map[i] = count
point_idx.append(i)
count += 1
return Mesh(
points=self.points[point_idx],
cells=[
(
cell_type,
np.array([[node_map[i] for i in cell] for cell in cell_data[idx]]),
)
for idx, (cell_type, cell_data) in zip(cell_idx, self.cells)
],
point_data={k: v[point_idx] for k, v in self.point_data.items()},
cell_data={k: v[islice] for k, v in self.cell_data.items()},
field_data={k: v for k, v in self.field_data.items()},
)
[docs]
def extrude_to_3d(self, height=1.0, axis=2, inplace=True):
"""
Convert a 2D mesh to 3D by extruding cells along given axis.
Parameters
----------
height : scalar or array_like, optional, default 1.0
Height of extrusion.
axis : int (0, 1 or 2), optional, default 2
Axis along which extrusion is performed.
inplace : bool, optional, default True
If `False`, return a new :class:`toughio.Mesh`.
Returns
-------
toughio.Mesh
Extruded mesh (only if ``inplace == False``).
"""
if axis not in [0, 1, 2]:
raise ValueError("axis must be 0, 1 or 2.")
mesh = self if inplace else deepcopy(self)
height = [height] if np.ndim(height) == 0 else height
npts, nh = len(mesh.points), len(height)
if mesh.points.shape[1] == 3:
if len(set(mesh.points[:, axis])) != 1:
raise ValueError(f"Cannot extrude mesh along axis {axis}.")
else:
mesh.points = np.column_stack((mesh.points, np.zeros(npts)))
if axis != 2:
mesh.points[:, [axis, 2]] = mesh.points[:, [2, axis]]
extra_points = np.array(mesh.points)
for h in height:
extra_points[:, axis] += h
mesh.points = np.vstack((mesh.points, extra_points))
for k, v in mesh.point_data.items():
mesh.point_data[k] = (
np.tile(v, nh + 1) if np.ndim(v) == 1 else np.tile(v, (nh + 1, 1))
)
extruded_types = {
"triangle": "wedge",
"quad": "hexahedron",
}
cells = []
cell_data = {k: mesh.split(v) for k, v in mesh.cell_data.items()}
for ic, c in enumerate(mesh.cells):
if c.type in extruded_types:
extruded_type = extruded_types[c.type]
nr, nc = c.data.shape
cell = CellBlock(extruded_type, np.tile(c.data, (nh, 2)))
for i in range(nh):
ibeg, iend = i * nr, (i + 1) * nr
cell.data[ibeg:iend, :nc] += i * npts
cell.data[ibeg:iend, nc:] += (i + 1) * npts
cells.append(cell)
for k, v in cell_data.items():
v[ic] = (
np.tile(v[ic], nh)
if np.ndim(v[ic]) == 1
else np.tile(v[ic], (nh, 1))
)
mesh.cells = cells
mesh.cell_data = {k: np.concatenate(v) for k, v in cell_data.items()}
if mesh.field_data:
for k in mesh.field_data:
mesh.field_data[k][1] = 3
if not inplace:
return mesh
[docs]
def prune_duplicates(self, inplace=True):
"""
Delete duplicate points and cells.
Parameters
----------
inplace : bool, optional, default True
If `False`, return a new :class:`toughio.Mesh`.
Returns
-------
toughio.Mesh
Pruned mesh (only if ``inplace == False``).
Note
----
Does not preserve points order from original array in mesh.
"""
mesh = self if inplace else deepcopy(self)
cells = [[c.type, c.data] for c in mesh.cells]
# Prune duplicate points
unique_points, pind, pinv = np.unique(
mesh.points,
axis=0,
return_index=True,
return_inverse=True,
)
if len(unique_points) < len(mesh.points):
mesh.points = unique_points
for k, v in mesh.point_data.items():
mesh.point_data[k] = v[pind]
for ic, (k, v) in enumerate(cells):
cells[ic][1] = pinv[v]
# Prune duplicate cells
cell_data = {k: mesh.split(v) for k, v in mesh.cell_data.items()}
for ic, (k, v) in enumerate(cells):
vsort = np.sort(v, axis=1)
_, order = np.unique(vsort, axis=0, return_index=True)
cells[ic][1] = v[order]
for kk, vv in cell_data.items():
cell_data[kk][ic] = vv[ic][order]
mesh.cells = cells
mesh.cell_data = {k: np.concatenate(v) for k, v in cell_data.items()}
if not inplace:
return mesh
[docs]
def split(self, arr):
"""
Split input array into subarrays for each cell block in mesh.
Parameters
----------
arr : array_like
Input array.
Returns
-------
list of array_like
List of subarrays.
"""
if len(arr) != self.n_cells:
raise ValueError()
sizes = np.cumsum([len(c.data) for c in self.cells])
return np.split(np.asarray(arr), sizes[:-1])
[docs]
def to_meshio(self):
"""
Convert mesh to :class:`meshio.Mesh`.
Returns
-------
meshio.Mesh
Output mesh.
"""
keys = ["points", "point_data", "field_data"]
kwargs = {key: getattr(self, key) for key in keys}
cell_data = {k: self.split(v) for k, v in self.cell_data.items()}
kwargs.update(
{
"cells": self.cells,
"cell_data": cell_data,
"point_sets": self.point_sets,
"cell_sets": self.cell_sets,
}
)
return meshio.Mesh(**kwargs)
[docs]
def to_pyvista(self):
"""
Convert mesh to :class:`pyvista.UnstructuredGrid`.
Returns
-------
pyvista.UnstructuredGrid
Output mesh.
"""
try:
import pyvista
import vtk
from ._common import meshio_to_vtk_type, vtk_type_to_numnodes
VTK9 = vtk.vtkVersion().GetVTKMajorVersion() >= 9
except ImportError:
raise ImportError(
"Converting to pyvista.UnstructuredGrid requires pyvista to be installed."
)
# Extract cells from toughio.Mesh object
offset = []
cells = []
cell_type = []
next_offset = 0
for c in self.cells:
vtk_type = meshio_to_vtk_type[c.type]
numnodes = vtk_type_to_numnodes[vtk_type]
cells.append(
np.hstack((np.full((len(c.data), 1), numnodes), c.data)).ravel()
)
cell_type += [vtk_type] * len(c.data)
if not VTK9:
offset += [next_offset + i * (numnodes + 1) for i in range(len(c.data))]
next_offset = offset[-1] + numnodes + 1
# Create pyvista.UnstructuredGrid object
points = self.points
if points.shape[1] == 2:
points = np.hstack((points, np.zeros((len(points), 1))))
if VTK9:
mesh = pyvista.UnstructuredGrid(
np.concatenate(cells).astype(np.int64, copy=False),
np.array(cell_type),
np.array(points, np.float64),
)
else:
mesh = pyvista.UnstructuredGrid(
np.array(offset),
np.concatenate(cells).astype(np.int64, copy=False),
np.array(cell_type),
np.array(points, np.float64),
)
# Set point data
mesh.point_data.update(
{k: np.array(v, np.float64) for k, v in self.point_data.items()}
)
# Set cell data
mesh.cell_data.update(self.cell_data)
return mesh
[docs]
def write_tough(self, filename="MESH", **kwargs):
"""
Write TOUGH `MESH` file.
Parameters
----------
filename : str, pathlike or buffer, optional, default 'MESH'
Output file name or buffer.
Other Parameters
----------------
nodal_distance : str ('line' or 'orthogonal'), optional, default 'line'
Method to calculate connection nodal distances:
- 'line': distance between node and common face along connecting
line (distance is not normal),
- 'orthogonal' : distance between node and its orthogonal
projection onto common face (shortest distance).
material_name : dict or None, default None
Rename cell material.
material_end : str, array_like or None, default None
Move cells to bottom of block 'ELEME' if their materials is in `material_end`.
incon : bool, optional, default False
If `True`, initial conditions will be written in file `INCON`.
eos : str or None, optional, default None
Equation of State.
gravity : array_like or None, optional, default None
Gravity direction vector.
"""
self.write(filename, file_format="tough", **kwargs)
[docs]
def write_incon(self, filename="INCON", eos=None):
"""
Write TOUGH `INCON` file.
Parameters
----------
filename : str, pathlike or buffer, optional, default 'INCON'
Output file name or buffer.
eos : str or None, optional, default None
Equation of State.
Note
----
Mostly useful to restart a simulation with other initial conditions but with the same
mesh.
"""
from .tough._tough import check_incon, init_incon, write_incon
primary_variables, porosities, permeabilities, phase_compositions = init_incon(
self
)
incon = check_incon(
True,
primary_variables,
porosities,
permeabilities,
phase_compositions,
self.n_cells,
eos,
)
if incon:
write_incon(
filename,
self.labels,
primary_variables,
porosities,
permeabilities,
phase_compositions,
eos,
)
[docs]
def read_output(
self, file_or_output, time_step=-1, labels_order=None, connection=False
):
"""
Import TOUGH results to the mesh.
Parameters
----------
file_or_output : str, pathlike, buffer, namedtuple or list of namedtuple
Input file name or buffer, or output data.
time_step : int, optional, default -1
Data for given time step to import. Default is last time step.
labels_order : list of array_like or None, optional, default None
List of labels. If None, output will be assumed ordered.
connection : bool, optional, default False
Only for standard TOUGH output file. If `True`, read data related to connections.
"""
from .. import read_output
from .._io.output._common import Output, reorder_labels
if not isinstance(time_step, int):
raise TypeError()
if isinstance(file_or_output, str):
out = read_output(file_or_output, connection=connection)
else:
out = file_or_output
if not isinstance(out, Output):
if not (-len(out) <= time_step < len(out)):
raise ValueError()
out = out[time_step]
if out.type == "element":
if labels_order is not None:
out = reorder_labels(out, labels_order)
self.cell_data.update(out.data)
elif out.type == "connection":
centers = self.centers
labels_map = {k: v for v, k in enumerate(self.labels)}
data = {
k: [[[0.0, 0.0, 0.0]] for _ in range(self.n_cells)] for k in out.data
}
for i, (label1, label2) in enumerate(out.labels):
i1, i2 = labels_map[label1], labels_map[label2]
line = centers[i1] - centers[i2]
line /= np.linalg.norm(line)
for k, v in out.data.items():
iv = i1 if v[i] > 0.0 else i2
data[k][iv].append(v[i] * line)
data = {
k: np.array([np.sum(vv, axis=0) for vv in v]) for k, v in data.items()
}
self.cell_data.update(data)
[docs]
def write(self, filename, file_format=None, **kwargs):
"""
Write mesh to file.
Parameters
----------
filename : str, pathlike or buffer
Output file name or buffer.
file_format : str or None, optional, default None
Output file format. If `None`, it will be guessed from file's
extension. To write TOUGH MESH, `file_format` must be specified
as 'tough' (no specific extension exists for TOUGH MESH).
Other Parameters
----------------
nodal_distance : str ('line' or 'orthogonal'), optional, default 'line'
Only if ``file_format = "tough"``. Method to calculate connection
nodal distances:
- 'line': distance between node and common face along connecting
line (distance is not normal),
- 'orthogonal' : distance between node and its orthogonal
projection onto common face (shortest distance).
material_name : dict or None, default None
Only if ``file_format = "tough"``. Rename cell material.
material_end : str, array_like or None, default None
Only if ``file_format = "tough"``. Move cells to bottom of block
'ELEME' if their materials is in `material_end`.
incon : bool, optional, default False
Only if ``file_format = "tough"``. If `True`, initial conditions will be written in file `INCON`.
protocol : integer, optional, default `pickle.HIGHEST_PROTOCOL`
Only if ``file_format = "pickle"``. :mod:`pickle` protocol version.
"""
from ._helpers import write
write(filename, self, file_format, **kwargs)
[docs]
def plot(self, *args, **kwargs):
"""Display mesh using :meth:`pyvista.UnstructuredGrid.plot`."""
mesh = self.to_pyvista()
mesh.plot(*args, **kwargs)
[docs]
def add_material(self, label, imat):
"""
Add a material name.
Parameters
----------
label : str
Material name.
imat : int
Material ID.
"""
self.field_data[label] = np.array([imat, 3])
[docs]
def add_point_data(self, label, data):
"""
Add a new point data array.
Parameters
----------
label : str
Point data array name.
data : array_like
Point data array.
"""
if not isinstance(label, str):
raise TypeError()
if not isinstance(data, (list, tuple, np.ndarray)):
raise TypeError()
if len(data) != self.n_points:
raise ValueError()
self.point_data[label] = np.asarray(data)
[docs]
def add_cell_data(self, label, data):
"""
Add a new cell data array.
Parameters
----------
label : str
Cell data array name.
data : array_like
Cell data array.
"""
if not isinstance(label, str):
raise TypeError()
if not isinstance(data, (list, tuple, np.ndarray)):
raise TypeError()
if len(data) != self.n_cells:
raise ValueError()
self.cell_data[label] = np.asarray(data)
[docs]
def set_cell_labels(self, labels):
"""
Set labels to cells.
Parameters
----------
labels : array_like
Cell labels.
Note
----
If labels are not set, a custom labeler with a naming convention different that of TOUGH is used.
"""
if not isinstance(labels, (list, tuple, np.ndarray)):
raise TypeError()
if len(labels) != self.n_cells:
raise ValueError()
if not all(isinstance(label, str) for label in labels):
raise ValueError()
self._labels = labels
[docs]
def set_material(self, material, cells):
"""
Set material to cells.
Parameters
----------
material : str
Material name.
cells : array_like
Indices of cells or array of booleans.
"""
ints = (int, np.int8, np.int16, np.int32, np.int64)
cond_int = all(isinstance(c, ints) for c in cells)
cond_bool = all(isinstance(c, (bool, np.bool_)) for c in cells)
if cond_int:
if np.min(cells) < 0 or np.max(cells) >= self.n_cells:
raise ValueError()
elif cond_bool:
if len(cells) != self.n_cells:
raise ValueError()
cells = np.arange(self.n_cells)[cells]
else:
raise TypeError()
if len(cells):
data = self.cell_data["material"]
imat = (
self.field_data[material][0]
if material in self.field_data
else data.max() + 1
)
data[cells] = imat
self.add_cell_data("material", data)
self.add_material(material, imat)
[docs]
def cell_data_to_point_data(self):
"""Interpolate cell data to point data."""
from ._common import interpolate_data
points = [[] for _ in range(self.n_points)]
weights = [[] for _ in range(self.n_points)]
i = 0
volumes = self.volumes
for cellblock in self.cells:
for cell in cellblock.data:
for c in cell:
points[c].append(i)
weights[c].append(volumes[i])
i += 1
point_data = interpolate_data(self._cell_data, points, weights)
self._point_data.update(point_data)
self._cell_data = {}
[docs]
def point_data_to_cell_data(self):
"""Interpolate point data to cell data."""
from ._common import interpolate_data
cells = [c for cell in self.cells for c in cell.data]
cell_data = interpolate_data(self._point_data, cells)
self._cell_data.update(cell_data)
self._point_data = {}
[docs]
def near(self, points):
"""
Return indices of cells nearest to query points.
Parameters
----------
points : array_like
Coordinates of points to query.
Returns
-------
int or array_like
Index of cell or indices of cells.
"""
def distance(point, points):
"""Distance between point to list of points."""
dp = points - point
return np.einsum("ij,ij->i", dp, dp)
ndim = np.ndim(points)
if ndim == 0:
raise TypeError()
elif ndim == 1:
points = np.array([points])
else:
points = np.asarray(points)
if points.shape[1] != self.points.shape[1]:
raise ValueError()
centers = self.centers
idx = np.argmin([distance(point, centers) for point in points], axis=1)
return idx[0] if ndim == 1 else idx
@property
def points(self):
"""Return coordinates of points."""
return self._points
@points.setter
def points(self, value):
self._points = value
@property
def cells(self):
"""Return connectivity of cells."""
return self._cells
@cells.setter
def cells(self, value):
self._cells = [
CellBlock(*c) if isinstance(c, (list, tuple)) else CellBlock(c.type, c.data)
for c in value
]
@property
def point_data(self):
"""Return point data arrays."""
return self._point_data
@point_data.setter
def point_data(self, value):
self._point_data = value
@property
def cell_data(self):
"""Return cell data arrays."""
return self._cell_data
@cell_data.setter
def cell_data(self, value):
self._cell_data = value
@property
def field_data(self):
"""Return field data names."""
return self._field_data
@field_data.setter
def field_data(self, value):
self._field_data = value
@property
def point_sets(self):
"""Return sets of points."""
return self._point_sets
@point_sets.setter
def point_sets(self, value):
self._point_sets = value
@property
def cell_sets(self):
"""Return sets of cells."""
return self._cell_sets
@cell_sets.setter
def cell_sets(self, value):
self._cell_sets = value
@property
def n_points(self):
"""Return number of points."""
return len(self.points)
@property
def n_cells(self):
"""Return number of cells."""
return sum(len(c.data) for c in self.cells)
@property
def label_length(self):
"""Return length of labels."""
return self._label_length
@label_length.setter
def label_length(self, value):
self._label_length = value
@property
def labels(self):
"""Return labels of cell in mesh."""
from ._common import labeler
return (
self._labels
if hasattr(self, "_labels") and self._labels
else labeler(self.n_cells, self.label_length)
)
@property
def centers(self):
"""Return node centers of cell in mesh."""
return np.concatenate([self.points[c.data].mean(axis=1) for c in self.cells])
@property
def materials(self):
"""Return materials of cell in mesh."""
return _materials(self)
@property
def faces(self):
"""Return connectivity of faces of cell in mesh."""
out = _faces(self)
arr = np.full((self.n_cells, 6, 4), -1)
for i, x in enumerate(out):
arr[i, : len(x[0]), : x[0].shape[1]] = x[0]
if len(x) > 1:
arr[i, len(x[0]) : len(x[0]) + len(x[1]), : x[1].shape[1]] = x[1]
return arr
@property
def face_normals(self):
"""Return normal vectors of faces in mesh."""
return _face_normals(self)
@property
def face_areas(self):
"""Return areas of faces in mesh."""
return _face_areas(self)
@property
def volumes(self):
"""Return volumes of cell in mesh."""
return _volumes(self)
@property
def connections(self):
"""
Return mesh connections.
Assume conformity and that points and cells are uniquely defined in mesh.
Note
----
Only for 3D meshes and first order cells.
"""
return _connections(self)
@property
def qualities(self):
"""
Return qualities of cells in mesh.
The quality of a cell is measured as the minimum cosine angle between the
connection line and the interface normal vectors.
"""
return np.array([np.min(out) for out in _qualities(self)])
@property
def dim(self):
"""Return mesh dimension."""
celltypes = set([cell.type for cell in self.cells])
if celltypes.intersection({"tetra", "pyramid", "wedge", "hexahedron"}):
return 3
elif celltypes.intersection({"quad", "triangle"}):
return 2
else:
raise ValueError()
@property
def filter(self):
"""
Filter mesh.
Parameters
----------
filter_ : str, optional, default 'box'
Filter method.
Returns
-------
array_like
Indices of cells filtered.
"""
return MeshFilter(self)
[docs]
def from_meshio(mesh, material="dfalt"):
"""
Convert a :class:`meshio.Mesh` to :class:`toughio.Mesh`.
Parameters
----------
mesh : meshio.Mesh
Input mesh.
material : str, optional, default 'dfalt'
Default material name.
Returns
-------
toughio.Mesh
Output mesh.
"""
from ._helpers import get_material_key
if not isinstance(mesh, meshio.Mesh):
raise TypeError()
if not isinstance(material, str):
raise TypeError()
if mesh.cell_data:
cells = mesh.cells
cell_data = mesh.cell_data
key = get_material_key(cell_data)
if key:
cell_data["material"] = cell_data.pop(key)
cell_data = {k: np.concatenate(v) for k, v in cell_data.items()}
else:
cells = mesh.cells
cell_data = {}
out = Mesh(
points=mesh.points,
cells=cells,
point_data=mesh.point_data,
cell_data=cell_data,
field_data=mesh.field_data,
point_sets=(
mesh.point_sets
if hasattr(mesh, "point_sets")
else mesh.node_sets
if hasattr(mesh, "node_sets")
else None
),
cell_sets=mesh.cell_sets if hasattr(mesh, "cell_sets") else None,
)
if "material" not in out.cell_data:
imat = (
np.max([v[0] for v in mesh.field_data.values() if v[1] == 3]) + 1
if mesh.field_data
else 1
)
out.cell_data["material"] = np.full(out.n_cells, imat, dtype=np.int64)
out.field_data[material] = np.array([imat, 3])
return out
[docs]
def from_pyvista(mesh, material="dfalt"):
"""
Convert a :class:`pyvista.UnstructuredGrid` to :class:`toughio.Mesh`.
Parameters
----------
mesh : pyvista.UnstructuredGrid
Input mesh.
material : str, optional, default 'dfalt'
Default material name.
Returns
-------
toughio.Mesh
Output mesh.
"""
try:
import pyvista
import vtk
from ._common import vtk_to_meshio_type
VTK9 = vtk.vtkVersion().GetVTKMajorVersion() >= 9
except ImportError:
raise ImportError(
"Converting pyvista.UnstructuredGrid requires pyvista to be installed."
)
if not isinstance(mesh, pyvista.UnstructuredGrid):
raise TypeError()
if not isinstance(material, str):
raise TypeError()
# Copy useful arrays to avoid repeated calls to properties
vtk_offset = mesh.offset
vtk_cells = mesh.cells
vtk_cell_type = mesh.celltypes
# Check that meshio supports all cell types in input mesh
pixel_voxel = {8, 11} # Handle pixels and voxels
for cell_type in np.unique(vtk_cell_type):
if not (cell_type in vtk_to_meshio_type or cell_type in pixel_voxel):
raise ValueError(f"toughio does not support VTK type {cell_type}.")
# Get cells
cells = []
c = 0
for offset, cell_type in zip(vtk_offset, vtk_cell_type):
numnodes = vtk_cells[offset]
if VTK9:
cell = vtk_cells[offset + 1 + c : offset + 1 + c + numnodes]
c += 1
else:
cell = vtk_cells[offset + 1 : offset + 1 + numnodes]
cell = (
cell
if cell_type not in pixel_voxel
else cell[[0, 1, 3, 2]]
if cell_type == 8
else cell[[0, 1, 3, 2, 4, 5, 7, 6]]
)
cell_type = cell_type if cell_type not in pixel_voxel else cell_type + 1
cell_type = (
vtk_to_meshio_type[cell_type] if cell_type != 7 else f"polygon{numnodes}"
)
if len(cells) > 0 and cells[-1][0] == cell_type:
cells[-1][1].append(cell)
else:
cells.append((cell_type, [cell]))
for k, c in enumerate(cells):
cells[k] = (c[0], np.array(c[1]))
# Get point data
point_data = {k.replace(" ", "_"): v for k, v in mesh.point_data.items()}
# Get cell data
cell_data = {k.replace(" ", "_"): v for k, v in mesh.cell_data.items()}
# Create toughio.Mesh
out = Mesh(
points=np.array(mesh.points),
cells=cells,
point_data=point_data,
cell_data=cell_data,
)
if "material" not in out.cell_data:
imat = 1
out.cell_data["material"] = np.full(out.n_cells, imat, dtype=np.int64)
out.field_data[material] = np.array([imat, 3])
return out
