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