Coverage for src/beamme/four_c/dbc_monitor.py: 90%
89 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 function converts the DBC monitor log files to Neumann boundary
23conditions in a mesh."""
25import numpy as _np
26import yaml as _yaml
28from beamme.core.boundary_condition import BoundaryCondition as _BoundaryCondition
29from beamme.core.conf import bme as _bme
30from beamme.core.function import Function as _Function
31from beamme.core.geometry_set import GeometrySet as _GeometrySet
32from beamme.core.mesh import Mesh as _Mesh
33from beamme.four_c.function_utility import (
34 create_linear_interpolation_function as _create_linear_interpolation_function,
35)
36from beamme.four_c.function_utility import (
37 ensure_length_of_function_array as _ensure_length_of_function_array,
38)
41def linear_time_transformation(
42 time, values, time_span, *, flip=False, valid_start_and_end_point=False
43):
44 """Performs a transformation of the time to a new interval with an
45 appropriate value vector.
47 Args
48 ----
49 time: _np.array
50 array with time values
51 values: _np.array
52 corresponding values to time
53 time_span: [list] with 2 or 3 entries:
54 time_span[0:2] defines the time interval to which the initial time interval should be scaled.
55 time_span[3] optional timepoint to repeat last value
56 flip: Bool
57 Flag if the values should be reversed
58 valid_start_and_end_point: Bool
59 optionally adds a valid starting point at t=0 and timespan[3] if provided
60 """
62 # flip values if desired and adjust time
63 if flip is True:
64 values = _np.flipud(values)
65 time = _np.flip(-time) + time[-1]
67 # transform time to interval
68 min_t = _np.min(time)
69 max_t = _np.max(time)
71 # scaling/transforming the time into the user defined time
72 time = time_span[0] + (time - min_t) * (time_span[1] - time_span[0]) / (
73 max_t - min_t
74 )
76 # ensure that start time is valid
77 if valid_start_and_end_point and time[0] > 0.0:
78 # add starting time 0
79 time = _np.append(0.0, time)
81 # add first coordinate again at the beginning of the array
82 if len(values.shape) == 1:
83 values = _np.append(values[0], values)
84 else:
85 values = _np.append(values[0], values).reshape(
86 values.shape[0] + 1, values.shape[1]
87 )
89 # repeat last value at provided time point
90 if valid_start_and_end_point and len(time_span) > 2:
91 if time_span[2] > time_span[1]:
92 time = _np.append(time, time_span[2])
93 values = _np.append(values, values[-1]).reshape(
94 values.shape[0] + 1, values.shape[1]
95 )
96 if not valid_start_and_end_point and len(time_span) > 2:
97 raise Warning("You specified unnecessarily a third component of time_span.")
99 return time, values
102def read_dbc_monitor_file(file_path):
103 """Load the Dirichlet boundary condition monitor log and return the data as
104 well as the nodes of this boundary condition.
106 Args
107 ----
108 file_path: str
109 Path to the Dirichlet boundary condition monitor log.
111 Return
112 ----
113 [node_ids], [time], [force], [moment]
114 """
116 with open(file_path, "r") as f:
117 dbc_monitor_file = _yaml.safe_load(f)
119 nodes = dbc_monitor_file["dbc monitor condition"]["node gids"]
121 time = []
122 force = []
123 moment = []
124 for time_step_data in dbc_monitor_file["dbc monitor condition data"]:
125 time.append(time_step_data["time"])
126 force.append(time_step_data["f"])
127 moment.append(time_step_data["m"])
129 return (
130 nodes,
131 _np.asarray(time, dtype=float),
132 _np.asarray(force, dtype=float),
133 _np.asarray(moment, dtype=float),
134 )
137def add_point_neuman_condition_to_mesh(
138 mesh: _Mesh,
139 nodes: list[int],
140 function_array: list[_Function],
141 force: _np.ndarray,
142 *,
143 n_dof: int = 3,
144):
145 """Adds a Neumann boundary condition to a mesh for the given node_ids with
146 the function_array and force values by creating a new geometry set.
148 Args
149 ----
150 mesh: Mesh
151 Mesh where the boundary conditions are added to
152 nodes: [node_id]
153 list containing the ids of the nodes for the condition
154 function_array: [function]
155 list with functions
156 force: [_np.ndarray]
157 values to scale the function array with
158 n_dof: int
159 Number of DOFs per node.
160 """
162 # check if the dimensions of force and functions match
163 if force.size != 3:
164 raise ValueError(
165 f"The forces vector must have dimensions [3x1] not [{force.size}x1]"
166 )
168 function_array = _ensure_length_of_function_array(function_array, 3)
170 # Add the function to the mesh, if they are not previously added.
171 for function in function_array:
172 mesh.add(function)
174 # Create GeometrySet with nodes.
175 mesh_nodes = [mesh.nodes[i_node] for i_node in nodes]
176 geo = _GeometrySet(mesh_nodes)
178 # Create the Boundary Condition.
179 bc = _BoundaryCondition(
180 geo,
181 {
182 "NUMDOF": n_dof,
183 "ONOFF": [1, 1, 1] + [0] * (n_dof - 3),
184 "VAL": force.tolist() + [0] * (n_dof - 3),
185 "FUNCT": function_array + [0] * (n_dof - 3),
186 },
187 bc_type=_bme.bc.neumann,
188 )
189 mesh.add(bc)
192def dbc_monitor_to_mesh_all_values(
193 mesh: _Mesh,
194 file_path: str,
195 *,
196 steps: list[int] = [],
197 time_span: list[int] = [0, 1, 2],
198 type: str | None = "linear",
199 flip_time_values: bool = False,
200 functions: list[_Function] = [],
201 **kwargs,
202):
203 """Extracts all the force values of the monitored Dirichlet boundary
204 condition and converts them into a Function with a Neumann boundary
205 condition for a given mesh. The monitor log force values must be obtained
206 from a previous simulation with constant step size. The discretization of
207 the previous simulation must be identical to the one within the mesh. The
208 extracted force values are passed to a linear interpolation 4C-function. It
209 is advisable to only call this function once all nodes have been added to
210 the mesh.
212 Args
213 ----
214 mesh: Mesh
215 The mesh where the created Neumann boundary condition is added
216 to. The nodes(e.g., discretization) referred to in the log file must match with the ones
217 in the mesh.
218 file_path: str
219 Path to the Dirichlet boundary condition log file.
220 steps: [int,int]
221 Index range of which the force values are extracted. Default 0 and -1 extracts every point from the array.
222 time_span: [t1, t2, t3] in float
223 Transforms the given time array into this specific format.
224 The time array always starts at 0 and ends at t3 to ensure a valid simulation.
225 type: str or linear
226 two types are available:
227 1) "linear": not specified simply extract all values and apply them between time interval t1 and t2.
228 2) "hat": puts the values first until last value is reached and then decays them back to first value.
229 Interpolation starting from t1 going to the last value at (t1+t2)/2 and going back to the value at time t2.
230 flip_time_values: bool
231 indicates, if the extracted forces should be flipped or rearranged wrt. to the time
232 For flip_time_values=true, the forces at the final time are applied at t_start.
233 functions: [Function, Function, Function]
234 Array consisting of 3 custom functions(x,y,z). The value for boundary condition is selected from the last steps.
235 """
237 nodes, time, force, _ = read_dbc_monitor_file(file_path)
239 # The forces are the negative reactions at the Dirichlet boundaries.
240 force *= -1.0
242 # if special index range is provided use it
243 if steps:
244 time = time[steps[0] : steps[1] + 1]
245 force = force[steps[0] : steps[1] + 1, :]
246 else:
247 # otherwise define steps from start to end
248 steps = [0, -1]
250 # apply transformations to time and forces according to the schema
251 if type == "linear":
252 if not len(time_span) == 2:
253 raise ValueError(
254 f"Please provide a time_span with size 1x2 not {len(time_span)}"
255 )
257 time, force = linear_time_transformation(
258 time, force, time_span, flip=flip_time_values
259 )
260 if len(functions) != 3:
261 print("Please provide a list with three valid Functions.")
263 elif type == "hat":
264 if not len(time_span) == 3:
265 raise ValueError(
266 f"Please provide a time_span with size 1x3 not {len(time_span)}"
267 )
269 if functions:
270 print(
271 "You selected type",
272 type,
273 ", however the provided functions ",
274 functions,
275 " are overwritten.",
276 )
277 functions = []
279 # create the two intervals
280 time1, force1 = linear_time_transformation(
281 time, force, time_span[0:2], flip=flip_time_values
282 )
283 time2, force2 = linear_time_transformation(
284 time, force, time_span[1:3], flip=(not flip_time_values)
285 )
287 # remove first element since it is duplicated zero
288 _np.delete(time2, 0)
289 _np.delete(force2, 0)
291 # add the respective force
292 time = _np.concatenate((time1, time2[1:]))
293 force = _np.concatenate((force1, force2[1:]), axis=0)
295 else:
296 raise ValueError(
297 "The selected type: "
298 + str(type)
299 + " is currently not supported. Feel free to add it here."
300 )
302 # overwrite the function, if one is provided since for the specific types the function is generated
303 if not type == "linear":
304 for dim in range(force.shape[1]):
305 # create a linear function with the force values per dimension
306 fun = _create_linear_interpolation_function(
307 time, force[:, dim], function_type="SYMBOLIC_FUNCTION_OF_TIME"
308 )
310 # add the function to the mesh
311 mesh.add(fun)
313 # store function
314 functions.append(fun)
316 # now set forces to 1 since the force values are already extracted in the function's values
317 force = _np.zeros_like(force) + 1.0
319 elif len(functions) != 3:
320 raise ValueError("Please provide functions with ")
322 # Create condition in mesh
323 add_point_neuman_condition_to_mesh(
324 mesh, nodes, functions, force[steps[1]], **kwargs
325 )
328def dbc_monitor_to_mesh(
329 mesh: _Mesh,
330 file_path: str,
331 *,
332 step: int = -1,
333 function: _Function,
334 **kwargs,
335):
336 """Converts the last value of a Dirichlet boundary condition monitor log to
337 a Neumann boundary condition in the mesh.
339 Args
340 ----
341 mesh: Mesh
342 The mesh where the created Neumann boundary condition is added to.
343 The nodes referred to in the log file have to match with the ones
344 in the mesh. It is advisable to only call this function once
345 all nodes have been added to the mesh.
346 file_path: str
347 Path to the Dirichlet boundary condition log file.
348 step: int
349 Step values to be used. Default is -1, i.e. the last step.
350 function: Function
351 Function for the Neumann boundary condition.
352 """
354 # read the force
355 nodes, _, force, _ = read_dbc_monitor_file(file_path)
357 # The forces are the negative reactions at the Dirichlet boundaries.
358 force *= -1.0
360 # Create condition in mesh
361 add_point_neuman_condition_to_mesh(
362 mesh, nodes, [function] * 3, force[step], **kwargs
363 )