# How do I calculate stress along a path of desired length normal to a surface at a node?

Options
Member, Moderator, Employee Posts: 414
edited June 2023

How do I calculate stress along a path of desired length normal to a surface at a node?

Tagged:

• Member, Moderator, Employee Posts: 414
Options

With ACT you can only define path through the complete length of the body of the chosen node and interpolate to get value a certain point. With PyDPF can define path of desired length.

ACT - define a named selection `face` and `node`. Node where you want the path to start and the face to which the node belongs.

```# Create face and node named selections
geom_id = ExtAPI.DataModel.GetObjectsByName("face")[0].Ids[0]
node_id = ExtAPI.DataModel.GetObjectsByName("node")[0].Ids[0]

mesh_data = DataModel.MeshDataByName(DataModel.MeshDataNames[0])
my_face = DataModel.GeoData.GeoEntityById(geom_id)
mesh_node = mesh_data.NodeById(node_id)
base_point = mesh_node.X, mesh_node.Y, mesh_node.Z

# Get Normal vector at the Node
u, v = my_face.ParamAtPoint(base_point)
vec_normal_out = my_face.NormalAtParam(u, v)  # Normal Vector to the face at given node
vec_normal_in = [-vec_normal_out[0], -vec_normal_out[1], -vec_normal_out[2]]

vec_global = [1, 0, 0]  # Global X axis

# Cross Product function definition
def cross_product(vec_1, vec_2):
vec_cross_product = [vec_1[1] * vec_2[2] - vec_1[2] * vec_2[1],
vec_1[2] * vec_2[0] - vec_1[0] * vec_2[2],
vec_1[0] * vec_2[1] - vec_1[1] * vec_2[0]]
return vec_cross_product

# Dot product of two vectors
dot_product = lambda vec_1, vec_2: sum(map(lambda x, y: x * y, vec_1, vec_2))

# Vector magnitude
def magnitude(vec):
return (vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2) ** 0.5

# Euler axis - angle
rot_rad = acos(dot_product(vec_global, vec_normal_in) / magnitude(vec_global) * magnitude(vec_normal_in))
# Euler axis - Axis of rotation
axis = cross_product(vec_global, vec_normal_in)
# Unit vector normal to Surface and Global X axis vector
axis_unit_vec = [axis[0] / magnitude(axis), axis[1] / magnitude(axis), axis[2] / magnitude(axis)]

def get_euler_angles_from_axis(euler_axis, rot_angle):
# Euler axis–angle ↔ quaternions
q1 = euler_axis[0] * sin(rot_angle / 2)
q2 = euler_axis[1] * sin(rot_angle / 2)
q3 = euler_axis[2] * sin(rot_angle / 2)
q4 = cos(rot_angle / 2)

# Quaternion → Euler angles (z-y′-x″ intrinsic)
roll = atan2((2 * q4 * q1 + 2 * q2 * q3), (1 - 2 * (q1 ** 2 + q2 ** 2)))
pitch = asin(2 * (q4 * q2 - q1 * q3))
yaw = atan2((2 * q4 * q3 + 2 * q1 * q2), (1 - 2 * (q2 ** 2 + q3 ** 2)))

# Euler angles (degrees)
roll_deg = roll * 180 / pi
pitch_deg = pitch * 180 / pi
yaw_deg = yaw * 180 / pi

return yaw_deg, pitch_deg, roll_deg

yaw_deg, pitch_deg, roll_deg = get_euler_angles_from_axis(axis_unit_vec, rot_rad)

# New coordinate system creation and orientation change to match surface normal vector
selection = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.MeshNodes)
selection.Ids = [mesh_node.Id]
new_cs.OriginLocation = selection
new_cs.RotateZ()
new_cs.RotateY()
new_cs.RotateX()
new_cs.SetTransformationValue(1, yaw_deg)
new_cs.SetTransformationValue(2, pitch_deg)
new_cs.SetTransformationValue(3, roll_deg)

# Creation of new path based on oriented coordinate system & X-axis intersection method
new_path.PathType = PathScopingType.Ray
new_path.PathCoordinateSystem = new_cs
new_path.NumberOfSamplingPoints = 100

# Add Stress object
vms.ScopingMethod = GeometryDefineByType.Path
vms.Location = new_path
vms.EvaluateAllResults()

vms.Activate()

paneTabular=ExtAPI.UserInterface.GetPane(MechanicalPanelEnum.TabularData)
control = paneTabular.ControlUnknown
col = control.ColumnsCount
row = control.RowsCount
abcissa = [float(control.cell(row_nr, 2).Text) for row_nr in range(2, row + 1)]
ordinate = [float(control.cell(row_nr, 3).Text) for row_nr in range(2, row + 1)]

def linear_interpolation(p1, p2, x):
slope = (x - p1[0]) / (p2[0] - p1[0])
p_int = p1[1] + (p2[1] - p1[1]) * slope
return p_int

at_x = 0.6
temp_val = filter(lambda x: (x > at_x), abcissa)[0]
index = abcissa.IndexOf(temp_val)

p1 = abcissa[index - 1], ordinate[index - 1]
p2 = abcissa[index], ordinate[index]
stress_at_x = linear_interpolation(p1, p2, at_x)

```

