Coverage for src/beamme/four_c/dbc_monitor.py: 90%
89 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 function converts the DBC monitor log files to Neumann boundary conditions in a
23mesh."""
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 appropriate value
45 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 """
61 # flip values if desired and adjust time
62 if flip is True:
63 values = _np.flipud(values)
64 time = _np.flip(-time) + time[-1]
66 # transform time to interval
67 min_t = _np.min(time)
68 max_t = _np.max(time)
70 # scaling/transforming the time into the user defined time
71 time = time_span[0] + (time - min_t) * (time_span[1] - time_span[0]) / (
72 max_t - min_t
73 )
75 # ensure that start time is valid
76 if valid_start_and_end_point and time[0] > 0.0:
77 # add starting time 0
78 time = _np.append(0.0, time)
80 # add first coordinate again at the beginning of the array
81 if len(values.shape) == 1:
82 values = _np.append(values[0], values)
83 else:
84 values = _np.append(values[0], values).reshape(
85 values.shape[0] + 1, values.shape[1]
86 )
88 # repeat last value at provided time point
89 if valid_start_and_end_point and len(time_span) > 2:
90 if time_span[2] > time_span[1]:
91 time = _np.append(time, time_span[2])
92 values = _np.append(values, values[-1]).reshape(
93 values.shape[0] + 1, values.shape[1]
94 )
95 if not valid_start_and_end_point and len(time_span) > 2:
96 raise Warning("You specified unnecessarily a third component of time_span.")
98 return time, values
101def read_dbc_monitor_file(file_path):
102 """Load the Dirichlet boundary condition monitor log and return the data as well as
103 the nodes of this boundary condition.
105 Args
106 ----
107 file_path: str
108 Path to the Dirichlet boundary condition monitor log.
110 Return
111 ----
112 [node_ids], [time], [force], [moment]
113 """
114 with open(file_path, "r") as f:
115 dbc_monitor_file = _yaml.safe_load(f)
117 nodes = dbc_monitor_file["dbc monitor condition"]["node gids"]
119 time = []
120 force = []
121 moment = []
122 for time_step_data in dbc_monitor_file["dbc monitor condition data"]:
123 time.append(time_step_data["time"])
124 force.append(time_step_data["f"])
125 moment.append(time_step_data["m"])
127 return (
128 nodes,
129 _np.asarray(time, dtype=float),
130 _np.asarray(force, dtype=float),
131 _np.asarray(moment, dtype=float),
132 )
135def add_point_neuman_condition_to_mesh(
136 mesh: _Mesh,
137 nodes: list[int],
138 function_array: list[_Function],
139 force: _np.ndarray,
140 *,
141 n_dof: int = 3,
142):
143 """Adds a Neumann boundary condition to a mesh for the given node_ids with the
144 function_array and force values by creating a new geometry set.
146 Args
147 ----
148 mesh: Mesh
149 Mesh where the boundary conditions are added to
150 nodes: [node_id]
151 list containing the ids of the nodes for the condition
152 function_array: [function]
153 list with functions
154 force: [_np.ndarray]
155 values to scale the function array with
156 n_dof: int
157 Number of DOFs per node.
158 """
159 # check if the dimensions of force and functions match
160 if force.size != 3:
161 raise ValueError(
162 f"The forces vector must have dimensions [3x1] not [{force.size}x1]"
163 )
165 function_array = _ensure_length_of_function_array(function_array, 3)
167 # Add the function to the mesh, if they are not previously added.
168 for function in function_array:
169 mesh.add(function)
171 # Create GeometrySet with nodes.
172 mesh_nodes = [mesh.nodes[i_node] for i_node in nodes]
173 geo = _GeometrySet(mesh_nodes)
175 # Create the Boundary Condition.
176 bc = _BoundaryCondition(
177 geo,
178 {
179 "NUMDOF": n_dof,
180 "ONOFF": [1, 1, 1] + [0] * (n_dof - 3),
181 "VAL": force.tolist() + [0] * (n_dof - 3),
182 "FUNCT": function_array + [0] * (n_dof - 3),
183 },
184 bc_type=_bme.bc.neumann,
185 )
186 mesh.add(bc)
189def dbc_monitor_to_mesh_all_values(
190 mesh: _Mesh,
191 file_path: str,
192 *,
193 steps: list[int] = [],
194 time_span: list[int] = [0, 1, 2],
195 type: str | None = "linear",
196 flip_time_values: bool = False,
197 functions: list[_Function] = [],
198 **kwargs,
199):
200 """Extracts all the force values of the monitored Dirichlet boundary condition and
201 converts them into a Function with a Neumann boundary condition for a given mesh.
202 The monitor log force values must be obtained from a previous simulation with
203 constant step size. The discretization of the previous simulation must be identical
204 to the one within the mesh. The extracted force values are passed to a linear
205 interpolation 4C-function. It is advisable to only call this function once all nodes
206 have been added to the mesh.
208 Args
209 ----
210 mesh: Mesh
211 The mesh where the created Neumann boundary condition is added
212 to. The nodes(e.g., discretization) referred to in the log file must match with the ones
213 in the mesh.
214 file_path: str
215 Path to the Dirichlet boundary condition log file.
216 steps: [int,int]
217 Index range of which the force values are extracted. Default 0 and -1 extracts every point from the array.
218 time_span: [t1, t2, t3] in float
219 Transforms the given time array into this specific format.
220 The time array always starts at 0 and ends at t3 to ensure a valid simulation.
221 type: str or linear
222 two types are available:
223 1) "linear": not specified simply extract all values and apply them between time interval t1 and t2.
224 2) "hat": puts the values first until last value is reached and then decays them back to first value.
225 Interpolation starting from t1 going to the last value at (t1+t2)/2 and going back to the value at time t2.
226 flip_time_values: bool
227 indicates, if the extracted forces should be flipped or rearranged wrt. to the time
228 For flip_time_values=true, the forces at the final time are applied at t_start.
229 functions: [Function, Function, Function]
230 Array consisting of 3 custom functions(x,y,z). The value for boundary condition is selected from the last steps.
231 """
232 nodes, time, force, _ = read_dbc_monitor_file(file_path)
234 # The forces are the negative reactions at the Dirichlet boundaries.
235 force *= -1.0
237 # if special index range is provided use it
238 if steps:
239 time = time[steps[0] : steps[1] + 1]
240 force = force[steps[0] : steps[1] + 1, :]
241 else:
242 # otherwise define steps from start to end
243 steps = [0, -1]
245 # apply transformations to time and forces according to the schema
246 if type == "linear":
247 if not len(time_span) == 2:
248 raise ValueError(
249 f"Please provide a time_span with size 1x2 not {len(time_span)}"
250 )
252 time, force = linear_time_transformation(
253 time, force, time_span, flip=flip_time_values
254 )
255 if len(functions) != 3:
256 print("Please provide a list with three valid Functions.")
258 elif type == "hat":
259 if not len(time_span) == 3:
260 raise ValueError(
261 f"Please provide a time_span with size 1x3 not {len(time_span)}"
262 )
264 if functions:
265 print(
266 "You selected type",
267 type,
268 ", however the provided functions ",
269 functions,
270 " are overwritten.",
271 )
272 functions = []
274 # create the two intervals
275 time1, force1 = linear_time_transformation(
276 time, force, time_span[0:2], flip=flip_time_values
277 )
278 time2, force2 = linear_time_transformation(
279 time, force, time_span[1:3], flip=(not flip_time_values)
280 )
282 # remove first element since it is duplicated zero
283 _np.delete(time2, 0)
284 _np.delete(force2, 0)
286 # add the respective force
287 time = _np.concatenate((time1, time2[1:]))
288 force = _np.concatenate((force1, force2[1:]), axis=0)
290 else:
291 raise ValueError(
292 "The selected type: "
293 + str(type)
294 + " is currently not supported. Feel free to add it here."
295 )
297 # overwrite the function, if one is provided since for the specific types the function is generated
298 if not type == "linear":
299 for dim in range(force.shape[1]):
300 # create a linear function with the force values per dimension
301 fun = _create_linear_interpolation_function(
302 time, force[:, dim], function_type="SYMBOLIC_FUNCTION_OF_TIME"
303 )
305 # add the function to the mesh
306 mesh.add(fun)
308 # store function
309 functions.append(fun)
311 # now set forces to 1 since the force values are already extracted in the function's values
312 force = _np.zeros_like(force) + 1.0
314 elif len(functions) != 3:
315 raise ValueError("Please provide functions with ")
317 # Create condition in mesh
318 add_point_neuman_condition_to_mesh(
319 mesh, nodes, functions, force[steps[1]], **kwargs
320 )
323def dbc_monitor_to_mesh(
324 mesh: _Mesh,
325 file_path: str,
326 *,
327 step: int = -1,
328 function: _Function,
329 **kwargs,
330):
331 """Converts the last value of a Dirichlet boundary condition monitor log to a
332 Neumann boundary condition in the mesh.
334 Args
335 ----
336 mesh: Mesh
337 The mesh where the created Neumann boundary condition is added to.
338 The nodes referred to in the log file have to match with the ones
339 in the mesh. It is advisable to only call this function once
340 all nodes have been added to the mesh.
341 file_path: str
342 Path to the Dirichlet boundary condition log file.
343 step: int
344 Step values to be used. Default is -1, i.e. the last step.
345 function: Function
346 Function for the Neumann boundary condition.
347 """
348 # read the force
349 nodes, _, force, _ = read_dbc_monitor_file(file_path)
351 # The forces are the negative reactions at the Dirichlet boundaries.
352 force *= -1.0
354 # Create condition in mesh
355 add_point_neuman_condition_to_mesh(
356 mesh, nodes, [function] * 3, force[step], **kwargs
357 )