Feat stress elm local

[normanrichardson]

This is draft for some of basic elements of issue #175. It gives the ability for the Tri6 element to determine its stress at a point within its domain.

Known limitations:

1. The mappings in Tri6.local_coord that take a global coordinate and maps it to a local coordinate (or physical coordinate to a natural coordinate) are based on alphine transformations and are not general enough for isoparametric triangular elements. The mappings are sufficient for "straight line" Tri6 elements. I believe the underlying mesher triangle is limited to "straight line" Tri6's.

In it's current state this can be used as follows:

import sectionproperties.pre.library.steel_sections as steel_sections
from sectionproperties.analysis.section import Section
import sectionproperties.analysis.fea as fea
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

# Create a 150x100x6 RHS on its side
geometry = steel_sections.rectangular_hollow_section(d=100, b=150, t=6, r_out=15, n_r=8)

# Create a mesh and section object. For the mesh, use a maximum area of 2
geometry.create_mesh(mesh_sizes=[100])
section = Section(geometry)
section.plot_mesh()

# Perform a geometry and warping analysis
section.calculate_geometric_properties()
section.calculate_warping_properties()

# Section properties needed for local stress
sect_prop = {
    'ea':  section.section_props.ea,
    'cx':  section.section_props.cx,
    'cy':  section.section_props.cy,
    'ixx':  section.section_props.ixx_c,
    'iyy':  section.section_props.iyy_c,
    'ixy':  section.section_props.ixy_c,
    'i11':  section.section_props.i11_c,
    'i22':  section.section_props.i22_c,
    'phi':  section.section_props.phi,
    'j':  section.section_props.j,
    'Delta_s':  section.section_props.Delta_s,
    'nu':  section.section_props.nu_eff,
}

# Perform a stress analysis with Mx = 5 kN.m; Vx = 10 kN and Mzz = 3 kN.m
state = {
    'N': 0,
    'Mxx': 5e6,
    'Myy': 0,
    'M11': 0,
    'M22': 0,
    'Mzz': 3e6,
    'Vx': 10e3,
    'Vy': 0
}
case1 = section.calculate_stress(**state)

# set up plots
fig = plt.figure()

gs = GridSpec(2, 2, figure=fig)
ax1 = fig.add_subplot(gs[0, :])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[1, 1])

# Plot the bending stress for case1
case1.plot_stress_v_zx(pausr=False, render=False, ax=ax1)

def getElementStess(p):

    # find the Tri6 element containing the point
    tri_id = None
    for (idx , tri) in enumerate(section.elements):
        if tri.point_within_element(p):
            tri_id = idx
    tri = section.elements[tri_id]

    # get the stess at the point
    sigs = tri.local_element_stress(
        p=p,
        **state,
        **sect_prop,
        omega=section.section_props.omega[tri.node_ids],
        psi_shear=section.section_props.psi_shear[tri.node_ids],
        phi_shear=section.section_props.phi_shear[tri.node_ids]
        )

    return sigs

# get and plot the stress at points
s = []
x = []
y = []
for ii in np.arange(5.015,148.99,4):
    # global point of interest
    p = np.array([ii,4.5])
    y.append(p[1])
    x.append(p[0])
    s.append(getElementStess(p)[7])

ax1.plot(x,y, c="green", marker='o', alpha=.6)
ax2.plot(x,s, c="green", marker='o')

s = []
x = []
y = []
for ii in np.arange(5.015,94.99,3):
    # global point of interest
    p = np.array([5, ii])
    y.append(p[1])
    x.append(p[0])
    s.append(getElementStess(p)[7])
ax1.plot(x,y, c="blue", marker='o',alpha=.6)
ax3.plot(y,s, c="blue", marker='o')

plt.show()

[robbievanleeuwen]

This looks great @normanrichardson! I think we can get it to work based off what you've done here, and I'm happy with the linear interpolation results!

Design-wise, I think we can store the section actions in the StressPost class and this could be passed to a Section.get_stress_at_point(pt, stress_post) method. We could also look at the stress_post (or similarly named) variable being either a StressPost object or a dictionary of actions if the user doesn't want to run a full stress analysis on the section. I think this is maybe a good idea, but open to your thoughts?

All the section properties could be retrieved from the SectionProperties class attribute and then we should all have all the relevant information to pass to Tri6.local_element_stress(). I think it makes sense to combine the stresses in Tri6.local_element_stress() to only output (sigma_zz, tau_xz, tau_yz) at previously discussed.

Moving forward, do you have the time/enough information to finish this off? If not, I can have a look into it over the next month or so. Thanks again for all your work on this!

[robbievanleeuwen]