PyDPF

```"""

Stress gradient normal to a defined node.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This example shows how to plot a stress gradient normal to a selected node.
A path is created of a defined length.

"""

###############################################################################
# First, import the DPF-Core module as ``dpf_core`` and import the
# included examples file and ``DpfPlotter``
#
import matplotlib.pyplot as plt
from ansys.dpf import core as dpf
from ansys.dpf.core import operators as ops
from ansys.dpf.core.plotter import DpfPlotter
from ansys.dpf.core import examples

###############################################################################
# Next, open an example and print out the ``model`` object.  The
# :class:`Model <ansys.dpf.core.model.Model> class helps to organize access
# methods for the result by keeping track of the operators and data sources
# used by the result
# file.
#
# Printing the model displays:
#
# - Analysis type
# - Available results
# - Size of the mesh
# - Number of results
# - Unit
#
model = dpf.Model(examples.hemisphere)
print(model)
###############################################################################
# Define the `node_id` normal to which a stress gradient should be plotted.
#
node_id = 1928
###############################################################################
# The following command prints the mesh unit
#
print("Unit: %s" % unit)
###############################################################################
# `depth` defines the path length / depth to which the path will penetrate.
# While defining `depth` make sure you use the correct mesh unit.
# `delta` defines distance between consecutive points on the path.
depth = 10  # in mm
delta = 0.1  # in mm
###############################################################################
# Get the meshed region
#
###############################################################################
# Get Equivalent stress fields container.
#
stress_fc = model.results.stress().eqv().outputs.fields_container()
###############################################################################
# Define Nodal scoping.
# Make sure to define ``"Nodal"`` as the requested location, important for the
# `normals` operator.
#
nodal_scoping = dpf.Scoping(location="Nodal")
nodal_scoping.ids = [node_id]
###############################################################################
# Get Skin Mesh because `normals` operator requires Shells as input.
#
skin_mesh = ops.mesh.skin(mesh=mesh)
skin_meshed_region = skin_mesh.outputs.mesh.get_data()
###############################################################################
# Get normal at a node using `normals` operator.
#
normal = ops.geo.normals()
normal.inputs.mesh.connect(skin_meshed_region)
normal.inputs.mesh_scoping.connect(nodal_scoping)
normal_vec_out_field = normal.outputs.field.get_data()
###############################################################################
# Normal vector is along the surface normal. We need to invert the vector
# using `math.scale` operator inwards in the geometry, to get the path
# direction.
#
normal_vec_in_field = ops.math.scale(field=normal_vec_out_field,
ponderation=-1.0)
normal_vec_in = normal_vec_in_field.outputs.field.get_data().data[0]
###############################################################################
# Get Nodal coordinates, they serve as the first point on the line.
#
node = mesh.nodes.node_by_id(node_id)
line_fp = node.coordinates
###############################################################################
# Create 3D line equation.
#
fx = lambda t: line_fp[0] + normal_vec_in[0] * t
fy = lambda t: line_fp[1] + normal_vec_in[1] * t
fz = lambda t: line_fp[2] + normal_vec_in[2] * t
###############################################################################
# Create coordinates using 3D line equation-
#
coordinates = [[fx(t * delta), fy(t * delta), fz(t * delta)] for t in
range(int(depth / delta))]
flat_coordinates = [entry for data in coordinates for entry in data]
###############################################################################
# Create Field for coordinates of the path.
#
field_coord = dpf.fields_factory.create_3d_vector_field(len(coordinates))
field_coord.data = flat_coordinates
field_coord.scoping.ids = list(range(1, len(coordinates) + 1))
###############################################################################
# Let's now map results on the path.
mapping_operator = ops.mapping.on_coordinates(
fields_container=stress_fc,
coordinates=field_coord,
create_support=True,
mesh=mesh)
fields_mapped = mapping_operator.outputs.fields_container()
###############################################################################
# Here, we request the mapped field data and its mesh
field_m = fields_mapped[0]
mesh_m = field_m.meshed_region
###############################################################################
# Create stress vs length chart.
#
x_initial = 0.0
length = [x_initial + delta * index for index in range(len(field_m.data))]
plt.plot(length, field_m.data, "r")
plt.xlabel("Length (%s)" % mesh.unit)
plt.ylabel("Stress (%s)" % field_m.unit)
plt.show()
###############################################################################
# To create a plot we need to add both the meshes
# `mesh_m` - mapped mesh
# `mesh` - original mesh
pl = DpfPlotter()
color="w", opacity=0.3)
pl.show_figure(show_axes=True, cpos=[
(62.687, 50.119, 67.247),
(5.135, 6.458, -0.355),
(-0.286, 0.897, -0.336)])
```

• Member, Employee Posts: 16
Options

Ayush,

Very useful script, thanks for posting.

Could you show where this method is documented:

stress_fc = model.results.stress().eqv().outputs.fields_container()

If I look at the API I don't see how to do this.

• Member, Moderator, Employee Posts: 414
Options

Jim,

Not sure if this method is documented separately, I found it the generic example: