How do I calculate linearized stress (Membrane, Bending and Peak) using DPF?

Ayush Kumar

How do I calculate linearized stress (Membrane, Bending and Peak) using DPF?

Ayush Kumar

import mech_dpf
import Ans.DataProcessing as dpf


path = ExtAPI.DataModel.GetObjectsByName("Path")[0]
x_0 = path.StartXCoordinate.Value
y_0 = path.StartYCoordinate.Value
z_0 = path.StartZCoordinate.Value


x_1 = path.EndXCoordinate.Value
y_1 = path.EndYCoordinate.Value
z_1 = path.EndZCoordinate.Value


path_length = ((x_1 - x_0) ** 2 + (y_1 - y_0) ** 2 + (z_1 - z_0) ** 2) ** 0.5


n_points = path.DiscretizationPoints + 2
delta = path_length / (n_points - 1)


line_unit_vector = [(x_1 - x_0) / path_length, (y_1 - y_0) / path_length, (z_1 - z_0) / path_length]


# Line equation
fx = lambda t: x_0 + line_unit_vector[0] * t
fy = lambda t: y_0 + line_unit_vector[1] * t
fz = lambda t: z_0 + line_unit_vector[2] * t


coordinates = [[fx(t * delta), fy(t * delta), fz(t * delta)] for t in range(n_points)]
flat_coordinates = [entry for data in coordinates for entry in data]


field_coord = dpf.FieldsFactory.Create3DVectorField(n_points)
field_coord.Data = flat_coordinates
field_coord.Scoping.Ids = list(range(1, n_points + 1))


model = dpf.Model(Model.Analyses[0].ResultFileName)


ts = dpf.Scoping()
ts.Ids = [2]


# SX
s = model.CreateOperator("SX")
s.inputs.requested_location.Connect("Nodal")
s.inputs.time_scoping.Connect(ts)
s2_f = s.outputs.fields_container.GetData()


mapping_operator = dpf.operators.mapping.on_coordinates(
    fields_container=s2_f,
    coordinates=field_coord,
    create_support=True,
    mesh=model.Mesh
)
fields_mapped = mapping_operator.outputs.fields_container.GetData()


ls = list(fields_mapped[0].Data)
# Membrane Stress
membrane_stress = (ls[0] / 2 + ls[-1] / 2 + sum(ls[1:-1])) / (n_points - 1)


# Bending stress
path_1 = -1 * path_length / 2
path_n = path_length / 2
path_range = [path_1 + delta * i for i in range(n_points)]
path_range_field = dpf.FieldsFactory.CreateScalarField(numEntities=n_points, location=dpf.locations.nodal)
path_range_field.Data = path_range
path_range_field.ScopingIds = range(1, 50)


# Function to be integrated 
stress_scaled = dpf.operators.math.scale_by_field(fieldA=fields_mapped[0], fieldB=path_range_field)
stress_scaled_data = list(stress_scaled.outputs.field.GetData().Data)


# Use extended Simpson's rule for Numerical Integration of `stress_scaled_data`
stress_scaled_integral = (17 * stress_scaled_data[0] + 59 * stress_scaled_data[1] + 43 * stress_scaled_data[2] + 49 *
                          stress_scaled_data[3] + 48 * sum(stress_scaled_data[4:-4]) + 49 * stress_scaled_data[
                              n_points - 4] + 43 * stress_scaled_data[n_points - 3] + 59 * stress_scaled_data[
                              n_points - 2] + 17 * stress_scaled_data[n_points - 1]) / 48.0


b1 = stress_scaled_integral * (-6.0 / path_length) / 48.0
b2 = (-1.0) * b1
<?WORD FLAGGED?>_bending = (b2 - b1) / (path_length)
fb = lambda t: b1 + <?WORD FLAGGED?>_bending * t
bending_stress = [fb(t * delta) for t in range(n_points)]


bend_f = dpf.FieldsFactory.CreateScalarField(numEntities=49, location=dpf.locations.nodal)
bend_f.Data = bending_stress
bend_f.Scoping.Ids = list(range(1, n_points + 1))


# Add Membrane and Bending Stress
bend_mem_f = dpf.operators.math.add_constant(field=bend_f, ponderation=membrane_stress)
bend_mem_f_neg = dpf.operators.math.scale(field=bend_mem_f, ponderation=-1.0)


# Peak stress
ps = dpf.operators.math.add(fieldA=fields_mapped[0], fieldB=bend_mem_f_neg)

dafedin

Hello Ayush,

Thanks for this comprehensive example. I wonder if there is a way to solve the reversed task: find elements of the model whose results were used for mapping to path points. Maybe there is another option rather than screening over the whole mesh and finding elements whose centroids are situated in the path direction.

Regards

Ayush Kumar

Hi @dafedin,

That information is not provided as an output of the operator at the moment.

Jim Kosloski

Is there a more efficient way than looping where we could get these results (bending and membrane) for ALL timepoints? Most of our DPF examples focus on a single time point, but one of the big benefits of DPF is the ability to process multiple timepoints simultaneously. All of these calculations use fields_mapped[0].Data how can we use all the data in the fields_mapped container?

Also how can I use this in a Python Result and plot the path results on the path. If I send bend_f to the dpf_workflow.Add I get an error.

Jim Kosloski

Is there a more efficient way than looping where we could get these results (bending and membrane) for ALL timepoints? Most of our DPF examples focus on a single time point, but one of the big benefits of DPF is the ability to process multiple timepoints simultaneously. All of these calculations use fields_mapped[0].Data how can we use all the data in the fields_mapped container?

Also how can I use this in a Python Result and plot the path results on the path. If I send bend_f to the dpf_workflow.Add I get an error.

Jim Kosloski

@Ayush Kumar a few more questions on this. 1. What does the math.scale_by_field operator do? I do not see it documented anywhere.
2. Do the field coordinates input into the mapping operator need to be in SI?

Ayush Kumar

@Jim Kosloski
1. As most of the calculation is done on the Client side, I don't see a DPF way to calculate for all time-steps at once (in Mechanical DPF). For PyDPF, NumPy can be used for vectorized operations.
2. You can try forwarding the field to an operator (dpf.operators.utility.forward()) and then pass that operator to workflow. Also make sure you have the new mesh defined using SetOutputMesh(...).
3. scale_by_field: https://dpf.docs.pyansys.com/version/stable/api/ansys.dpf.core.operators.math.scale_by_field.html
4. Coordinate field input in the mapping operator is not unit dependent.