from __future__ import annotations import copy import numpy as np from numpy.typing import ArrayLike from ._common import num_nodes_per_cell, warn topological_dimension = { "line": 1, "polygon": 2, "triangle": 2, "quad": 2, "tetra": 3, "hexahedron": 3, "wedge": 3, "pyramid": 3, "line3": 1, "triangle6": 2, "quad9": 2, "tetra10": 3, "hexahedron27": 3, "wedge18": 3, "pyramid14": 3, "vertex": 0, "quad8": 2, "hexahedron20": 3, "triangle10": 2, "triangle15": 2, "triangle21": 2, "line4": 1, "line5": 1, "line6": 1, "tetra20": 3, "tetra35": 3, "tetra56": 3, "quad16": 2, "quad25": 2, "quad36": 2, "triangle28": 2, "triangle36": 2, "triangle45": 2, "triangle55": 2, "triangle66": 2, "quad49": 2, "quad64": 2, "quad81": 2, "quad100": 2, "quad121": 2, "line7": 1, "line8": 1, "line9": 1, "line10": 1, "line11": 1, "tetra84": 3, "tetra120": 3, "tetra165": 3, "tetra220": 3, "tetra286": 3, "wedge40": 3, "wedge75": 3, "hexahedron64": 3, "hexahedron125": 3, "hexahedron216": 3, "hexahedron343": 3, "hexahedron512": 3, "hexahedron729": 3, "hexahedron1000": 3, "wedge126": 3, "wedge196": 3, "wedge288": 3, "wedge405": 3, "wedge550": 3, "VTK_LAGRANGE_CURVE": 1, "VTK_LAGRANGE_TRIANGLE": 2, "VTK_LAGRANGE_QUADRILATERAL": 2, "VTK_LAGRANGE_TETRAHEDRON": 3, "VTK_LAGRANGE_HEXAHEDRON": 3, "VTK_LAGRANGE_WEDGE": 3, "VTK_LAGRANGE_PYRAMID": 3, } class CellBlock: def __init__( self, cell_type: str, data: list | np.ndarray, tags: list[str] | None = None, ): self.type = cell_type self.data = data if cell_type.startswith("polyhedron"): self.dim = 3 else: self.data = np.asarray(self.data) self.dim = topological_dimension[cell_type] self.tags = [] if tags is None else tags def __repr__(self): items = [ "meshio CellBlock", f"type: {self.type}", f"num cells: {len(self.data)}", f"tags: {self.tags}", ] return "<" + ", ".join(items) + ">" def __len__(self): return len(self.data) class Mesh: def __init__( self, points: ArrayLike, cells: dict[str, ArrayLike] | list[tuple[str, ArrayLike] | CellBlock], point_data: dict[str, ArrayLike] | None = None, cell_data: dict[str, list[ArrayLike]] | None = None, field_data=None, point_sets: dict[str, ArrayLike] | None = None, cell_sets: dict[str, list[ArrayLike]] | None = None, gmsh_periodic=None, info=None, ): self.points = np.asarray(points) if isinstance(cells, dict): # Let's not deprecate this for now. # warn( # "cell dictionaries are deprecated, use list of tuples, e.g., " # '[("triangle", [[0, 1, 2], ...])]', # DeprecationWarning, # ) # old dict, deprecated # # convert dict to list of tuples cells = list(cells.items()) self.cells = [] for cell_block in cells: if isinstance(cell_block, tuple): cell_type, data = cell_block cell_block = CellBlock( cell_type, # polyhedron data cannot be converted to numpy arrays # because the sublists don't all have the same length data if cell_type.startswith("polyhedron") else np.asarray(data), ) self.cells.append(cell_block) self.point_data = {} if point_data is None else point_data self.cell_data = {} if cell_data is None else cell_data self.field_data = {} if field_data is None else field_data self.point_sets = {} if point_sets is None else point_sets self.cell_sets = {} if cell_sets is None else cell_sets self.gmsh_periodic = gmsh_periodic self.info = info # assert point data consistency and convert to numpy arrays for key, item in self.point_data.items(): self.point_data[key] = np.asarray(item) if len(self.point_data[key]) != len(self.points): raise ValueError( f"len(points) = {len(self.points)}, " f'but len(point_data["{key}"]) = {len(self.point_data[key])}' ) # assert cell data consistency and convert to numpy arrays for key, data in self.cell_data.items(): if len(data) != len(cells): raise ValueError( f"Incompatible cell data '{key}'. " f"{len(cells)} cell blocks, but '{key}' has {len(data)} blocks." ) for k in range(len(data)): data[k] = np.asarray(data[k]) if len(data[k]) != len(self.cells[k]): raise ValueError( "Incompatible cell data. " + f"Cell block {k} ('{self.cells[k].type}') " + f"has length {len(self.cells[k])}, but " + f"corresponding cell data item has length {len(data[k])}." ) def __repr__(self): lines = ["", f" Number of points: {len(self.points)}"] special_cells = [ "polygon", "polyhedron", "VTK_LAGRANGE_CURVE", "VTK_LAGRANGE_TRIANGLE", "VTK_LAGRANGE_QUADRILATERAL", "VTK_LAGRANGE_TETRAHEDRON", "VTK_LAGRANGE_HEXAHEDRON", "VTK_LAGRANGE_WEDGE", "VTK_LAGRANGE_PYRAMID", ] if len(self.cells) > 0: lines.append(" Number of cells:") for cell_block in self.cells: string = cell_block.type if cell_block.type in special_cells: string += f"({cell_block.data.shape[1]})" lines.append(f" {string}: {len(cell_block)}") else: lines.append(" No cells.") if self.point_sets: names = ", ".join(self.point_sets.keys()) lines.append(f" Point sets: {names}") if self.cell_sets: names = ", ".join(self.cell_sets.keys()) lines.append(f" Cell sets: {names}") if self.point_data: names = ", ".join(self.point_data.keys()) lines.append(f" Point data: {names}") if self.cell_data: names = ", ".join(self.cell_data.keys()) lines.append(f" Cell data: {names}") if self.field_data: names = ", ".join(self.field_data.keys()) lines.append(f" Field data: {names}") return "\n".join(lines) def copy(self): return copy.deepcopy(self) def write(self, path_or_buf, file_format: str | None = None, **kwargs): # avoid circular import from ._helpers import write write(path_or_buf, self, file_format, **kwargs) def get_cells_type(self, cell_type: str): if not any(c.type == cell_type for c in self.cells): return np.empty((0, num_nodes_per_cell[cell_type]), dtype=int) return np.concatenate([c.data for c in self.cells if c.type == cell_type]) def get_cell_data(self, name: str, cell_type: str): return np.concatenate( [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type] ) @property def cells_dict(self): cells_dict = {} for cell_block in self.cells: if cell_block.type not in cells_dict: cells_dict[cell_block.type] = [] cells_dict[cell_block.type].append(cell_block.data) # concatenate for key, value in cells_dict.items(): cells_dict[key] = np.concatenate(value) return cells_dict @property def cell_data_dict(self): cell_data_dict = {} for key, value_list in self.cell_data.items(): cell_data_dict[key] = {} for value, cell_block in zip(value_list, self.cells): if cell_block.type not in cell_data_dict[key]: cell_data_dict[key][cell_block.type] = [] cell_data_dict[key][cell_block.type].append(value) for cell_type, val in cell_data_dict[key].items(): cell_data_dict[key][cell_type] = np.concatenate(val) return cell_data_dict @property def cell_sets_dict(self): sets_dict = {} for key, member_list in self.cell_sets.items(): sets_dict[key] = {} offsets = {} for members, cells in zip(member_list, self.cells): if members is None: continue if cells.type in offsets: offset = offsets[cells.type] offsets[cells.type] += cells.data.shape[0] else: offset = 0 offsets[cells.type] = cells.data.shape[0] if cells.type in sets_dict[key]: sets_dict[key][cells.type].append(members + offset) else: sets_dict[key][cells.type] = [members + offset] return { key: { cell_type: np.concatenate(members) for cell_type, members in sets.items() if sum(map(np.size, members)) } for key, sets in sets_dict.items() } @classmethod def read(cls, path_or_buf, file_format=None): # avoid circular import from ._helpers import read # 2021-02-21 warn("meshio.Mesh.read is deprecated, use meshio.read instead") return read(path_or_buf, file_format) def cell_sets_to_data(self, data_name: str | None = None): # If possible, convert cell sets to integer cell data. This is possible if all # cells appear exactly in one group. default_value = -1 if len(self.cell_sets) > 0: intfun = [] for k, c in enumerate(zip(*self.cell_sets.values())): # Go for -1 as the default value. (NaN is not int.) arr = np.full(len(self.cells[k]), default_value, dtype=int) for i, cc in enumerate(c): if cc is None: continue arr[cc] = i intfun.append(arr) for item in intfun: num_default = np.sum(item == default_value) if num_default > 0: warn( f"{num_default} cells are not part of any cell set. " f"Using default value {default_value}." ) break if data_name is None: data_name = "-".join(self.cell_sets.keys()) self.cell_data[data_name] = intfun self.cell_sets = {} def point_sets_to_data(self, join_char: str = "-") -> None: # now for the point sets # Go for -1 as the default value. (NaN is not int.) default_value = -1 if len(self.point_sets) > 0: intfun = np.full(len(self.points), default_value, dtype=int) for i, cc in enumerate(self.point_sets.values()): intfun[cc] = i if np.any(intfun == default_value): warn( "Not all points are part of a point set. " f"Using default value {default_value}." ) data_name = join_char.join(self.point_sets.keys()) self.point_data[data_name] = intfun self.point_sets = {} # This used to be int_data_to_sets(), converting _all_ cell and point data. # This is not useful in many cases, as one usually only wants one # particular data array (e.g., "MaterialIDs") converted to sets. def cell_data_to_sets(self, key: str): """Convert point_data to cell_sets.""" data = self.cell_data[key] # handle all int and uint data if not all(v.dtype.kind in ["i", "u"] for v in data): raise RuntimeError(f"cell_data['{key}'] is not int data.") tags = np.unique(np.concatenate(data)) # try and get the names by splitting the key along "-" (this is how # sets_to_int_data() forms the key) names = key.split("-") # remove duplicates and preserve order # : names = list(dict.fromkeys(names)) if len(names) != len(tags): # alternative names names = [f"set-{key}-{tag}" for tag in tags] # TODO there's probably a better way besides np.where, something from # np.unique or np.sort for name, tag in zip(names, tags): self.cell_sets[name] = [np.where(d == tag)[0] for d in data] # remove the cell data del self.cell_data[key] def point_data_to_sets(self, key: str): """Convert point_data to point_sets.""" data = self.point_data[key] # handle all int and uint data if not all(v.dtype.kind in ["i", "u"] for v in data): raise RuntimeError(f"point_data['{key}'] is not int data.") tags = np.unique(data) # try and get the names by splitting the key along "-" (this is how # sets_to_int_data() forms the key names = key.split("-") # remove duplicates and preserve order # : names = list(dict.fromkeys(names)) if len(names) != len(tags): # alternative names names = [f"set-key-{tag}" for tag in tags] # TODO there's probably a better way besides np.where, something from # np.unique or np.sort for name, tag in zip(names, tags): self.point_sets[name] = np.where(data == tag)[0] # remove the cell data del self.point_data[key]