Coverage for src/meshpy/four_c/dbc_monitor.py: 89%
99 statements
« prev ^ index » next coverage.py v7.9.0, created at 2025-06-13 04:26 +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 boundary
23conditions in a mesh."""
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.function import Function as _Function
32from meshpy.core.geometry_set import GeometrySet as _GeometrySet
33from meshpy.core.mesh import Mesh as _Mesh
34from meshpy.four_c.function_utility import (
35 create_linear_interpolation_function as _create_linear_interpolation_function,
36)
37from meshpy.four_c.function_utility import (
38 ensure_length_of_function_array as _ensure_length_of_function_array,
39)
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_mesh(
151 mesh: _Mesh,
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 a mesh for the given node_ids with
159 the function_array and force values by creating a new geometry set.
161 Args
162 ----
163 mesh: Mesh
164 Mesh where the boundary conditions are added to
165 nodes: [node_id]
166 list containing the ids of the nodes for the condition
167 function_array: [function]
168 list with functions
169 force: [_np.ndarray]
170 values to scale the function array with
171 n_dof: int
172 Number of DOFs per node.
173 """
175 # check if the dimensions of force and functions match
176 if force.size != 3:
177 raise ValueError(
178 f"The forces vector must have dimensions [3x1] not [{force.size}x1]"
179 )
181 function_array = _ensure_length_of_function_array(function_array, 3)
183 # Add the function to the mesh, if they are not previously added.
184 for function in function_array:
185 mesh.add(function)
187 # Create GeometrySet with nodes.
188 mesh_nodes = [mesh.nodes[i_node] for i_node in nodes]
189 geo = _GeometrySet(mesh_nodes)
191 # Create the Boundary Condition.
192 bc = _BoundaryCondition(
193 geo,
194 {
195 "NUMDOF": n_dof,
196 "ONOFF": [1, 1, 1] + [0] * (n_dof - 3),
197 "VAL": force.tolist() + [0] * (n_dof - 3),
198 "FUNCT": function_array + [0] * (n_dof - 3),
199 },
200 bc_type=_mpy.bc.neumann,
201 )
202 mesh.add(bc)
205def dbc_monitor_to_mesh_all_values(
206 mesh: _Mesh,
207 file_path: str,
208 *,
209 steps: list[int] = [],
210 time_span: list[int] = [0, 1, 2],
211 type: _Optional[str] = "linear",
212 flip_time_values: bool = False,
213 functions: list[_Function] = [],
214 **kwargs,
215):
216 """Extracts all the force values of the monitored Dirichlet boundary
217 condition and converts them into a Function with a Neumann boundary
218 condition for a given mesh. The monitor log force values must be obtained
219 from a previous simulation with constant step size. The discretization of
220 the previous simulation must be identical to the one within the mesh. The
221 extracted force values are passed to a linear interpolation 4C-function. It
222 is advisable to only call this function once all nodes have been added to
223 the mesh.
225 Args
226 ----
227 mesh: Mesh
228 The mesh where the created Neumann boundary condition is added
229 to. The nodes(e.g., discretization) referred to in the log file must match with the ones
230 in the mesh.
231 file_path: str
232 Path to the Dirichlet boundary condition log file.
233 steps: [int,int]
234 Index range of which the force values are extracted. Default 0 and -1 extracts every point from the array.
235 time_span: [t1, t2, t3] in float
236 Transforms the given time array into this specific format.
237 The time array always starts at 0 and ends at t3 to ensure a valid simulation.
238 type: str or linear
239 two types are available:
240 1) "linear": not specified simply extract all values and apply them between time interval t1 and t2.
241 2) "hat": puts the values first until last value is reached and then decays them back to first value.
242 Interpolation starting from t1 going to the last value at (t1+t2)/2 and going back to the value at time t2.
243 flip_time_values: bool
244 indicates, if the extracted forces should be flipped or rearranged wrt. to the time
245 For flip_time_values=true, the forces at the final time are applied at t_start.
246 functions: [Function, Function, Function]
247 Array consisting of 3 custom functions(x,y,z). The value for boundary condition is selected from the last steps.
248 """
250 nodes, time, force, _ = read_dbc_monitor_file(file_path)
252 # The forces are the negative reactions at the Dirichlet boundaries.
253 force *= -1.0
255 # if special index range is provided use it
256 if steps:
257 time = time[steps[0] : steps[1] + 1]
258 force = force[steps[0] : steps[1] + 1, :]
259 else:
260 # otherwise define steps from start to end
261 steps = [0, -1]
263 # apply transformations to time and forces according to the schema
264 if type == "linear":
265 if not len(time_span) == 2:
266 raise ValueError(
267 f"Please provide a time_span with size 1x2 not {len(time_span)}"
268 )
270 time, force = linear_time_transformation(
271 time, force, time_span, flip=flip_time_values
272 )
273 if len(functions) != 3:
274 print("Please provide a list with three valid Functions.")
276 elif type == "hat":
277 if not len(time_span) == 3:
278 raise ValueError(
279 f"Please provide a time_span with size 1x3 not {len(time_span)}"
280 )
282 if functions:
283 print(
284 "You selected type",
285 type,
286 ", however the provided functions ",
287 functions,
288 " are overwritten.",
289 )
290 functions = []
292 # create the two intervals
293 time1, force1 = linear_time_transformation(
294 time, force, time_span[0:2], flip=flip_time_values
295 )
296 time2, force2 = linear_time_transformation(
297 time, force, time_span[1:3], flip=(not flip_time_values)
298 )
300 # remove first element since it is duplicated zero
301 _np.delete(time2, 0)
302 _np.delete(force2, 0)
304 # add the respective force
305 time = _np.concatenate((time1, time2[1:]))
306 force = _np.concatenate((force1, force2[1:]), axis=0)
308 else:
309 raise ValueError(
310 "The selected type: "
311 + str(type)
312 + " is currently not supported. Feel free to add it here."
313 )
315 # overwrite the function, if one is provided since for the specific types the function is generated
316 if not type == "linear":
317 for dim in range(force.shape[1]):
318 # create a linear function with the force values per dimension
319 fun = _create_linear_interpolation_function(
320 time, force[:, dim], function_type="SYMBOLIC_FUNCTION_OF_TIME"
321 )
323 # add the function to the mesh
324 mesh.add(fun)
326 # store function
327 functions.append(fun)
329 # now set forces to 1 since the force values are already extracted in the function's values
330 force = _np.zeros_like(force) + 1.0
332 elif len(functions) != 3:
333 raise ValueError("Please provide functions with ")
335 # Create condition in mesh
336 add_point_neuman_condition_to_mesh(
337 mesh, nodes, functions, force[steps[1]], **kwargs
338 )
341def dbc_monitor_to_mesh(
342 mesh: _Mesh,
343 file_path: str,
344 *,
345 step: int = -1,
346 function: _Function,
347 **kwargs,
348):
349 """Converts the last value of a Dirichlet boundary condition monitor log to
350 a Neumann boundary condition in the mesh.
352 Args
353 ----
354 mesh: Mesh
355 The mesh where the created Neumann boundary condition is added to.
356 The nodes referred to in the log file have to match with the ones
357 in the mesh. It is advisable to only call this function once
358 all nodes have been added to the mesh.
359 file_path: str
360 Path to the Dirichlet boundary condition log file.
361 step: int
362 Step values to be used. Default is -1, i.e. the last step.
363 function: Function
364 Function for the Neumann boundary condition.
365 """
367 # read the force
368 nodes, _, force, _ = read_dbc_monitor_file(file_path)
370 # The forces are the negative reactions at the Dirichlet boundaries.
371 force *= -1.0
373 # Create condition in mesh
374 add_point_neuman_condition_to_mesh(
375 mesh, nodes, [function] * 3, force[step], **kwargs
376 )