Coverage for src/meshpy/four_c/header

1# The MIT License (MIT)

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 List as _List

26from typing import Union as _Union

28from meshpy.core.conf import mpy as _mpy

29from meshpy.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 )

101

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 )

115

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 )

129

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 )

145

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 )

162

163 if output_energy:

164 input_file["STRUCTURAL DYNAMIC"]["RESEVERYERGY"] = 1

165

166

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.

183

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 """

212

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"}})

220

221 if (

222 "BEAM CONTACT" not in input_file

223 or input_file["BEAM CONTACT"].get("MODELEVALUATOR") != "Standard"

224 ):

225 input_file.add({"BEAM CONTACT": {"MODELEVALUATOR": "Standard"}})

226

227 set_binning_strategy_section(

228 input_file,

229 **binning_parameters,

230 )

231

232 # Add the beam to solid volume mesh tying options.

233 bts_parameters = {}

234 if interaction_type == _mpy.bc.beam_to_solid_volume_meshtying:

235 bts_section_name = "BEAM INTERACTION/BEAM TO SOLID VOLUME MESHTYING"

236 elif interaction_type == _mpy.bc.beam_to_solid_surface_meshtying:

237 bts_section_name = "BEAM INTERACTION/BEAM TO SOLID SURFACE MESHTYING"

238 if coupling_type is not None:

239 bts_parameters["COUPLING_TYPE"] = coupling_type

240 else:

241 raise ValueError(

242 "Got wrong beam-to-solid mesh tying type. "

243 f"Got {interaction_type} of type {type(interaction_type)}."

244 )

245 bts_parameters["CONSTRAINT_STRATEGY"] = "penalty"

246 if penalty_parameter is not None:

247 bts_parameters["PENALTY_PARAMETER"] = penalty_parameter

248 bts_parameters["GAUSS_POINTS"] = n_gauss_points

249

250 if contact_discretization == "mortar":

251 bts_parameters["CONTACT_DISCRETIZATION"] = "mortar"

252 bts_parameters["MORTAR_SHAPE_FUNCTION"] = mortar_shape

253 segmentation_strategy = _get_segmentation_strategy(segmentation)

254 elif contact_discretization == "gp":

255 bts_parameters["CONTACT_DISCRETIZATION"] = "gauss_point_to_segment"

256 segmentation_strategy = _get_segmentation_strategy(segmentation)

257 elif contact_discretization == "circ":

258 bts_parameters["CONTACT_DISCRETIZATION"] = "gauss_point_cross_section"

259 bts_parameters["INTEGRATION_POINTS_CIRCUMFERENCE"] = n_integration_points_circ

260 segmentation_strategy = "gauss_point_projection_cross_section"

261 else:

262 raise ValueError(

263 f'Wrong contact_discretization "{contact_discretization}" given!'

264 )

265

266 bts_parameters["GEOMETRY_PAIR_STRATEGY"] = segmentation_strategy

267 bts_parameters["GEOMETRY_PAIR_SEGMENTATION_SEARCH_POINTS"] = (

268 segmentation_search_points

269 )

270 if interaction_type == _mpy.bc.beam_to_solid_volume_meshtying:

271 bts_parameters["COUPLE_RESTART_STATE"] = couple_restart

272

273 input_file.add({bts_section_name: bts_parameters})

274

275

276def set_header_static(

277 input_file,

278 *,

279 time_step=None,

280 n_steps=None,

281 total_time=None,

282 max_iter=20,

283 tol_residuum=1e-8,

284 tol_increment=1e-10,

285 load_lin=False,

286 write_bin=False,

287 write_stress="no",

288 write_strain="no",

289 prestress="None",

290 prestress_time=0,

291):

292 """Set the default parameters for a static structure analysis.

293

294 At least two of the three time stepping keyword arguments ["time_step",

295 "n_steps", "total_time"] have to be set.

296

297 Args

298 ----

299 input_file:

300 Input file that the options will be added to.

301 time_step: float

302 Time increment per step.

303 n_steps: int

304 Number of time steps.

305 total_time: float

306 Total time of simulation

307 max_iter: int

308 Maximal number of Newton iterations.

309 tol_residuum: float

310 Tolerance for the convergence of the residuum.

311 tol_increment: int

312 Tolerance for the convergence of the displacement increment.

313 load_lin: bool

314 If the load_lin option should be set.

315 write_bin: bool

316 If binary output should be written.

317 write_stress: string

318 If and which stress output to write

319 write_strain: string

320 If and which strain output to write

321 prestress: string

322 Type of prestressing strategy to be used

323 presetrss_time: int

324 Prestress Time

325 """

326

327 # Set the parameters for a static analysis.

328 input_file.add(

329 {

330 "PROBLEM TYPE": {

331 "PROBLEMTYPE": "Structure",

332 }

333 }

334 )

335

336 input_file.add(

337 {

338 "IO": {

339 "OUTPUT_BIN": write_bin,

340 "STRUCT_DISP": False,

341 "STRUCT_STRESS": write_stress,

342 "STRUCT_STRAIN": write_strain,

343 "VERBOSITY": "Standard",

344 }

345 }

346 )

347

348 # Set the time step parameters

349 given_time_arguments = sum(

350 1 for arg in (time_step, n_steps, total_time) if arg is not None

351 )

352 if given_time_arguments < 2:

353 raise ValueError(

354 'At least two of the following arguments "time_step", "n_steps" or '

355 '"total_time" are required'

356 )

357 if time_step is None:

358 time_step = total_time / n_steps

359 elif n_steps is None:

360 n_steps = round(total_time / time_step)

361 elif total_time is None:

362 total_time = time_step * n_steps

363

364 input_file.add(

365 {

366 "STRUCTURAL DYNAMIC": {

367 "LINEAR_SOLVER": 1,

368 "INT_STRATEGY": "Standard",

369 "DYNAMICTYPE": "Statics",

370 "PREDICT": "TangDis",

371 "PRESTRESS": prestress,

372 "PRESTRESSTIME": prestress_time,

373 "TIMESTEP": time_step,

374 "NUMSTEP": n_steps,

375 "MAXTIME": total_time,

376 "LOADLIN": load_lin,

377 }

378 }

379 )

380

381 input_file.add(

382 {

383 "SOLVER 1": {

384 "NAME": "Structure_Solver",

385 "SOLVER": "Superlu",

386 }

387 }

388 )

389

390 # Set the contents of the NOX xml file.

391 nox_xml_contents = f"""

392 <ParameterList name="Status Test">

393

394 <ParameterList name="Outer Status Test">

395 <Parameter name="Test Type" type="string" value="Combo"/>

396 <Parameter name="Combo Type" type="string" value="OR" />

397

398 <ParameterList name="Test 0">

399 <Parameter name="Test Type" type="string" value="Combo" />

400 <Parameter name="Combo Type" type="string" value="AND" />

401

402 <ParameterList name="Test 0">

403 <Parameter name="Test Type" type="string" value="NormF" />

404

405 <ParameterList name="Quantity 0">

406 <Parameter name="Quantity Type" type="string" value="Structure" />

407 <Parameter name="Tolerance Type" type="string" value="Absolute" />

408 <Parameter name="Tolerance" type="double" value="{tol_residuum}" />

409 <Parameter name="Norm Type" type="string" value="Two Norm" />

410 <Parameter name="Scale Type" type="string" value="Scaled" />

411 </ParameterList>

412 </ParameterList>

413

414

415 <ParameterList name="Test 1">

416 <Parameter name="Test Type" type="string" value="NormUpdate" />

417

418 <ParameterList name="Quantity 0">

419 <Parameter name="Quantity Type" type="string" value="Structure" />

420 <Parameter name="Tolerance Type" type="string" value="Absolute" />

421 <Parameter name="Tolerance" type="double" value="{tol_increment}" />

422 <Parameter name="Norm Type" type="string" value="Two Norm" />

423 <Parameter name="Scale Type" type="string" value="Scaled" />

424 </ParameterList>

425 </ParameterList>

426

427 </ParameterList>

428

429

430 <ParameterList name="Test 1">

431 <Parameter name="Test Type" type="string" value="MaxIters" />

432 <Parameter name="Maximum Iterations" type="int" value="{max_iter}" />

433 </ParameterList>

434 </ParameterList>

435 </ParameterList>

