diff --git a/examples/14-lsdyna/01-lsdyna_beam.py b/examples/14-lsdyna/01-lsdyna_beam.py new file mode 100644 index 00000000000..1f588fe541a --- /dev/null +++ b/examples/14-lsdyna/01-lsdyna_beam.py @@ -0,0 +1,359 @@ +# Copyright (C) 2020 - 2025 ANSYS, Inc. and/or its affiliates. +# SPDX-License-Identifier: MIT +# +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +""" +.. _lsdyna_operators: + +Beam results manipulations +-------------------------- + +This example provides an overview of the LS-DYNA beam results manipulations. + +.. note:: + This example requires DPF 6.1 (ansys-dpf-server-2023-2-pre0) or above. + For more information, see :ref:`ref_compatibility`. + +""" + +import matplotlib.pyplot as plt + +from ansys.dpf import core as dpf +from ansys.dpf.core import examples, operators as ops + +############################################################################### +# d3plot file data extraction +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Create the model and print its contents. This LS-DYNA d3plot file contains +# several individual results, each at different times. The d3plot file does not +# contain information related to Units. +# +# In this case, as the simulation was run through Mechanical, a ''file.actunits'' +# file is produced. If this file is supplemented in the data_sources, the units +# will be correctly fetched for all results in the file as well as for the mesh. + +d3plot = examples.download_d3plot_beam() +my_data_sources = dpf.DataSources() +my_data_sources.set_result_file_path(d3plot[0], key="d3plot") +my_data_sources.add_file_path(d3plot[3], key="actunits") +my_model = dpf.Model(my_data_sources) +print(my_model) + +############################################################################### +# Exploring the mesh +# ~~~~~~~~~~~~~~~~~~ +# +# The model has solid (3D) elements and beam (1D) elements. Some of the results +# only apply to one type of elements (such as the stress tensor for solids, or +# the axial force for beams, for example). +# +# By splitting the mesh by element shape we see that the ball is made by the solid +# 3D elements and the plate by the beam 1D elements +# +# - Define the analysis mesh +my_meshed_region = my_model.metadata.meshed_region + +# - Get separate meshes for each body +my_meshes = ops.mesh.split_mesh( + mesh=my_meshed_region, property=dpf.common.elemental_properties.element_shape +).eval() + +# - Define the meshes for each body in separate variables +ball_mesh = my_meshes.get_mesh(label_space_or_index={"body": 1, "elshape": 1}) +plate_mesh = my_meshes.get_mesh(label_space_or_index={"body": 2, "elshape": 2}) + +print(my_meshes) + +############################################################################### +# Plate mesh + +print("Plate mesh", "\n", plate_mesh) +plate_mesh.plot(title="Plate mesh", text="Plate mesh") + +############################################################################### +# Ball mesh + +print("Ball mesh", "\n", ball_mesh, "\n") +ball_mesh.plot(title="Ball mesh", text="Ball mesh") + +############################################################################### +# Scoping +# ~~~~~~~ +# +# - Define the mesh scoping to use it with the operators +my_meshes_scoping = ops.scoping.split_on_property_type(mesh=my_meshed_region).eval() + +############################################################################### +# - Define the mesh scoping for each body/element shape in separate variables +ball_scoping = my_meshes_scoping.get_scoping(label_space_or_index={"elshape": 1}) +plate_scoping = my_meshes_scoping.get_scoping(label_space_or_index={"elshape": 2}) + +############################################################################### +# - We will plot the results in a mesh deformed by the displacement. +# The displacement is in a nodal location, so we need to define a nodal scoping for the plate +plate_scoping_nodal = dpf.operators.scoping.transpose( + mesh_scoping=plate_scoping, meshed_region=my_meshed_region +).eval() + +############################################################################### +# Beam results +# ~~~~~~~~~~~~ +# The next manipulations can be applied to the following beam operators +# that handle the correspondent results : +# +# - beam_axial_force: Beam Axial Force +# - beam_s_shear_force: Beam S Shear Force +# - beam_t_shear_force: Beam T Shear Force +# - beam_s_bending_moment: Beam S Bending Moment +# - beam_t_bending_moment: Beam T Bending Moment +# - beam_torsional_moment: Beam Torsional Moment +# - beam_axial_stress: Beam Axial Stress +# - beam_rs_shear_stress: Beam Rs Shear Stress +# - beam_tr_shear_stress: Beam Tr Shear Stress +# - beam_axial_plastic_strain: Beam Axial Plastic Strain +# - beam_axial_total_strain: Beam Axial Total Strain +# +# We do not demonstrate separately how to use each of them in this example +# once they have similar methods. +# +# So, if you want to operate on other operator, uou just need to change their +# scripting name in the code lines. + +############################################################################### +# Comparing results in different time steps +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# 1) Define the time steps set +time_steps_set = [2, 6, 12] + +# 2) Prepare the collections to store the results for each time step + +# a. To compare the results in the same image you have to copy the mesh for each plot +plate_meshes = dpf.MeshesContainer() +plate_meshes.add_label("time") + +# b. The displacements for each time steps to deform the mesh accordingly +plate_displacements = dpf.FieldsContainer() +plate_displacements.add_label(label="time") + +# c. The axial force results for each time steps. Here +plate_axial_force = dpf.FieldsContainer() +plate_axial_force.add_label(label="time") + +# 3) Use the Plotter class to add the plots in the same image +comparison_plot = dpf.plotter.DpfPlotter() + +# Side bar arguments definition +side_bar_args = dict( + title="Beam axial force (N)", fmt="%.2e", title_font_size=15, label_font_size=15 +) + +# 4) As we want to compare the results in the same plot we will need this variable. +# It represents the distance between the meshes +j = -400 + +# 5) Copy the mesh of interest. Here it is the plate mesh that we copy along the X axis +# Here we use a loop where each iteration correspond to the manipulations for a given time step + +for i in time_steps_set: # Loop through the time steps + # Copy the mesh + plate_meshes.add_mesh(label_space={"time": i}, mesh=plate_mesh.deep_copy()) + + # 6) Get the plot coordinates that will be changed (so we can compare the results side by side) + coords_to_update = plate_meshes.get_mesh( + label_space_or_index={"time": i} + ).nodes.coordinates_field + + # 7) Define the coordinates where the new mesh will be placed + overall_field = dpf.fields_factory.create_3d_vector_field( + num_entities=1, location=dpf.locations.overall + ) + overall_field.append(data=[j, 0.0, 0.0], scopingid=1) + + # 8) Define the updated coordinates + new_coordinates = ops.math.add(fieldA=coords_to_update, fieldB=overall_field).eval() + coords_to_update.data = new_coordinates.data + + # 9) Extract the result, here we start by getting the beam_rs_shear_stress + plate_axial_force.add_field( + label_space={"time": i}, + field=my_model.results.beam_axial_force( + time_scoping=i, mesh_scoping=plate_scoping_nodal + ).eval()[0], + ) + # 10) We will also get the displacement to deform the mesh + plate_displacements.add_field( + label_space={"time": i}, + field=my_model.results.displacement( + time_scoping=i, mesh_scoping=plate_scoping_nodal + ).eval()[0], + ) + # 11) Add the result and the mesh to the plot + comparison_plot.add_field( + field=plate_axial_force.get_field(label_space_or_index={"time": i}), + meshed_region=plate_meshes.get_mesh(label_space_or_index={"time": i}), + deform_by=plate_displacements.get_field(label_space_or_index={"time": i}), + scalar_bar_args=side_bar_args, + ) + comparison_plot.add_node_labels( + nodes=[289], + labels=[f"Time step = {i}"], + meshed_region=plate_meshes.get_mesh(label_space_or_index={"time": i}), + font_size=10, + ) + # 12) Increment the coordinate value for the loop + j = j - 400 + + +# Visualise the plot +comparison_plot.show_figure() + +############################################################################### +# Plot a graph over time for the elements with max and min results values +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Here we make a workflow with a more verbose approach. This is useful because we use operators +# having several matching inputs or outputs. So the connexions are more clear, and it is +# easier to use and reuse the workflow. +# +# The following workflow finds the element with the max values over all the time steps and return its ID + +# Define the workflow object +max_workflow = dpf.Workflow() +max_workflow.progress_bar = False +# Define the norm operator +max_norm = ops.math.norm_fc() +# Define the max of each entity with the evaluated norm as an input +max_per_ent = ops.min_max.min_max_by_entity(fields_container=max_norm.outputs.fields_container) +# Define the max over all entities +global_max = ops.min_max.min_max(field=max_per_ent.outputs.field_max) +# Get the scoping +max_scop = ops.utility.extract_scoping(field_or_fields_container=global_max.outputs.field_max) +# Get the id +max_id = ops.scoping.scoping_get_attribute( + scoping=max_scop.outputs.mesh_scoping_as_scoping, property_name="ids" +) + +# Add the operators to the workflow +max_workflow.add_operators(operators=[max_norm, max_per_ent, global_max, max_scop, max_id]) +max_workflow.set_input_name("fields_container", max_norm.inputs.fields_container) +max_workflow.set_output_name("max_id", max_id.outputs.property_as_vector_int32_) +max_workflow.set_output_name("max_entity_scoping", max_scop.outputs.mesh_scoping_as_scoping) + +############################################################################### +# Using the workflow to the stresses results on the plate: +# +# - Extract the results + +# Get all the time steps +time_all = my_model.metadata.time_freq_support.time_frequencies + +# Extract all the stresses results on the plate +plate_beam_axial_stress = my_model.results.beam_axial_stress( + time_scoping=time_all, mesh_scoping=plate_scoping +).eval() +plate_beam_rs_shear_stress = my_model.results.beam_rs_shear_stress( + time_scoping=time_all, mesh_scoping=plate_scoping +).eval() +plate_beam_tr_shear_stress = my_model.results.beam_tr_shear_stress( + time_scoping=time_all, mesh_scoping=plate_scoping +).eval() + +############################################################################### +# - As we will use the workflow for different results operators we group them and +# use a loop through the group. Here we prepare where the workflow outputs will be stored + +# List of operators to be used in the workflow +beam_stresses = [plate_beam_axial_stress, plate_beam_rs_shear_stress, plate_beam_tr_shear_stress] +graph_labels = [ + "Beam axial stress", + "Beam rs shear stress", + "Beam tr shear stress", +] + +# List of elements ids that we will get from the workflow +max_stress_elements_ids = [] + +# Scopings container +max_stress_elements_scopings = dpf.ScopingsContainer() +max_stress_elements_scopings.add_label("stress_result") + +############################################################################### +# - The following loop: +# a) Goes through each stress result and get the element id with maximum solicitation +# b) Re-escope the fields container to keep only the data for this element +# c) Plot a stress x time graph + +for j in range(0, len(beam_stresses)): # Loop through each stress result + # Use the pre-defined workflow to define the element with maximum solicitation + max_workflow.connect(pin_name="fields_container", inpt=beam_stresses[j]) + max_stress_elements_ids.append( + max_workflow.get_output(pin_name="max_id", output_type=dpf.types.vec_int) + ) + max_stress_elements_scopings.add_scoping( + label_space={"stress_result": j}, + scoping=max_workflow.get_output( + pin_name="max_entity_scoping", output_type=dpf.types.scoping + ), + ) + + # Re-scope the results to keep only the data for the identified element + beam_stresses[j] = ops.scoping.rescope_fc( + fields_container=beam_stresses[j], + mesh_scoping=max_stress_elements_scopings.get_scoping( + label_space_or_index={"stress_result": j} + ), + ).eval() + + # The d3plot file gives us fields containers labeled by time. So in each field we have the stress value in a + # given time for the chosen element. We need to rearrange the fields container into fields. + + beam_stresses[j] = ops.utility.merge_to_field_matrix(fields1=beam_stresses[j]).eval() + plt.plot( + time_all.data, + beam_stresses[j].data[0], + label=f"{graph_labels[j]}, element id:{max_stress_elements_ids[j][0]}", + ) + +# Graph formatting +plt.title("Beam stresses evolution") +plt.xlabel("Time (s)") +plt.ylabel("Beam stresses (MPa)") +plt.legend() +plt.show() + +############################################################################### +# Results coordinates system +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The general results are given in the Cartesian coordinates system by default. +# +# The beam results are given directly in the local directions as scalars. +# For example the beam stresses we have: +# +# - The axial stress, given in the beam axis +# - The stresses defined in the cross-section directions: tr stress in the transverse +# direction (t) and rs stress perpendicular to the tr direction (s). +# +# +# Unfortunately there are no operators for LS-DYNA files that directly allows you to: +# - Rotate results from local coordinate system to global coordinate system; +# - Extract the rotation matrix between the local and global coordinate systems;