Coverage for src/beamme/core/mesh_representation.py: 93%
152 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-08 11:03 +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 dataclasses import dataclass as _dataclass
25from itertools import repeat as _repeat
26from typing import Any as _Any
27from typing import Iterable as _Iterable
29import numpy as _np
30from numpy.typing import NDArray as _NDArray
32from beamme.core.conf import Geometry as _Geometry
33from beamme.core.conf import bme as _bme
35MESH_REPRESENTATION_MAPPINGS: dict[str, _Any] = {}
36# fmt: off
37MESH_REPRESENTATION_MAPPINGS[
38 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
39] = {
41 # Only list the non-standard mappings
42 (_bme.element_type.solid, 20):
43 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15],
44 (_bme.element_type.solid, 27):
45 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15,
46 24, 22, 21, 23, 20, 25, 26],
47}
48# fmt: on
49MESH_REPRESENTATION_MAPPINGS[
50 "element_type_and_n_nodes_to_connectivity_mapping_vtk_to_beamme"
51] = {
52 # Only list the non-standard mappings
53 (_bme.element_type.solid, 20): _np.argsort(
54 MESH_REPRESENTATION_MAPPINGS[
55 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
56 ][(_bme.element_type.solid, 20)]
57 ),
58 (_bme.element_type.solid, 27): _np.argsort(
59 MESH_REPRESENTATION_MAPPINGS[
60 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
61 ][(_bme.element_type.solid, 27)]
62 ),
63}
65GEOMETRY_SET_INFO_PREFIX = "geometry_set"
68@_dataclass
69class GeometrySetInfo:
70 """Data class that contains information of a geometry set."""
72 geometry_type: _Geometry
73 i_global: int
74 name: str | None = None
75 point_flag_vector: _NDArray | None = None
76 cell_flag_vector: _NDArray | None = None
78 def __str__(self):
79 """Return a string representation of this geometry set info - the node and cell flags are not included."""
80 name = self.name
81 return (
82 f"{GEOMETRY_SET_INFO_PREFIX}_{self.i_global}_{self.geometry_type.name}"
83 + (f"_{name}" if name is not None else "")
84 )
87def string_to_geometry_set_info(name: str) -> GeometrySetInfo | None:
88 """Extract the geometry set information from a given string."""
89 if not name.startswith(GEOMETRY_SET_INFO_PREFIX):
90 return None
92 name_without_prefix = name[len(GEOMETRY_SET_INFO_PREFIX) + 1 :]
93 split = name_without_prefix.split("_", 2)
94 i_global = int(split[0])
95 geometry_type = _Geometry[split[1]]
96 return GeometrySetInfo(
97 i_global=i_global,
98 geometry_type=geometry_type,
99 name=split[2] if len(split) == 3 else None,
100 )
103class MeshRepresentation:
104 """Class representing a generic mesh."""
106 def __init__(
107 self,
108 cell_connectivity: _NDArray[_np.integer] | None = None,
109 cell_types: _NDArray[_np.integer] | None = None,
110 points: _NDArray[_np.floating] | None = None,
111 geometry_sets: list[GeometrySetInfo] | None = None,
112 cell_data: dict[str, _NDArray | None] | None = None,
113 point_data: dict[str, _NDArray | None] | None = None,
114 ):
115 def _convert_argument_numpy(
116 argument: _NDArray | None,
117 default_shape: tuple[int, ...],
118 dtype: type,
119 ) -> _NDArray:
120 """Convert given array arguments so we can store them in this
121 object."""
122 if argument is None:
123 return _np.empty(default_shape, dtype=dtype)
124 else:
125 return _np.asarray(argument, dtype=dtype)
127 def _filter_none_entries(
128 argument: dict[str, _NDArray | None] | None, expected_size: int
129 ) -> dict[str, _NDArray]:
130 """Check if a dictionary is given, and if so, filter None entries
131 from it."""
132 if argument is None:
133 return {}
134 else:
135 data = {
136 key: _np.asarray(value)
137 for key, value in argument.items()
138 if value is not None
139 }
140 for key, value in data.items():
141 if len(value) != expected_size:
142 raise ValueError(
143 f"Data field {key} has size {len(value)}, but expected size is {expected_size}."
144 )
145 return data
147 self.cell_connectivity = _convert_argument_numpy(
148 cell_connectivity, default_shape=(0,), dtype=int
149 )
150 self.cell_types = _convert_argument_numpy(
151 cell_types, default_shape=(0,), dtype=int
152 )
153 self.points = _convert_argument_numpy(points, default_shape=(0, 3), dtype=float)
155 self.cell_data = _filter_none_entries(cell_data, self.n_cells)
156 self.point_data = _filter_none_entries(point_data, self.n_points)
158 def _store_geometry_set_data(
159 data_field: str,
160 name: str,
161 field_vector: _NDArray | None,
162 expected_size: int,
163 ) -> None:
164 """Store point or cell data defining a geometry set."""
165 if field_vector is not None:
166 if len(field_vector) != expected_size:
167 raise ValueError(
168 f"Geometry set {name} has size {len(field_vector)},"
169 f" but expected size is {expected_size}."
170 )
171 getattr(self, data_field)[str(name)] = field_vector
173 if geometry_sets is not None:
174 for geometry_set in geometry_sets:
175 geometry_set_name = str(geometry_set)
176 _store_geometry_set_data(
177 "point_data",
178 geometry_set_name,
179 geometry_set.point_flag_vector,
180 self.n_points,
181 )
182 _store_geometry_set_data(
183 "cell_data",
184 geometry_set_name,
185 geometry_set.cell_flag_vector,
186 self.n_cells,
187 )
189 # Get the offset array so we can iterate over the connectivity in a
190 # performant manner.
191 self.cell_connectivity_offsets = _np.zeros(self.n_cells + 1, dtype=int)
192 for i_cell in range(self.n_cells):
193 last_offset = self.cell_connectivity_offsets[i_cell]
194 if last_offset >= len(self.cell_connectivity):
195 raise ValueError(
196 "Invalid offset, is larger than the size of the connectivity."
197 )
198 cell_size = self.cell_connectivity[last_offset]
199 if cell_size < 1:
200 raise ValueError("Invalid offset, has to be at least 1.")
201 self.cell_connectivity_offsets[i_cell + 1] = (
202 cell_size + self.cell_connectivity_offsets[i_cell] + 1
203 )
204 if self.cell_connectivity_offsets[-1] != len(self.cell_connectivity):
205 raise ValueError("Invalid cell connectivity offsets.")
207 @property
208 def n_cells(self) -> int:
209 """Return the number of cells in this mesh representation.
211 Returns:
212 Number of cells in this mesh representation.
213 """
214 return len(self.cell_types)
216 @property
217 def n_points(self) -> int:
218 """Return the number of points in this mesh representation.
220 Returns:
221 Number of points in this mesh representation.
222 """
223 return len(self.points)
225 def connectivity_iterator(self) -> _Iterable:
226 """Return an iterator over the cell connectivity.
228 Returns:
229 An iterator that returns the connectivity array for each cell.
230 """
231 for i in range(self.n_cells):
232 start = self.cell_connectivity_offsets[i] + 1
233 end = self.cell_connectivity_offsets[i + 1]
234 yield self.cell_connectivity[start:end]
236 def data_iterator(self, data_field: str, data_name: str) -> _Iterable:
237 """This method returns an iterator for the given data field and data
238 name.
240 This is useful, when looping over the data, as accessing the data field in each
241 loop iteration can be expensive. If the data field is not present, a iterator
242 is returned that will always return `None`.
244 Args:
245 data_field: The data field to get the iterator for. This can be either
246 "point_data" or "cell_data".
247 data_name: The name of the data to get the iterator for.
249 Returns:
250 An iterator to iterate through the data. If the data field is not present,
251 a iterator is returned that will always return `None`.
252 """
253 data_dict = getattr(self, data_field)
254 if data_name in data_dict:
255 return data_dict[data_name]
256 else:
257 match data_field:
258 case "point_data":
259 size = self.n_points
260 case "cell_data":
261 size = self.n_cells
262 case _:
263 raise ValueError(f"Invalid data field: {data_field}")
264 return _repeat(None, size)
266 def offset_indices(
267 self,
268 element_type_id_offset: int | None = None,
269 material_offset: int | None = None,
270 geometry_set_offset: int | None = None,
271 ) -> None:
272 """Add offsets to the internal indices of this mesh representation.
274 This is required when merging multiple mesh representations together,
275 to ensure that there are no index conflicts.
277 Args:
278 element_type_id_offset: The offset to add to the element type IDs.
279 material_offset: The offset to add to the material IDs.
280 geometry_set_offset: The offset to add to the geometry set IDs.
281 """
283 if element_type_id_offset is not None:
284 if "element_type_id" in self.cell_data:
285 self.cell_data["element_type_id"] += element_type_id_offset
287 if material_offset is not None:
288 if "material_id" in self.cell_data:
289 # Add the offset to all material entries that are not -1. We use
290 # -1 to indicate no assigned material and those values should not change.
291 self.cell_data["material_id"][self.cell_data["material_id"] >= 0] += (
292 material_offset
293 )
295 if geometry_set_offset is not None:
296 for field_type in ("point_data", "cell_data"):
297 new_dict = {}
298 data_field = getattr(self, field_type)
299 # We change the keys in the loop, so we iterate over a copy of them.
300 for name in list(data_field.keys()):
301 info = string_to_geometry_set_info(name)
302 if info is not None:
303 new_name = str(
304 GeometrySetInfo(
305 geometry_type=info.geometry_type,
306 i_global=info.i_global + geometry_set_offset,
307 name=info.name,
308 )
309 )
310 new_dict[new_name] = data_field.pop(name)
311 # Add the renamed geometry sets back to the data field.
312 setattr(self, field_type, {**data_field, **new_dict})
315def merge_mesh_representations(
316 mesh_representation_a: MeshRepresentation, mesh_representation_b: MeshRepresentation
317) -> MeshRepresentation:
318 """Merge two mesh representations.
320 Args:
321 mesh_representation_a: First mesh representation.
322 mesh_representation_b: Second mesh representation.
324 Returns:
325 A merged mesh representation. Depending on the input data, this might
326 be a reference to one of the input mesh representations.
327 """
329 def _is_empty(mesh_representation: MeshRepresentation) -> bool:
330 """Check if the mesh representation contains any points or cells."""
331 return mesh_representation.n_points == 0 and mesh_representation.n_cells == 0
333 # If one of the two mesh representations is empty, we can simply return the other (which
334 # could be empty as well).
335 if _is_empty(mesh_representation_a):
336 return mesh_representation_b
337 elif _is_empty(mesh_representation_b):
338 return mesh_representation_a
340 # Merge the points
341 merged_points = _np.vstack(
342 (mesh_representation_a.points, mesh_representation_b.points)
343 )
345 # Merge the connectivity
346 # Before the connectivity can be merged, the IDs of the second mesh have
347 # to be increased by the size of the first mesh. To do so, we create a
348 # mask to identify the entries that have to be incremented.
349 cell_connectivity_a = mesh_representation_a.cell_connectivity
350 cell_connectivity_b = mesh_representation_b.cell_connectivity
351 offsets_b = mesh_representation_b.cell_connectivity_offsets
353 merged_cell_connectivity = _np.empty(
354 cell_connectivity_a.size + cell_connectivity_b.size,
355 dtype=cell_connectivity_a.dtype,
356 )
357 merged_cell_connectivity[: cell_connectivity_a.size] = cell_connectivity_a
358 merged_cell_connectivity[cell_connectivity_a.size :] = cell_connectivity_b
360 mask = _np.ones(cell_connectivity_b.size, dtype=bool)
361 mask[offsets_b[:-1]] = False
363 merged_cell_connectivity_view_part_b = merged_cell_connectivity[
364 cell_connectivity_a.size :
365 ]
366 merged_cell_connectivity_view_part_b[mask] += mesh_representation_a.n_points
368 # Merge the cell types.
369 merged_cell_types = _np.concatenate(
370 (
371 mesh_representation_a.cell_types,
372 mesh_representation_b.cell_types,
373 )
374 )
376 # Merge the contained data dictionaries.
377 def _merge_data_dicts(
378 dict_a: dict[str, _NDArray],
379 size_a: int,
380 dict_b: dict[str, _NDArray],
381 size_b: int,
382 ) -> dict[str, _NDArray]:
383 """Merge the given data dictionaries and fill in non-existing data."""
385 def _ensure_array_size(size: int, reference_array: _NDArray) -> _NDArray:
386 """Create an empty array that matches the columns of the reference
387 array and has the given size, i.e., number of rows."""
388 new_shape = reference_array.shape
389 if len(new_shape) == 1:
390 new_shape = (size,)
391 elif len(new_shape) == 2:
392 new_shape = (size, new_shape[1])
393 else:
394 raise ValueError(f"Got unexpected array shape {new_shape}.")
395 return _np.zeros(new_shape, dtype=reference_array.dtype)
397 all_keys = set(dict_a.keys()) | set(dict_b.keys())
398 merged_data = {}
399 for key in all_keys:
400 data_a = dict_a.get(key, None)
401 data_b = dict_b.get(key, None)
402 if data_a is None:
403 data_a = _ensure_array_size(size_a, data_b)
404 elif data_b is None:
405 data_b = _ensure_array_size(size_b, data_a)
406 merged_data[key] = _np.concatenate((data_a, data_b))
407 return merged_data
409 merged_cell_data = _merge_data_dicts(
410 mesh_representation_a.cell_data,
411 mesh_representation_a.n_cells,
412 mesh_representation_b.cell_data,
413 mesh_representation_b.n_cells,
414 )
415 merged_point_data = _merge_data_dicts(
416 mesh_representation_a.point_data,
417 mesh_representation_a.n_points,
418 mesh_representation_b.point_data,
419 mesh_representation_b.n_points,
420 )
422 return MeshRepresentation(
423 cell_connectivity=merged_cell_connectivity,
424 cell_types=merged_cell_types,
425 points=merged_points,
426 cell_data=merged_cell_data,
427 point_data=merged_point_data,
428 )