436 """

437

438 input_file.add(

439 {

440 "STRUCT NOX/Printing": {

441 "Error": True,

442 "Inner Iteration": False,

443 "Details": True,

444 "Linear Solver Details": True,

445 "Test Details": True,

446 }

447 }

448 )

449

450 # Set the xml content in the input file.

451 input_file.nox_xml_contents = nox_xml_contents

452

453

454def set_binning_strategy_section(

455 input_file: _InputFile,

456 binning_bounding_box: _Union[_List[int], None] = None,

457 binning_cutoff_radius: _Union[float, None] = None,

458):

459 """Set binning strategy in section of the input file.

460

461 Args

462 ----

463 input_file:

464 Input file that the options will be added to.

465 binning_bounding_box:

466 List with the limits of the bounding box.

467 binning_cutoff_radius:

468 Maximal influence radius of pair elements.

469 """

470

471 if binning_bounding_box is not None and binning_cutoff_radius is not None:

472 binning_bounding_box_string = " ".join(

473 [str(val) for val in binning_bounding_box]

474 )

475

476 input_file.add(

477 {

478 "BINNING STRATEGY": {

479 "BIN_SIZE_LOWER_BOUND": binning_cutoff_radius,

480 "DOMAINBOUNDINGBOX": binning_bounding_box_string,

481 }

482 }

483 )

484 elif [binning_bounding_box, binning_cutoff_radius].count(None) == 2:

485 return

486 else:

487 raise ValueError(

488 f"The variables binning_bounding_box {binning_bounding_box} and binning_cutoff_radius {binning_cutoff_radius} must both be set."

489 )

490

491

492def set_beam_interaction_section(

493 input_file: _InputFile,

494 *,

495 repartition_strategy: str = "everydt",

496 search_strategy: str = "bounding_volume_hierarchy",

497):

498 """Set beam interaction section in input file.

499

500 Args

501 ----

502 input_file:

503 Input file that the options will be added to.

504 repartition_strategy:

505 Type of employed repartitioning strategy

506 Options: "adaptive" or "everydt"

507 search_strategy:

508 Type of search strategy used for finding coupling pairs.

509 Options: "bruteforce_with_binning", "bounding_volume_hierarchy"

510 """

511

512 input_file.add(

513 {

514 "BEAM INTERACTION": {

515 "REPARTITIONSTRATEGY": repartition_strategy,

516 "SEARCH_STRATEGY": search_strategy,

517 }

518 }

519 )

520

521

522def set_beam_contact_runtime_output(

523 input_file: _InputFile, *, every_iteration: bool = False

524):

525 """Output the beam-to-beam contact forces and gaps with runtime output.

526

527 input_file:

528 Input file that the options will be added to.

529 every_iteration:

530 If output at every Newton iteration should be written.

531 """

532

533 input_file.add(

534 {

535 "BEAM CONTACT/RUNTIME VTK OUTPUT": {

536 "VTK_OUTPUT_BEAM_CONTACT": True,

537 "EVERY_ITERATION": every_iteration,

538 "INTERVAL_STEPS": 1,

539 "CONTACT_FORCES": True,

540 "GAPS": True,

541 }

542 }

543 )

544

545

546def set_beam_contact_section(

547 input_file: _InputFile,

548 *,

549 interaction_strategy: str = "penalty",

550 btb_penalty: float = 0,

551 btb_line_penalty: float = 0,

552 per_shift_angle: list[float] = [70, 80],

553 par_shift_angle: list[float] = [70, 80],

554 b_seg_angle: float = 12,

555 num_integration: int = 5,

556 penalty_law: str = "LinPosQuadPen",

557 penalty_regularization_g0: float = 0,

558 penalty_regularization_f0: float = 0,

559 penalty_regularization_c0: float = 0,

560 binning_parameters: dict = {},

561 beam_interaction_parameters: dict = {},

562):

563 """Set default beam contact section, for more and updated details see

564 respective input file within 4C. Parameters for set_binning_strategy and

565 set_beam_interaction may be forwarded as keyword arguments.

566

567 Args

568 ----

569 input_file:

570 Input file that the options will be added to.

571 interaction_strategy:

572 Type of employed solving strategy

573 Options: "none", "penalty" or "gmshonly"

574 btb_penalty: double

