Coverage for src/beamme/core/mesh_representation.py: 95%
168 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-28 15:20 +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_pyvista_grid(
330 self,
331 *,
332 cell_data_fields: bool | list[str] = False,
333 point_data_fields: bool | list[str] = False,
334 add_geometry_sets: bool = False,
335 ) -> _pv.UnstructuredGrid:
336 """Return a PyVista UnstructuredGrid representation of this mesh representation.
338 Args:
339 cell_data_fields: The cell data fields to add to the grid. This can be
340 either a list of field names, or a boolean. If a list of field names
341 is given, only those fields are added to the grid. If it is True, all
342 cell data fields are added to the grid. If it is False, no cell data
343 fields are added to the grid.
344 point_data_fields: The point data fields to add to the grid. This can be
345 either a list of field names, or a boolean. If a list of field names
346 is given, only those fields are added to the grid. If it is True, all
347 point data fields are added to the grid. If it is False, no point data
348 fields are added to the grid.
349 add_geometry_sets: If this is True, all geometry set information is added to
350 the grid, even if it is not included in the cell_data_fields or
351 point_data_fields arguments.
353 Returns:
354 A PyVista UnstructuredGrid representation of this mesh representation.
355 """
356 grid = _pv.UnstructuredGrid(
357 self.cell_connectivity,
358 self.cell_types,
359 self.points,
360 )
362 for names, mesh_representation_field, grid_field in zip(
363 (cell_data_fields, point_data_fields),
364 (self.cell_data, self.point_data),
365 (grid.cell_data, grid.point_data),
366 ):
367 for name, data in mesh_representation_field.items():
368 add_field = False
369 if type(names) is bool:
370 add_field = names
371 elif name in names:
372 add_field = True
373 elif add_geometry_sets:
374 add_field = string_to_geometry_set_info(name) is not None
375 if add_field:
376 grid_field[name] = data
378 return grid
381def merge_mesh_representations(
382 mesh_representation_a: MeshRepresentation, mesh_representation_b: MeshRepresentation
383) -> MeshRepresentation:
384 """Merge two mesh representations.
386 Args:
387 mesh_representation_a: First mesh representation.
388 mesh_representation_b: Second mesh representation.
390 Returns:
391 A merged mesh representation. Depending on the input data, this might
392 be a reference to one of the input mesh representations.
393 """
395 def _is_empty(mesh_representation: MeshRepresentation) -> bool:
396 """Check if the mesh representation contains any points or cells."""
397 return mesh_representation.n_points == 0 and mesh_representation.n_cells == 0
399 # If one of the two mesh representations is empty, we can simply return the other (which
400 # could be empty as well).
401 if _is_empty(mesh_representation_a):
402 return mesh_representation_b
403 elif _is_empty(mesh_representation_b):
404 return mesh_representation_a
406 # Merge the points
407 merged_points = _np.vstack(
408 (mesh_representation_a.points, mesh_representation_b.points)
409 )
411 # Merge the connectivity
412 # Before the connectivity can be merged, the IDs of the second mesh have
413 # to be increased by the size of the first mesh. To do so, we create a
414 # mask to identify the entries that have to be incremented.
415 cell_connectivity_a = mesh_representation_a.cell_connectivity
416 cell_connectivity_b = mesh_representation_b.cell_connectivity
417 offsets_b = mesh_representation_b.cell_connectivity_offsets
419 merged_cell_connectivity = _np.empty(
420 cell_connectivity_a.size + cell_connectivity_b.size,
421 dtype=cell_connectivity_a.dtype,
422 )
423 merged_cell_connectivity[: cell_connectivity_a.size] = cell_connectivity_a
424 merged_cell_connectivity[cell_connectivity_a.size :] = cell_connectivity_b
426 mask = _np.ones(cell_connectivity_b.size, dtype=bool)
427 mask[offsets_b[:-1]] = False
429 merged_cell_connectivity_view_part_b = merged_cell_connectivity[
430 cell_connectivity_a.size :
431 ]
432 merged_cell_connectivity_view_part_b[mask] += mesh_representation_a.n_points
434 # Merge the cell types.
435 merged_cell_types = _np.concatenate(
436 (
437 mesh_representation_a.cell_types,
438 mesh_representation_b.cell_types,
439 )
440 )
442 # Merge the contained data dictionaries.
443 def _merge_data_dicts(
444 dict_a: dict[str, _NDArray],
445 size_a: int,
446 dict_b: dict[str, _NDArray],
447 size_b: int,
448 ) -> dict[str, _NDArray]:
449 """Merge the given data dictionaries and fill in non-existing data."""
451 def _ensure_array_size(size: int, reference_array: _NDArray) -> _NDArray:
452 """Create an empty array that matches the columns of the reference array and
453 has the given size, i.e., number of rows."""
454 new_shape = reference_array.shape
455 if len(new_shape) == 1:
456 new_shape = (size,)
457 elif len(new_shape) == 2:
458 new_shape = (size, new_shape[1])
459 else:
460 raise ValueError(f"Got unexpected array shape {new_shape}.")
461 return _np.zeros(new_shape, dtype=reference_array.dtype)
463 all_keys = set(dict_a.keys()) | set(dict_b.keys())
464 merged_data = {}
465 for key in all_keys:
466 data_a = dict_a.get(key, None)
467 data_b = dict_b.get(key, None)
468 if data_a is None:
469 data_a = _ensure_array_size(size_a, data_b)
470 elif data_b is None:
471 data_b = _ensure_array_size(size_b, data_a)
472 merged_data[key] = _np.concatenate((data_a, data_b))
473 return merged_data
475 merged_cell_data = _merge_data_dicts(
476 mesh_representation_a.cell_data,
477 mesh_representation_a.n_cells,
478 mesh_representation_b.cell_data,
479 mesh_representation_b.n_cells,
480 )
481 merged_point_data = _merge_data_dicts(
482 mesh_representation_a.point_data,
483 mesh_representation_a.n_points,
484 mesh_representation_b.point_data,
485 mesh_representation_b.n_points,
486 )
488 return MeshRepresentation(
489 cell_connectivity=merged_cell_connectivity,
490 cell_types=merged_cell_types,
491 points=merged_points,
492 cell_data=merged_cell_data,
493 point_data=merged_point_data,
494 )