Coverage for src/beamme/core/mesh_representation.py: 95%
178 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-05 08:51 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2026 BeamMe Authors
4#
5# Permission is hereby granted, free of charge, to any person obtaining a copy
6# of this software and associated documentation files (the "Software"), to deal
7# in the Software without restriction, including without limitation the rights
8# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9# copies of the Software, and to permit persons to whom the Software is
10# furnished to do so, subject to the following conditions:
11#
12# The above copyright notice and this permission notice shall be included in
13# all copies or substantial portions of the Software.
14#
15# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21# THE SOFTWARE.
22"""This module defines the mesh representation data structure."""
24from collections.abc import Iterable as _Iterable
25from dataclasses import dataclass as _dataclass
26from itertools import repeat as _repeat
27from typing import Any as _Any
29import numpy as _np
30import pyvista as _pv
31from numpy.typing import NDArray as _NDArray
33from beamme.core.conf import Geometry as _Geometry
34from beamme.core.conf import bme as _bme
36MESH_REPRESENTATION_MAPPINGS: dict[str, _Any] = {}
37# fmt: off
38MESH_REPRESENTATION_MAPPINGS[
39 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
40] = {
42 # Only list the non-standard mappings
43 (_bme.element_type.solid, 20):
44 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15],
45 (_bme.element_type.solid, 27):
46 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15,
47 24, 22, 21, 23, 20, 25, 26],
48}
49# fmt: on
50MESH_REPRESENTATION_MAPPINGS[
51 "element_type_and_n_nodes_to_connectivity_mapping_vtk_to_beamme"
52] = {
53 # Only list the non-standard mappings
54 (_bme.element_type.solid, 20): _np.argsort(
55 MESH_REPRESENTATION_MAPPINGS[
56 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
57 ][(_bme.element_type.solid, 20)]
58 ),
59 (_bme.element_type.solid, 27): _np.argsort(
60 MESH_REPRESENTATION_MAPPINGS[
61 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
62 ][(_bme.element_type.solid, 27)]
63 ),
64}
66# fmt: off
67MESH_REPRESENTATION_MAPPINGS["connectivity_mapping_exodus_to_vtk"] = {
68 # Only list the non-standard mappings
69 20: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15],
70 27: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15, 23, 24, 25, 26, 21, 22, 20],
71}
72# fmt: on
75GEOMETRY_SET_INFO_PREFIX = "geometry_set"
78@_dataclass
79class GeometrySetInfo:
80 """Data class that contains information of a geometry set."""
82 geometry_type: _Geometry
83 i_global: int
84 name: str | None = None
85 point_flag_vector: _NDArray | None = None
86 cell_flag_vector: _NDArray | None = None
88 def __str__(self):
89 """Return a string representation of this geometry set info - the node and cell flags are not included."""
90 name = self.name
91 return (
92 f"{GEOMETRY_SET_INFO_PREFIX}_{self.i_global}_{self.geometry_type.name}"
93 + (f"_{name}" if name is not None else "")
94 )
97def string_to_geometry_set_info(name: str) -> GeometrySetInfo | None:
98 """Extract the geometry set information from a given string."""
99 if not name.startswith(GEOMETRY_SET_INFO_PREFIX):
100 return None
102 name_without_prefix = name[len(GEOMETRY_SET_INFO_PREFIX) + 1 :]
103 split = name_without_prefix.split("_", 2)
104 i_global = int(split[0])
105 geometry_type = _Geometry[split[1]]
106 return GeometrySetInfo(
107 i_global=i_global,
108 geometry_type=geometry_type,
109 name=split[2] if len(split) == 3 else None,
110 )
113class MeshRepresentation:
114 """Class representing a generic mesh."""
116 def __init__(
117 self,
118 cell_connectivity: _NDArray[_np.integer] | None = None,
119 cell_types: _NDArray[_np.integer] | None = None,
120 points: _NDArray[_np.floating] | None = None,
121 geometry_sets: list[GeometrySetInfo] | None = None,
122 cell_data: dict[str, _NDArray | None] | None = None,
123 point_data: dict[str, _NDArray | None] | None = None,
124 ):
125 def _convert_argument_numpy(
126 argument: _NDArray | None,
127 default_shape: tuple[int, ...],
128 dtype: type,
129 ) -> _NDArray:
130 """Convert given array arguments so we can store them in this object."""
131 if argument is None:
132 return _np.empty(default_shape, dtype=dtype)
133 else:
134 return _np.asarray(argument, dtype=dtype)
136 def _filter_none_entries(
137 argument: dict[str, _NDArray | None] | None, expected_size: int
138 ) -> dict[str, _NDArray]:
139 """Check if a dictionary is given, and if so, filter None entries from
140 it."""
141 if argument is None:
142 return {}
143 else:
144 data = {
145 key: _np.asarray(value)
146 for key, value in argument.items()
147 if value is not None
148 }
149 for key, value in data.items():
150 if len(value) != expected_size:
151 raise ValueError(
152 f"Data field {key} has size {len(value)}, but expected size is {expected_size}."
153 )
154 return data
156 self.cell_connectivity = _convert_argument_numpy(
157 cell_connectivity, default_shape=(0,), dtype=int
158 )
159 self.cell_types = _convert_argument_numpy(
160 cell_types, default_shape=(0,), dtype=int
161 )
162 self.points = _convert_argument_numpy(points, default_shape=(0, 3), dtype=float)
164 self.cell_data = _filter_none_entries(cell_data, self.n_cells)
165 self.point_data = _filter_none_entries(point_data, self.n_points)
167 def _store_geometry_set_data(
168 data_field: str,
169 name: str,
170 field_vector: _NDArray | None,
171 expected_size: int,
172 ) -> None:
173 """Store point or cell data defining a geometry set."""
174 if field_vector is not None:
175 if len(field_vector) != expected_size:
176 raise ValueError(
177 f"Geometry set {name} has size {len(field_vector)},"
178 f" but expected size is {expected_size}."
179 )
180 getattr(self, data_field)[str(name)] = field_vector
182 if geometry_sets is not None:
183 for geometry_set in geometry_sets:
184 geometry_set_name = str(geometry_set)
185 _store_geometry_set_data(
186 "point_data",
187 geometry_set_name,
188 geometry_set.point_flag_vector,
189 self.n_points,
190 )
191 _store_geometry_set_data(
192 "cell_data",
193 geometry_set_name,
194 geometry_set.cell_flag_vector,
195 self.n_cells,
196 )
198 # Get the offset array so we can iterate over the connectivity in a
199 # performant manner.
200 self.cell_connectivity_offsets = _np.zeros(self.n_cells + 1, dtype=int)
201 for i_cell in range(self.n_cells):
202 last_offset = self.cell_connectivity_offsets[i_cell]
203 if last_offset >= len(self.cell_connectivity):
204 raise ValueError(
205 "Invalid offset, is larger than the size of the connectivity."
206 )
207 cell_size = self.cell_connectivity[last_offset]
208 if cell_size < 1:
209 raise ValueError("Invalid offset, has to be at least 1.")
210 self.cell_connectivity_offsets[i_cell + 1] = (
211 cell_size + self.cell_connectivity_offsets[i_cell] + 1
212 )
213 if self.cell_connectivity_offsets[-1] != len(self.cell_connectivity):
214 raise ValueError("Invalid cell connectivity offsets.")
216 @property
217 def n_cells(self) -> int:
218 """Return the number of cells in this mesh representation.
220 Returns:
221 Number of cells in this mesh representation.
222 """
223 return len(self.cell_types)
225 @property
226 def n_points(self) -> int:
227 """Return the number of points in this mesh representation.
229 Returns:
230 Number of points in this mesh representation.
231 """
232 return len(self.points)
234 def connectivity_iterator(
235 self, *, element_indices: list[int] | _NDArray[_np.integer] | None = None
236 ) -> _Iterable:
237 """Return an iterator over the cell connectivity.
239 Args:
240 element_indices: If this argument is given, only iterate over the
241 given indices.
243 Returns:
244 An iterator that returns the connectivity array for each cell.
245 """
246 if element_indices is None:
247 element_indices = _np.arange(self.n_cells)
248 for i in element_indices:
249 start = self.cell_connectivity_offsets[i] + 1
250 end = self.cell_connectivity_offsets[i + 1]
251 yield self.cell_connectivity[start:end]
253 def data_iterator(self, data_field: str, data_name: str) -> _Iterable:
254 """This method returns an iterator for the given data field and data name.
256 This is useful, when looping over the data, as accessing the data field in each
257 loop iteration can be expensive. If the data field is not present, a iterator
258 is returned that will always return `None`.
260 Args:
261 data_field: The data field to get the iterator for. This can be either
262 "point_data" or "cell_data".
263 data_name: The name of the data to get the iterator for.
265 Returns:
266 An iterator to iterate through the data. If the data field is not present,
267 a iterator is returned that will always return `None`.
268 """
269 data_dict = getattr(self, data_field)
270 if data_name in data_dict:
271 return data_dict[data_name]
272 else:
273 match data_field:
274 case "point_data":
275 size = self.n_points
276 case "cell_data":
277 size = self.n_cells
278 case _:
279 raise ValueError(f"Invalid data field: {data_field}")
280 return _repeat(None, size)
282 def offset_indices(
283 self,
284 element_type_id_offset: int | None = None,
285 material_offset: int | None = None,
286 geometry_set_offset: int | None = None,
287 ) -> None:
288 """Add offsets to the internal indices of this mesh representation.
290 This is required when merging multiple mesh representations together,
291 to ensure that there are no index conflicts.
293 Args:
294 element_type_id_offset: The offset to add to the element type IDs.
295 material_offset: The offset to add to the material IDs.
296 geometry_set_offset: The offset to add to the geometry set IDs.
297 """
298 if element_type_id_offset is not None:
299 if "element_type_id" in self.cell_data:
300 self.cell_data["element_type_id"] += element_type_id_offset
302 if material_offset is not None:
303 if "material_id" in self.cell_data:
304 # Add the offset to all material entries that are not -1. We use
305 # -1 to indicate no assigned material and those values should not change.
306 self.cell_data["material_id"][self.cell_data["material_id"] >= 0] += (
307 material_offset
308 )
310 if geometry_set_offset is not None:
311 for field_type in ("point_data", "cell_data"):
312 new_dict = {}
313 data_field = getattr(self, field_type)
314 # We change the keys in the loop, so we iterate over a copy of them.
315 for name in list(data_field.keys()):
316 info = string_to_geometry_set_info(name)
317 if info is not None:
318 new_name = str(
319 GeometrySetInfo(
320 geometry_type=info.geometry_type,
321 i_global=info.i_global + geometry_set_offset,
322 name=info.name,
323 )
324 )
325 new_dict[new_name] = data_field.pop(name)
326 # Add the renamed geometry sets back to the data field.
327 setattr(self, field_type, {**data_field, **new_dict})
329 def get_geometry_set_infos(self) -> list[GeometrySetInfo]:
330 """Get a list of all the geometry set info objects contained in this mesh
331 representation.
333 Returns:
334 A list of GeometrySetInfo objects. A reference to the internal data is added to the info object.
335 """
336 geometry_sets = []
337 for field_type, flag_vector_name in (
338 ("point_data", "point_flag_vector"),
339 ("cell_data", "cell_flag_vector"),
340 ):
341 data_field = getattr(self, field_type)
342 for name in data_field.keys():
343 info = string_to_geometry_set_info(name)
344 if info is not None:
345 # Set the flag vector in the info object.
346 setattr(info, flag_vector_name, data_field[name])
347 geometry_sets.append(info)
348 return geometry_sets
350 def get_pyvista_grid(
351 self,
352 *,
353 cell_data_fields: bool | list[str] = False,
354 point_data_fields: bool | list[str] = False,
355 add_geometry_sets: bool = False,
356 ) -> _pv.UnstructuredGrid:
357 """Return a PyVista UnstructuredGrid representation of this mesh representation.
359 Args:
360 cell_data_fields: The cell data fields to add to the grid. This can be
361 either a list of field names, or a boolean. If a list of field names
362 is given, only those fields are added to the grid. If it is True, all
363 cell data fields are added to the grid. If it is False, no cell data
364 fields are added to the grid.
365 point_data_fields: The point data fields to add to the grid. This can be
366 either a list of field names, or a boolean. If a list of field names
367 is given, only those fields are added to the grid. If it is True, all
368 point data fields are added to the grid. If it is False, no point data
369 fields are added to the grid.
370 add_geometry_sets: If this is True, all geometry set information is added to
371 the grid, even if it is not included in the cell_data_fields or
372 point_data_fields arguments.
374 Returns:
375 A PyVista UnstructuredGrid representation of this mesh representation.
376 """
377 grid = _pv.UnstructuredGrid(
378 self.cell_connectivity,
379 self.cell_types,
380 self.points,
381 )
383 for names, mesh_representation_field, grid_field in zip(
384 (cell_data_fields, point_data_fields),
385 (self.cell_data, self.point_data),
386 (grid.cell_data, grid.point_data),
387 ):
388 for name, data in mesh_representation_field.items():
389 add_field = False
390 if type(names) is bool:
391 add_field = names
392 elif name in names:
393 add_field = True
394 elif add_geometry_sets:
395 add_field = string_to_geometry_set_info(name) is not None
396 if add_field:
397 grid_field[name] = data
399 return grid
402def merge_mesh_representations(
403 mesh_representation_a: MeshRepresentation, mesh_representation_b: MeshRepresentation
404) -> MeshRepresentation:
405 """Merge two mesh representations.
407 Args:
408 mesh_representation_a: First mesh representation.
409 mesh_representation_b: Second mesh representation.
411 Returns:
412 A merged mesh representation. Depending on the input data, this might
413 be a reference to one of the input mesh representations.
414 """
416 def _is_empty(mesh_representation: MeshRepresentation) -> bool:
417 """Check if the mesh representation contains any points or cells."""
418 return mesh_representation.n_points == 0 and mesh_representation.n_cells == 0
420 # If one of the two mesh representations is empty, we can simply return the other (which
421 # could be empty as well).
422 if _is_empty(mesh_representation_a):
423 return mesh_representation_b
424 elif _is_empty(mesh_representation_b):
425 return mesh_representation_a
427 # Merge the points
428 merged_points = _np.vstack(
429 (mesh_representation_a.points, mesh_representation_b.points)
430 )
432 # Merge the connectivity
433 # Before the connectivity can be merged, the IDs of the second mesh have
434 # to be increased by the size of the first mesh. To do so, we create a
435 # mask to identify the entries that have to be incremented.
436 cell_connectivity_a = mesh_representation_a.cell_connectivity
437 cell_connectivity_b = mesh_representation_b.cell_connectivity
438 offsets_b = mesh_representation_b.cell_connectivity_offsets
440 merged_cell_connectivity = _np.empty(
441 cell_connectivity_a.size + cell_connectivity_b.size,
442 dtype=cell_connectivity_a.dtype,
443 )
444 merged_cell_connectivity[: cell_connectivity_a.size] = cell_connectivity_a
445 merged_cell_connectivity[cell_connectivity_a.size :] = cell_connectivity_b
447 mask = _np.ones(cell_connectivity_b.size, dtype=bool)
448 mask[offsets_b[:-1]] = False
450 merged_cell_connectivity_view_part_b = merged_cell_connectivity[
451 cell_connectivity_a.size :
452 ]
453 merged_cell_connectivity_view_part_b[mask] += mesh_representation_a.n_points
455 # Merge the cell types.
456 merged_cell_types = _np.concatenate(
457 (
458 mesh_representation_a.cell_types,
459 mesh_representation_b.cell_types,
460 )
461 )
463 # Merge the contained data dictionaries.
464 def _merge_data_dicts(
465 dict_a: dict[str, _NDArray],
466 size_a: int,
467 dict_b: dict[str, _NDArray],
468 size_b: int,
469 ) -> dict[str, _NDArray]:
470 """Merge the given data dictionaries and fill in non-existing data."""
472 def _ensure_array_size(size: int, reference_array: _NDArray) -> _NDArray:
473 """Create an empty array that matches the columns of the reference array and
474 has the given size, i.e., number of rows."""
475 new_shape = reference_array.shape
476 if len(new_shape) == 1:
477 new_shape = (size,)
478 elif len(new_shape) == 2:
479 new_shape = (size, new_shape[1])
480 else:
481 raise ValueError(f"Got unexpected array shape {new_shape}.")
482 return _np.zeros(new_shape, dtype=reference_array.dtype)
484 all_keys = set(dict_a.keys()) | set(dict_b.keys())
485 merged_data = {}
486 for key in all_keys:
487 data_a = dict_a.get(key, None)
488 data_b = dict_b.get(key, None)
489 if data_a is None:
490 data_a = _ensure_array_size(size_a, data_b)
491 elif data_b is None:
492 data_b = _ensure_array_size(size_b, data_a)
493 merged_data[key] = _np.concatenate((data_a, data_b))
494 return merged_data
496 merged_cell_data = _merge_data_dicts(
497 mesh_representation_a.cell_data,
498 mesh_representation_a.n_cells,
499 mesh_representation_b.cell_data,
500 mesh_representation_b.n_cells,
501 )
502 merged_point_data = _merge_data_dicts(
503 mesh_representation_a.point_data,
504 mesh_representation_a.n_points,
505 mesh_representation_b.point_data,
506 mesh_representation_b.n_points,
507 )
509 return MeshRepresentation(
510 cell_connectivity=merged_cell_connectivity,
511 cell_types=merged_cell_types,
512 points=merged_points,
513 cell_data=merged_cell_data,
514 point_data=merged_point_data,
515 )