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