Coverage for src/beamme/four_c/header_functions.py: 90%
99 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 module defines functions that can be used to add header information to
23an input file."""
25from typing import Any as _Any
26from typing import List as _List
28from beamme.core.conf import bme as _bme
29from beamme.four_c.input_file import InputFile as _InputFile
32def _get_segmentation_strategy(segmentation):
33 """Get the 4C string for a geometry pair strategy."""
34 if segmentation:
35 return "segmentation"
36 else:
37 return "gauss_point_projection_without_boundary_segmentation"
40def set_runtime_output(
41 input_file,
42 *,
43 output_solid=True,
44 output_stress_strain=False,
45 btsvmt_output=True,
46 btss_output=True,
47 output_triad=True,
48 every_iteration=False,
49 absolute_beam_positions=True,
50 element_owner=True,
51 element_gid=True,
52 element_mat_id=True,
53 output_energy=False,
54 output_strains=True,
55):
56 """Set the basic runtime output options.
58 Args
59 ----
60 input_file:
61 Input file that the options will be added to.
62 output_solid: bool
63 If the solid output should be written at runtime.
64 output_stress_strain: bool
65 If stress and strain output should be written for the solid.
66 btsvmt_output: bool
67 If the output for btsvmt should be written.
68 btss_output: bool
69 If the output for beam-to-surface coupling should be written.
70 output_triad: bool
71 If the triads along the beam should be written.
72 every_iteration: int
73 If output at every Newton iteration should be written.
74 absolute_beam_positions: bool
75 If the beams should be written at the current position or always at
76 the reference position.
77 element_owner: bool
78 If the owing rank of each element should be output (currently
79 only affects the solid elements in 4C, beam element owners are
80 written by default).
81 element_gid: bool
82 If the 4C internal GID of each element should be output.
83 element_mat_id: bool
84 If the 4C internal material ID of each element should be output.
85 output_energy: bool
86 If the energy output from 4C should be activated.
87 output_strains: bool
88 If the strains in the Gauss points should be output.
89 """
91 # Set the basic runtime output options.
92 input_file.add(
93 {
94 "IO/RUNTIME VTK OUTPUT": {
95 "OUTPUT_DATA_FORMAT": "binary",
96 "INTERVAL_STEPS": 1,
97 "EVERY_ITERATION": every_iteration,
98 }
99 }
100 )
102 # Set the structure runtime output options
103 input_file.add(
104 {
105 "IO/RUNTIME VTK OUTPUT/STRUCTURE": {
106 "OUTPUT_STRUCTURE": output_solid,
107 "DISPLACEMENT": True,
108 "STRESS_STRAIN": output_stress_strain,
109 "ELEMENT_OWNER": element_owner,
110 "ELEMENT_GID": element_gid,
111 "ELEMENT_MAT_ID": element_mat_id,
112 }
113 }
114 )
116 # Set the beam runtime output options
117 input_file.add(
118 {
119 "IO/RUNTIME VTK OUTPUT/BEAMS": {
120 "OUTPUT_BEAMS": True,
121 "DISPLACEMENT": True,
122 "USE_ABSOLUTE_POSITIONS": absolute_beam_positions,
123 "TRIAD_VISUALIZATIONPOINT": output_triad,
124 "STRAINS_GAUSSPOINT": output_strains,
125 "ELEMENT_GID": element_gid,
126 }
127 }
128 )
130 if btsvmt_output:
131 # Set the beam to solid volume mesh tying runtime output options.
132 input_file.add(
133 {
134 "BEAM INTERACTION/BEAM TO SOLID VOLUME MESHTYING/RUNTIME VTK OUTPUT": {
135 "WRITE_OUTPUT": True,
136 "NODAL_FORCES": True,
137 "MORTAR_LAMBDA_DISCRET": True,
138 "MORTAR_LAMBDA_CONTINUOUS": True,
139 "MORTAR_LAMBDA_CONTINUOUS_SEGMENTS": 5,
140 "SEGMENTATION": True,
141 "INTEGRATION_POINTS": True,
142 }
143 }
144 )
146 if btss_output:
147 # Set the beam to solid surface coupling runtime output options.
148 input_file.add(
149 {
150 "BEAM INTERACTION/BEAM TO SOLID SURFACE/RUNTIME VTK OUTPUT": {
151 "WRITE_OUTPUT": True,
152 "NODAL_FORCES": True,
153 "MORTAR_LAMBDA_DISCRET": True,
154 "MORTAR_LAMBDA_CONTINUOUS": True,
155 "MORTAR_LAMBDA_CONTINUOUS_SEGMENTS": 5,
156 "SEGMENTATION": True,
157 "INTEGRATION_POINTS": True,
158 "AVERAGED_NORMALS": True,
159 }
160 }
161 )
163 if output_energy:
164 input_file["STRUCTURAL DYNAMIC"]["RESEVERYERGY"] = 1
167def set_beam_to_solid_meshtying(
168 input_file,
169 interaction_type,
170 *,
171 contact_discretization=None,
172 segmentation=True,
173 segmentation_search_points=2,
174 couple_restart=False,
175 mortar_shape="none",
176 n_gauss_points=6,
177 n_integration_points_circ=None,
178 penalty_parameter=None,
179 coupling_type=None,
180 binning_parameters: dict = {},
181):
182 """Set the beam to solid meshtying options.
184 Args
185 ----
186 input_file:
187 Input file that the options will be added to.
188 interaction_type: BoundaryCondition
189 Type of beam-to-solid interaction.
190 contact_discretization: str
191 Type of contact (mortar, Gauss point, ...)
192 segmentation: bool
193 If segmentation should be used in the numerical integration.
194 segmentation_search_points: int
195 Number of search points for segmentation.
196 couple_restart: bool
197 If the restart configuration should be used for the coupling
198 mortar_shape: str
199 Type of shape function for mortar discretization.
200 n_gauss_points: int
201 Number of Gauss points for numerical integration.
202 n_integration_points_circ: int
203 Number of integration points along the circumference of the cross
204 section.
205 penalty_parameter: float
206 Penalty parameter for contact enforcement.
207 coupling_type: str
208 Type of coupling for beam-to-surface coupling.
209 binning_parameters:
210 Keyword parameters for the binning section
211 """
213 # Set the beam contact options.
214 # check if these keys are already set, otherwise set them
215 if (
216 "BEAM INTERACTION" not in input_file
217 or input_file["BEAM INTERACTION"].get("REPARTITIONSTRATEGY") != "everydt"
218 ):
219 input_file.add({"BEAM INTERACTION": {"REPARTITIONSTRATEGY": "everydt"}})
221 set_binning_strategy_section(
222 input_file,
223 **binning_parameters,
224 )
226 # Add the beam to solid volume mesh tying options.
227 bts_parameters = {}
228 if interaction_type == _bme.bc.beam_to_solid_volume_meshtying:
229 bts_section_name = "BEAM INTERACTION/BEAM TO SOLID VOLUME MESHTYING"
230 elif interaction_type == _bme.bc.beam_to_solid_surface_meshtying:
231 bts_section_name = "BEAM INTERACTION/BEAM TO SOLID SURFACE MESHTYING"
232 if coupling_type is not None:
233 bts_parameters["COUPLING_TYPE"] = coupling_type
234 else:
235 raise ValueError(
236 "Got wrong beam-to-solid mesh tying type. "
237 f"Got {interaction_type} of type {type(interaction_type)}."
238 )
239 bts_parameters["CONSTRAINT_STRATEGY"] = "penalty"
240 if penalty_parameter is not None:
241 bts_parameters["PENALTY_PARAMETER"] = penalty_parameter
242 bts_parameters["GAUSS_POINTS"] = n_gauss_points
244 if contact_discretization == "mortar":
245 bts_parameters["CONTACT_DISCRETIZATION"] = "mortar"
246 bts_parameters["MORTAR_SHAPE_FUNCTION"] = mortar_shape
247 segmentation_strategy = _get_segmentation_strategy(segmentation)
248 elif contact_discretization == "gp":
249 bts_parameters["CONTACT_DISCRETIZATION"] = "gauss_point_to_segment"
250 segmentation_strategy = _get_segmentation_strategy(segmentation)
251 elif contact_discretization == "circ":
252 bts_parameters["CONTACT_DISCRETIZATION"] = "gauss_point_cross_section"
253 bts_parameters["INTEGRATION_POINTS_CIRCUMFERENCE"] = n_integration_points_circ
254 segmentation_strategy = "gauss_point_projection_cross_section"
255 else:
256 raise ValueError(
257 f'Wrong contact_discretization "{contact_discretization}" given!'
258 )
260 bts_parameters["GEOMETRY_PAIR_STRATEGY"] = segmentation_strategy
261 bts_parameters["GEOMETRY_PAIR_SEGMENTATION_SEARCH_POINTS"] = (
262 segmentation_search_points
263 )
264 if interaction_type == _bme.bc.beam_to_solid_volume_meshtying:
265 bts_parameters["COUPLE_RESTART_STATE"] = couple_restart
267 input_file.add({bts_section_name: bts_parameters})
270def set_header_static(
271 input_file: _InputFile,
272 *,
273 time_step: float | None = None,
274 n_steps: int | None = None,
275 total_time: float | None = None,
276 max_iter: int = 20,
277 tol_residuum: float = 1e-8,
278 tol_increment: float = 1e-10,
279 restart_every: int = 1,
280 load_lin: bool = False,
281 write_bin: bool = False,
282 write_stress: str = "no",
283 write_strain: str = "no",
284 predictor: str = "TangDis",
285 prestress: str = "None",
286 prestress_time: float = 0,
287 create_nox_file: bool = True,
288):
289 """Set the default parameters for a static structure analysis.
291 At least two of the three time stepping keyword arguments ["time_step",
292 "n_steps", "total_time"] have to be set.
294 Args:
295 input_file:
296 Input file that the options will be added to.
297 time_step:
298 Time increment per step.
299 n_steps:
300 Number of time steps.
301 total_time:
302 Total time of simulation
303 max_iter:
304 Maximal number of Newton iterations.
305 tol_residuum:
306 Tolerance for the convergence of the residuum.
307 tol_increment:
308 Tolerance for the convergence of the displacement increment.
309 load_lin:
310 If the load_lin option should be set.
311 write_bin:
312 If binary output should be written.
313 restart_every:
314 Frequency for writing restart output.
315 write_stress:
316 If and which stress output to write
317 write_strain:
318 If and which strain output to write
319 predictor:
320 Type of predictor to be used
321 prestress:
322 Type of prestressing strategy to be used
323 prestress_time:
324 Prestress Time
325 create_nox_file:
326 If the nonlinear solver parameters should be set via a NOX xml file or
327 directly in the input file.
328 """
330 input_file_parameters: dict[str, _Any] = {}
332 # Set the parameters for a static analysis.
333 input_file_parameters["PROBLEM TYPE"] = {"PROBLEMTYPE": "Structure"}
334 input_file_parameters["IO"] = {
335 "OUTPUT_BIN": write_bin,
336 "STRUCT_DISP": False,
337 "STRUCT_STRESS": write_stress,
338 "STRUCT_STRAIN": write_strain,
339 "VERBOSITY": "Standard",
340 }
342 # Set the time step parameters
343 given_time_arguments = sum(
344 1 for arg in (time_step, n_steps, total_time) if arg is not None
345 )
346 if given_time_arguments < 2:
347 raise ValueError(
348 'At least two of the following arguments "time_step", "n_steps" or '
349 '"total_time" are required'
350 )
351 elif time_step is None and total_time is not None and n_steps is not None:
352 time_step = total_time / n_steps
353 elif n_steps is None and total_time is not None and time_step is not None:
354 n_steps = round(total_time / time_step)
355 elif total_time is None and time_step is not None and n_steps is not None:
356 total_time = time_step * n_steps
358 input_file_parameters["STRUCTURAL DYNAMIC"] = {
359 "LINEAR_SOLVER": 1,
360 "INT_STRATEGY": "Standard",
361 "DYNAMICTYPE": "Statics",
362 "PREDICT": predictor,
363 "PRESTRESS": prestress,
364 "PRESTRESSTIME": prestress_time,
365 "TIMESTEP": time_step,
366 "NUMSTEP": n_steps,
367 "MAXTIME": total_time,
368 "LOADLIN": load_lin,
369 "RESTARTEVERY": restart_every,
370 }
371 input_file_parameters["SOLVER 1"] = {
372 "NAME": "Structure_Solver",
373 "SOLVER": "Superlu",
374 }
376 # Set the solver parameters.
377 if create_nox_file:
378 # Set the contents of the NOX xml file.
379 nox_xml_contents = f"""
380 <ParameterList name="Status Test">
381 <!-- Outer Status Test: This test is an OR combination of the structural convergence and the maximum number of iterations -->
382 <ParameterList name="Outer Status Test">
383 <Parameter name="Test Type" type="string" value="Combo"/>
384 <Parameter name="Combo Type" type="string" value="OR" />
385 <!-- Structural convergence is an AND combination of the residuum and step update -->
386 <ParameterList name="Test 0">
387 <Parameter name="Test Type" type="string" value="Combo" />
388 <Parameter name="Combo Type" type="string" value="AND" />
389 <!-- BEGIN: Combo AND - Test 0: "NormF" -->
390 <ParameterList name="Test 0">
391 <Parameter name="Test Type" type="string" value="NormF" />
392 <!-- NormF - Quantity 0: Check the right-hand-side norm of the structural quantities -->
393 <ParameterList name="Quantity 0">
394 <Parameter name="Quantity Type" type="string" value="Structure" />
395 <Parameter name="Tolerance Type" type="string" value="Absolute" />
396 <Parameter name="Tolerance" type="double" value="{tol_residuum}" />
397 <Parameter name="Norm Type" type="string" value="Two Norm" />
398 <Parameter name="Scale Type" type="string" value="Scaled" />
399 </ParameterList>
400 </ParameterList>
401 <!-- END: Combo AND - Test 0: "NormF" -->
402 <!-- BEGIN: Combo AND - Test 1: "NormWRMS" -->
403 <ParameterList name="Test 1">
404 <Parameter name="Test Type" type="string" value="NormUpdate" />
405 <!-- NormWRMS - Quantity 0: Check the increment of the structural displacements -->
406 <ParameterList name="Quantity 0">
407 <Parameter name="Quantity Type" type="string" value="Structure" />
408 <Parameter name="Tolerance Type" type="string" value="Absolute" />
409 <Parameter name="Tolerance" type="double" value="{tol_increment}" />
410 <Parameter name="Norm Type" type="string" value="Two Norm" />
411 <Parameter name="Scale Type" type="string" value="Scaled" />
412 </ParameterList>
413 </ParameterList>
414 <!-- END: Combo AND - Test 1: "NormWRMS" -->
415 </ParameterList>
416 <!-- END: Combo 0 - Test 0: "Combo" -->
417 <!-- BEGIN: Combo OR - Test 1: "MaxIters" -->
418 <ParameterList name="Test 1">
419 <Parameter name="Test Type" type="string" value="MaxIters" />
420 <Parameter name="Maximum Iterations" type="int" value="{max_iter}" />
421 </ParameterList> <!--END: "MaxIters" -->
422 </ParameterList>
423 </ParameterList>
424 """
426 input_file_parameters["STRUCT NOX/Printing"] = {
427 "Error": True,
428 "Inner Iteration": False,
429 "Details": True,
430 "Linear Solver Details": True,
431 "Test Details": True,
432 }
434 # Set the xml content in the input file.
435 input_file.nox_xml_contents = nox_xml_contents
437 else:
438 input_file_parameters["STRUCTURAL DYNAMIC"]["MAXITER"] = max_iter
439 input_file_parameters["STRUCTURAL DYNAMIC"]["TOLRES"] = tol_residuum
440 input_file_parameters["STRUCTURAL DYNAMIC"]["TOLDISP"] = tol_increment
442 input_file.add(input_file_parameters)
445def set_binning_strategy_section(
446 input_file: _InputFile,
447 binning_bounding_box: list[int] | None = None,
448 binning_cutoff_radius: float | None = None,
449):
450 """Set binning strategy in section of the input file.
452 Args
453 ----
454 input_file:
455 Input file that the options will be added to.
456 binning_bounding_box:
457 List with the limits of the bounding box.
458 binning_cutoff_radius:
459 Maximal influence radius of pair elements.
460 """
462 if binning_bounding_box is not None and binning_cutoff_radius is not None:
463 binning_bounding_box_string = " ".join(
464 [str(val) for val in binning_bounding_box]
465 )
467 input_file.add(
468 {
469 "BINNING STRATEGY": {
470 "BIN_SIZE_LOWER_BOUND": binning_cutoff_radius,
471 "DOMAINBOUNDINGBOX": binning_bounding_box_string,
472 }
473 }
474 )
475 elif [binning_bounding_box, binning_cutoff_radius].count(None) == 2:
476 return
477 else:
478 raise ValueError(
479 f"The variables binning_bounding_box {binning_bounding_box} and binning_cutoff_radius {binning_cutoff_radius} must both be set."
480 )
483def set_beam_interaction_section(
484 input_file: _InputFile, *, repartition_strategy: str = "everydt"
485):
486 """Set beam interaction section in input file.
488 Args
489 ----
490 input_file:
491 Input file that the options will be added to.
492 repartition_strategy:
493 Type of employed repartitioning strategy
494 Options: "adaptive" or "everydt"
495 search_strategy:
496 Type of search strategy used for finding coupling pairs.
497 """
499 input_file.add({"BEAM INTERACTION": {"REPARTITIONSTRATEGY": repartition_strategy}})
502def set_beam_contact_runtime_output(
503 input_file: _InputFile, *, every_iteration: bool = False
504):
505 """Output the beam-to-beam contact forces and gaps with runtime output.
507 input_file:
508 Input file that the options will be added to.
509 every_iteration:
510 If output at every Newton iteration should be written.
511 """
513 input_file.add(
514 {
515 "BEAM INTERACTION/BEAM TO BEAM CONTACT/RUNTIME VTK OUTPUT": {
516 "VTK_OUTPUT_BEAM_CONTACT": True,
517 "EVERY_ITERATION": every_iteration,
518 "INTERVAL_STEPS": 1,
519 "CONTACT_FORCES": True,
520 "GAPS": True,
521 }
522 }
523 )
526def set_beam_contact_section(
527 input_file: _InputFile,
528 *,
529 interaction_strategy: str = "penalty",
530 btb_penalty: float = 0,
531 btb_line_penalty: float = 0,
532 per_shift_angle: list[float] = [70, 80],
533 par_shift_angle: list[float] = [70, 80],
534 b_seg_angle: float = 12,
535 num_integration: int = 5,
536 penalty_law: str = "LinPosQuadPen",
537 penalty_regularization_g0: float = 0,
538 penalty_regularization_f0: float = 0,
539 penalty_regularization_c0: float = 0,
540 binning_parameters: dict = {},
541 beam_interaction_parameters: dict = {},
542):
543 """Set default beam contact section, for more and updated details see
544 respective input file within 4C. Parameters for set_binning_strategy and
545 set_beam_interaction may be forwarded as keyword arguments.
547 Args
548 ----
549 input_file:
550 Input file that the options will be added to.
551 interaction_strategy:
552 Type of employed solving strategy
553 Options: "none", "penalty" or "gmshonly"
554 btb_penalty: double
555 Penalty parameter for beam-to-beam point contact
556 btb_line_penalty:
557 Penalty parameter per unit length for beam-to-beam line contact
558 per_shift_angle:
559 Lower and upper shift angle (in degrees) for penalty scaling of large-angle-contact
560 par_shift_angle:
561 Lower and upper shift angle (in degrees) for penalty scaling of small-angle-contact
562 b_seg_angle:
563 Maximal angle deviation allowed for contact search segmentation
564 num_integration:
565 Number of integration intervals per element
566 penalty_law:
567 Penalty Law Options: "LinPen", "QuadPen", "LinNegQuadPen", "LinPosQuadPen", "LinPosCubPen", "LinPosDoubleQuadPen", "LinPosExpPen"
568 penalty_regularization_g0:
569 First penalty regularization parameter G0
570 penalty_regularization_f0:
571 Second penalty regularization parameter F0
572 penalty_regularization_c0:
573 Third penalty regularization parameter C0
574 binning_parameters:
575 Keyword parameters for the binning section
576 beam_interaction_parameters:
577 Keyword parameters for the beam-contact section
578 """
580 if len(per_shift_angle) != 2:
581 raise ValueError(
582 "Please provide lower and upper value of BEAMS_PERPSHIFTANGLE."
583 )
585 if len(par_shift_angle) != 2:
586 raise ValueError("Please provide lower and upper value of BEAMS_PARSHIFTANGLE.")
588 input_file.add(
589 {
590 "BEAM INTERACTION/BEAM TO BEAM CONTACT": {
591 "BEAMS_BTBPENALTYPARAM": btb_penalty,
592 "BEAMS_BTBLINEPENALTYPARAM": btb_line_penalty,
593 "BEAMS_SEGCON": True,
594 "BEAMS_PERPSHIFTANGLE1": per_shift_angle[0],
595 "BEAMS_PERPSHIFTANGLE2": per_shift_angle[1],
596 "BEAMS_PARSHIFTANGLE1": par_shift_angle[0],
597 "BEAMS_PARSHIFTANGLE2": par_shift_angle[1],
598 "BEAMS_SEGANGLE": b_seg_angle,
599 "BEAMS_NUMINTEGRATIONINTERVAL": num_integration,
600 "BEAMS_PENALTYLAW": penalty_law,
601 "BEAMS_PENREGPARAM_G0": penalty_regularization_g0,
602 "BEAMS_PENREGPARAM_F0": penalty_regularization_f0,
603 "BEAMS_PENREGPARAM_C0": penalty_regularization_c0,
604 }
605 }
606 )
608 # beam contact needs a binning strategy
609 set_binning_strategy_section(input_file, **binning_parameters)
611 # beam contact needs interaction strategy
612 set_beam_interaction_section(input_file, **beam_interaction_parameters)
615def add_result_description(
616 input_file: _InputFile,
617 displacements: _List,
618 node_ids: _List[int],
619 *,
620 tol: float = 1e-10,
621):
622 """Add result descriptions for structure problems to the input file.
624 Args:
625 input_file: Input file to add the result description to
626 displacements: Array with the displacements (n_nodes x 3)
627 node_ids: List with the IDs of the nodes to check
628 tol: Tolerance
629 """
630 result_descriptions = []
632 for i_node, node in enumerate(node_ids):
633 for i_dir, direction in enumerate(["x", "y", "z"]):
634 result_descriptions.append(
635 {
636 "STRUCTURE": {
637 "DIS": "structure",
638 "NODE": node,
639 "QUANTITY": f"disp{direction}",
640 "VALUE": displacements[i_node][i_dir],
641 "TOLERANCE": tol,
642 },
643 }
644 )
646 input_file.add({"RESULT DESCRIPTION": result_descriptions})