575 Penalty parameter for beam-to-beam point contact

576 btb_line_penalty:

577 Penalty parameter per unit length for beam-to-beam line contact

578 per_shift_angle:

579 Lower and upper shift angle (in degrees) for penalty scaling of large-angle-contact

580 par_shift_angle:

581 Lower and upper shift angle (in degrees) for penalty scaling of small-angle-contact

582 b_seg_angle:

583 Maximal angle deviation allowed for contact search segmentation

584 num_integration:

585 Number of integration intervals per element

586 penalty_law:

587 Penalty Law Options: "LinPen", "QuadPen", "LinNegQuadPen", "LinPosQuadPen", "LinPosCubPen", "LinPosDoubleQuadPen", "LinPosExpPen"

588 penalty_regularization_g0:

589 First penalty regularization parameter G0

590 penalty_regularization_f0:

591 Second penalty regularization parameter F0

592 penalty_regularization_c0:

593 Third penalty regularization parameter C0

594 binning_parameters:

595 Keyword parameters for the binning section

596 beam_interaction_parameters:

597 Keyword parameters for the beam-contact section

598 """

599

600 if len(per_shift_angle) != 2:

601 raise ValueError(

602 "Please provide lower and upper value of BEAMS_PERPSHIFTANGLE."

603 )

604

605 if len(par_shift_angle) != 2:

606 raise ValueError("Please provide lower and upper value of BEAMS_PARSHIFTANGLE.")

607

608 input_file.add(

609 {

610 "BEAM INTERACTION/BEAM TO BEAM CONTACT": {

611 "STRATEGY": interaction_strategy,

612 }

613 }

614 )

615

616 input_file.add(

617 {

618 "BEAM CONTACT": {

619 "MODELEVALUATOR": "standard",

620 "BEAMS_STRATEGY": "penalty",

621 "BEAMS_BTBPENALTYPARAM": btb_penalty,

622 "BEAMS_BTBLINEPENALTYPARAM": btb_line_penalty,

623 "BEAMS_SEGCON": True,

624 "BEAMS_PERPSHIFTANGLE1": per_shift_angle[0],

625 "BEAMS_PERPSHIFTANGLE2": per_shift_angle[1],

626 "BEAMS_PARSHIFTANGLE1": par_shift_angle[0],

627 "BEAMS_PARSHIFTANGLE2": par_shift_angle[1],

628 "BEAMS_SEGANGLE": b_seg_angle,

629 "BEAMS_NUMINTEGRATIONINTERVAL": num_integration,

630 "BEAMS_PENALTYLAW": penalty_law,

631 "BEAMS_PENREGPARAM_G0": penalty_regularization_g0,

632 "BEAMS_PENREGPARAM_F0": penalty_regularization_f0,

633 "BEAMS_PENREGPARAM_C0": penalty_regularization_c0,

634 "BEAMS_MAXDELTADISSCALEFAC": -1.0,

635 }

636 }

637 )

638

639 # beam contact needs a binning strategy

640 set_binning_strategy_section(input_file, **binning_parameters)

641

642 # beam contact needs interaction strategy

643 set_beam_interaction_section(input_file, **beam_interaction_parameters)

644

645

646def add_result_description(

647 input_file: _InputFile,

648 displacements: _List,

649 node_ids: _List[int],

650 *,

651 tol: float = 1e-10,

652):

653 """Add result descriptions for structure problems to the input file.

654

655 Args:

656 input_file: Input file to add the result description to

657 displacements: Array with the displacements (n_nodes x 3)

658 node_ids: List with the IDs of the nodes to check

659 tol: Tolerance

660 """

661 result_descriptions = []

662

663 for i_node, node in enumerate(node_ids):

664 for i_dir, direction in enumerate(["x", "y", "z"]):

665 result_descriptions.append(

666 {

667 "STRUCTURE": {

668 "DIS": "structure",

669 "NODE": node,

670 "QUANTITY": f"disp{direction}",

671 "VALUE": displacements[i_node][i_dir],

672 "TOLERANCE": tol,

673 },

674 }

675 )

676

677 input_file.add({"RESULT DESCRIPTION": result_descriptions})

Coverage for src/meshpy/four_c/header_functions.py: 90%

96 statements