Also love the line plots you've done, seems like it should be easy to build a plot_stress_along_line() function or similar!

[Czarified]

Thanks so much @normanrichardson ! Great work!

[normanrichardson]

@Czarified no worries, thanks for the prodding.

@robbievanleeuwen, yeah I think I can finish this up.

"store the section actions in the StressPost class"

In section.calculate_stress when the StressPost is created I will adjust this to take the in the forces (N, Vx, Vy, Mxx, Myy, M11, M22, Mzz) and store these in StressPost. Then basically just copy over the function getElementStess(p) defined above to StressPost. Thoughts beyond may become more clear as I go.

One thought is to change stress_post to calculate and hold the element stresses, and not the averaged nodal stresses, or maybe hold both. But I would hesitate to go that far for now.

"only output (sigma_zz, tau_xz, tau_yz)"

It has been a while, so I'm a little rusty on the details here. Are these StressResult.sig_zz, StressResult.sig_zx , and StressResult.sig_zy? And can find those combinations here?

[normanrichardson]

@robbievanleeuwen I think I miss understood.

Would Section.get_stress_at_point be a compliment to Section.calculate_stress?

The way I was describing it would work:

# the usual above
state = {
    'N': 0,
    'Mxx': 5e6,
    'Myy': 0,
    'M11': 0,
    'M22': 0,
    'Mzz': 3e6,
    'Vx': 10e3,
    'Vy': 0
}
case1 = section.calculate_stress(**state)
p = np.array([5, ii])
stress_p = case1.get_stress_at_point( p)

But what I'm now thinking you're saying:

# the usual above
state = {
    'N': 0,
    'Mxx': 5e6,
    'Myy': 0,
    'M11': 0,
    'M22': 0,
    'Mzz': 3e6,
    'Vx': 10e3,
    'Vy': 0
}
case1 = section.calculate_stress(**state)
stress_p = section.get_stress_at_point( p, **state)

Where the type of case1 and stress_p could differ.

Is this closer to what you are thinking?

[robbievanleeuwen]

Yes I was thinking the second option because if a user was only interested in the stress at certain points they wouldn't have to run the entire analysis using Section.calculate_stress(). If I understand the implementation correctly, we don't need a StressPost object at all, just access to the calculated section properties, and we are only performing a calculation on one element. I was thinking the Section method could look like:

def get_stress_at_point(
    self,
    pt: List[float],
    actions: Dict[str, float],
) -> Tuple[float]:
    """description...

    e.g. pt = [0.5, 1.5]
    e.g actions = state (from your comment above)

    :return: Resultant normal and shear stresses (sigma_zz, tau_xz, tau_yz)
    """

    # do stuff

Disregard the comment about saving the actions to the StressPost object, I think this is a lot cleaner and should be done independently of StressPost.

It has been a while, so I'm a little rusty on the details here. Are these StressResult.sig_zz, StressResult.sig_zx , and StressResult.sig_zy? And can find those combinations here?

Correct, these are just the resultant stresses, calculation as you've highlighted!

[normanrichardson]

Potential issue:
I'm not sure I like how points on the element boundaries are handled. Right now it takes the values from the last element it finds, here. If I remember correctly, the stress over element boundaries is not continuous. So some sort of averaging may make sense.

Other potential improvements:
The element search is not great. I believe there are ways to speed this up. Initially, I would try using shapely to generate the mesh polygons and perform STR-tree search. This should help if asking for more than one point (it may also help for just one point). Should I add a function get_stress_at_points that does this?

[robbievanleeuwen]

Correct, pretty sure the stresses are not continuous at element boundaries. This is a common theme in FEA with extrapolating stresses from gauss points. If you could average the boundary result based on the two adjacent elements I think that would be suitable?

Regarding the element search, if you think it is worth adding this feature I'm all for it 🙌

[Czarified]

From my experience, and referencing a "Practical FEA" textbook, element-nodal bilinear extrapolation with nodally averaged boundary results is common, and with a slight conservatism. I could look up the reference, if you're curious.

[normanrichardson]

@robbievanleeuwen let me know your thoughts. Thanks

[normanrichardson]

I realized that I need to update the docstring to include the return type of None, for the case when a point is outside the section. Will do this this evening

[robbievanleeuwen]

This is great, thanks @normanrichardson! A few small minor spelling/docs items but the rest of it looks good to me!

[Czarified]

Thanks so much, @normanrichardson ! Great work on this feature!

[robbievanleeuwen]

Agreed thanks again @normanrichardson! I think this is a really useful addition to sectionproperties! I'm hoping to overhaul the docs in the near future and will definitely be adding an example that uses this feature!

[normanrichardson]

No worries, any time. It was good to partially dust off the fea knowledge