Coverage for src/beamme/core/mesh_utils.py: 92%
63 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 utility functions for meshes."""
24import numpy as _np
26from beamme.core.conf import bme as _bme
27from beamme.core.coupling import Coupling as _Coupling
28from beamme.core.geometry_set import GeometrySetBase as _GeometrySetBase
29from beamme.core.mesh import Mesh as _Mesh
30from beamme.core.mesh_representation import MeshRepresentation as _MeshRepresentation
31from beamme.core.mesh_representation import (
32 string_to_geometry_set_info as _string_to_geometry_set_info,
33)
34from beamme.core.node import Node as _Node
37def get_coupled_nodes_to_master_map(
38 mesh: _Mesh, *, assign_i_global: bool = False
39) -> tuple[dict[_Node, _Node], list[_Node]]:
40 """Get a mapping of nodes in a mesh that should be "replaced" because they are
41 coupled via a joint.
43 In some finite element (FE) solvers, nodes coupled via joints are resolved
44 by assigning a "master" node to represent the joint. This function identifies
45 such nodes and creates a mapping where each coupled node is mapped to its
46 master node.
48 Args
49 ----
50 mesh:
51 Input mesh
52 assign_i_global:
53 If this flag is set, the global indices are set in the node objects.
55 Return
56 ----
57 replaced_node_to_master_map:
58 A dictionary mapping each "replaced" node to its "master" node.
59 unique_nodes:
60 A list containing all unique nodes in the mesh, i.e., all nodes which
61 are not coupled and the master nodes.
62 """
63 # Get a dictionary that maps the "replaced" nodes to the "master" ones
64 replaced_node_to_master_map = {}
65 for coupling in mesh.boundary_conditions[_bme.bc.point_coupling, _bme.geo.point]:
66 if coupling.data is not _bme.coupling_dof.fix:
67 raise ValueError(
68 "This function is only implemented for rigid joints at the DOFs"
69 )
70 coupling_nodes = coupling.geometry_set.get_points()
71 for node in coupling_nodes[1:]:
72 replaced_node_to_master_map[node] = coupling_nodes[0]
74 # Check that no "replaced" node is a "master" node
75 master_nodes = set(replaced_node_to_master_map.values())
76 for replaced_node in replaced_node_to_master_map.keys():
77 if replaced_node in master_nodes:
78 raise ValueError(
79 "A replaced node is also a master nodes. This is not supported"
80 )
82 # Get all unique nodes
83 unique_nodes = [
84 node for node in mesh.nodes if node not in replaced_node_to_master_map
85 ]
87 # Optionally number the nodes
88 if assign_i_global:
89 for i_node, node in enumerate(unique_nodes):
90 node.i_global = i_node
91 for replaced_node, master_node in replaced_node_to_master_map.items():
92 replaced_node.i_global = master_node.i_global
94 # Return the mapping
95 return replaced_node_to_master_map, unique_nodes
98def apply_nodal_coupling_to_mesh_representation(
99 mesh_representation: _MeshRepresentation,
100 geometry_sets_to_i_global: dict[_GeometrySetBase, int],
101 coupling_conditions: list[_Coupling],
102):
103 """Modify a mesh representation such that coupled nodes are represented by a single
104 node.
106 Args:
107 mesh_representation: The mesh representation where coupling nodes should
108 be merged. Will be modified in-place.
109 geometry_sets_to_i_global: A mapping from geometry sets to their global
110 IDs.
111 coupling_conditions: A list of coupling conditions that define which
112 nodes should be coupled.
113 """
114 # Get a dictionary that maps the "replaced" nodes to the "master" ones
115 replaced_node_to_master_map = {}
116 for coupling in coupling_conditions:
117 if not coupling.data == _bme.coupling_dof.fix:
118 raise ValueError(
119 "This function is only implemented for rigid joints at the DOFs"
120 )
121 geometry_set_i_global = geometry_sets_to_i_global[coupling.geometry_set]
122 for name in mesh_representation.point_data.keys():
123 info = _string_to_geometry_set_info(name)
124 if info is not None and info.i_global == geometry_set_i_global:
125 coupled_node_ids = mesh_representation.point_data[name].nonzero()[0]
126 for node in coupled_node_ids[1:]:
127 replaced_node_to_master_map[node] = coupled_node_ids[0]
128 break
129 else:
130 raise ValueError(
131 f"Could not find geometry set with i_global {geometry_set_i_global} in the mesh representation"
132 )
134 # Check that no "replaced" node is a "master" node
135 master_nodes = set(replaced_node_to_master_map.values())
136 replaced_nodes = set(replaced_node_to_master_map.keys())
137 if len(master_nodes & replaced_nodes) > 0:
138 raise ValueError(
139 "A replaced node is also a master nodes. This is not supported. "
140 f"Replaced nodes: {replaced_nodes}, master nodes: {master_nodes}, "
141 f"intersection: {master_nodes & replaced_nodes}"
142 )
144 # Get mask that filters out replaced nodes
145 point_mask = _np.ones(mesh_representation.n_points, dtype=bool)
146 point_mask[list(replaced_nodes)] = False
148 # Get the node mapping vector from the old IDs to the new ones
149 mapping_vector = _np.zeros(mesh_representation.n_points, dtype=int)
150 mapping_vector[point_mask] = _np.arange(
151 mesh_representation.n_points - len(replaced_nodes)
152 )
153 for replaced_node, master_node in replaced_node_to_master_map.items():
154 # Map the replaced node to the (new) ID of the master node
155 mapping_vector[replaced_node] = mapping_vector[master_node]
157 # Merge geometry-set membership of replaced nodes into their master nodes before
158 # removing points (otherwise geometry sets / BC references can be lost).
159 for name, values in mesh_representation.point_data.items():
160 if _string_to_geometry_set_info(name) is not None:
161 for replaced_node, master_node in replaced_node_to_master_map.items():
162 values[master_node] = max(values[master_node], values[replaced_node])
164 # Filter the point data vectors in the mesh representation.
165 # For now, we simply drop the "replaced" points.
166 mesh_representation.points = mesh_representation.points[point_mask]
167 for name in mesh_representation.point_data.keys():
168 mesh_representation.point_data[name] = mesh_representation.point_data[name][
169 point_mask
170 ]
172 # Adapt the connectivity such that the "replaced" nodes are mapped to the "master" ones.
173 connectivity_mask = _np.ones(mesh_representation.cell_connectivity.size, dtype=bool)
174 connectivity_mask[mesh_representation.cell_connectivity_offsets[:-1]] = False
175 point_indices_old = mesh_representation.cell_connectivity[connectivity_mask]
176 point_indices_new = mapping_vector[point_indices_old]
177 mesh_representation.cell_connectivity[connectivity_mask] = point_indices_new