Extract stress at specific location

[robbievanleeuwen]

Discussed in #174

^{Originally posted by scote89 January 18, 2022}
Could we please have a function or example code on how to output the stress at a specific coordinate? I do not know how to get position values linked with the stress result to then interpolate.

[connorferster]

Could this play off of the method that generated the Mohr's circle? Does that not get the stress at a point?

[Czarified]

The original discussion thread was #68. I believe the issue is that the Mohr's Circle function doesn't return values, just the graph. If you want to pull the values and pass them on for anything else, you need to write your own function.

But I think you're right, we could take colin's method and break it up into 2 different methods.

[robbievanleeuwen]

There is also a bit of work in performing the correct interpolation based on the Tri6 shape functions. Currently the Mohr's circle method performs a linear interpolation. Once implemented, the Mohr's circle method will call this new get stress method to obtain the stress at a point.

[ccaprani]

@robbievanleeuwen I'll give this a go in the next few days. I'll use the shape function interpolations and expose a public API for it, and then update the Mohr's Ciricles accordingly. It would be a good addition indeed!

[normanrichardson]

@robbievanleeuwen and @ccaprani I had a look at this problem as I wanted to look at my FEM theory again. I had some issues getting the mapping (for a point in (x,y) to a point in (eta, xi, zeta)) right, but I think it is correct. I think I've made some progress that may be useful: a method in Tri6 that matches Tri6.element_stress. As is, it is a bit clunky to use, but it seems to work.

To make it easier to use I think some architectural decisions are needed. In particular, how to integrate this into StressPost.

Would this be a valuable contribution at this point?

One technical question I have is if it is correct to use the shape functions to interpolate the element's nodal stresses? As opposed to interpolating the average nodal stresses.

[ccaprani]

@normanrichardson those stress profiles are really nice to see! Slicing a section would be a super feature I think.

One technical question I have is if it is correct to use the shape functions to interpolate the element's nodal stresses? As opposed to interpolating the average nodal stresses.

Yes, it's tricky. The actual FE solution solves for the stresses at the Gauss Points. The shape functions are then used to extrapolate these results to the nodal points (https://sectionproperties.readthedocs.io/en/latest/rst/theory.html#extrapolation-to-nodes). So to recover stresses from the nodal points for an arbitrary location in the element there are two choices really:

• Work off the nodal stress, but make sure to use the shape functions as defined for the gauss points. What this means is that for the node near a particular gauss point, the relevant shape function will have an ordinate bigger than 1.0
• Convert the nodal stress back to the gauss point stress (or retrieve the GP results if they are already stored), and then use the shape functions without any need to start from an extrapolated value (i.e. the relevant shape function will have an ordinate of 1.0 at that GP). The stress at an arbitrary location in the element will then be found in the natural coordinates nicely.

Both approaches should give the same result.

At a more minor point, I don't think you need to recalculate the stresses within the elements from the external forces - we already have those from https://github.com/normanrichardson/section-properties/blob/deca66eea2008d1a7d7a25fb1fb7f5cbba05f3b2/sectionproperties/analysis/fea.py#L391

I meant to get to this sooner, but other things popped up. I'm happy to still do it, or if you want to do it, that's fine too; happy to give any input!

[normanrichardson]

@ccaprani great thanks for this, I think I have done something closer to the second approach.

The starting point was looking at Extrapolation to Nodes. Here is constructed as the mapping of the stress values from the nodes to the gauss points. Then it follows that maps the stress values from the gauss point to the node. Separately I confirmed this by constructing , inverted it, and verifying this with fea.extrapolate_to_nodes.

Looking at I saw row 1 of , which is , maps the stress from the nodes to gauss point 1, then at an arbitrary must map the stress over the domain. So if I want a stress at an arbitrary , then a similar process is:

So what I have put together so far is a method in Tri6 that does this. It assumes the provided point is in the triangle, finds the local coordinates , computes the shape function at the local coordinates , and multiplies this by the stress at the nodes (which I got from Tri6.element_stress).

@ccaprani Does this seem close to the second approach?

[ccaprani]

@normanrichardson lovely - you're right; once you have the nodal stresses, and your using the natural coordinates (where nodes are at (e.g.) (1,0,0), then you can just use the natural coordinate of the point of interest in the shape functions. No need to go back to the GP stresses at all!

FYI - which took me embarrassingly long to figure out - is that hardcoded H_inv comes from inverting the shape function matrix as evaluated at the gauss point coordinates.

In terms of the sectionproperties API, we already have the nodal stresses in a StressPost object and so a function like (sigma_zz, tau_xz, tau_yz) = stress_post.stress_at_point(x,y) is all that's needed I think (@robbievanleeuwen ?). It's low-level, but other things could build from it easily. For example:

• The Mohr's Circle plot interpolation code L4213-4221
• A stress_post.stress_along_slice() functionality, like you've indicated above.

And just because they are so lovely, here is a plot of the shape functions

import numpy as np
import matplotlib.pyplot as plt

def shape_funcs(xi, eta, zeta):
    return np.array(
        [
            eta * (2 * eta - 1),
            xi * (2 * xi - 1),
            zeta * (2 * zeta - 1),
            4 * eta * xi,
            4 * xi * zeta,
            4 * eta * zeta,
        ]
    )

# Cartesian coords
npts = 51
xi, eta = np.meshgrid(np.linspace(0, 1, npts), np.linspace(0, 1, npts))
xi = xi.flatten()
eta = eta.flatten()
zeta = 1 - xi - eta

indx = np.where(zeta < 0)
xi = np.delete(xi, indx)
eta = np.delete(eta, indx)
zeta = np.delete(zeta, indx)

N = shape_funcs(xi, eta, zeta)

fig = plt.figure(figsize=(8, 4.5))
for i, n in enumerate(N):
    ax = fig.add_subplot(2, 3, i + 1, projection="3d")
    ax.view_init(azim=-120, elev=32)
    ax.plot_trisurf(xi, eta, n, cmap=cm.jet, linewidth=0.2)
    ax.set_xlabel(r"$\xi$")
    ax.set_ylabel(r"$\eta$")
    ax.set_title(fr"$N_{i+1}(\xi,\eta)$")

fig.tight_layout()

[robbievanleeuwen]

Hi @ccaprani & @normanrichardson,

Yes I think this method is valid. I does seem like a bit of extra work to extrapolate to the node then interpolate back to an abritrary point but I like the simplicity of maintaining the use of the H matrix.

In terms of the implementation, I agree with @ccaprani, a decision just has to be made about whether we give the user all computed stresses, as in .get_stress(), or if we just give the 3 resultant components. The former will involve computing some of the combined stresses, currently done in section.py.

Looking forward to this PR :)

[normanrichardson]

@ccaprani thanks for this. It illustrates something that took me a while to figure out, that the elements local origin (0,0,1) is associated with node . That is, since:

then at the local origin :

as and .

This was important to the global to local mapping formulation. I initially tried to invert the above expression but gave up on that. Then I settled on using a kind of affine transformation:

is used for the above reason, and M is solved for each element.

I convinced myself that this was equivalent by: taking a global point and applying the affine transformation (to make its local point ) followed by produces .

This approach is implemented here and here. Let me know if anyone has some insight into this, or how it could be done better.

@robbievanleeuwen I have to leave this for another week or two and it seems like a good point to leave it for now. Should I submit this as a PR (or draft PR) for anyone who has time and sees an opportunity in building out the StressPost components?

[ccaprani]

@normanrichardson I should have looked at this a bit closer previously, because it's not so straightforward after all - my apologies! The great benefit of isoparametric elements is that the same shape functions are used for both mapping from physical to natural element coordinates, and for interpolating strain/stress fields within the element. This means that a quadratic triangle element in natural coordinates can have curved sides in physical coordinates. The mapping is almost always done from the natural to physical, since the element is developed in the natural coordinate system.

Our goal here is to map coordinates from the physical to the natural system, which is an inverse parametric mapping. For quadratic elements this is then a nonlinear transformation, and implementations of this use a Newton iteration scheme to solve the equations (e.g. Section 3 of this paper). A basic example is shown in Slides 5 & 6 here.

The affine transformation you used is linear, and so equivalent to the linear interpolation I used in physical coordinates for the Mohr's Circle L4213-4221.

With deference to @robbievanleeuwen of course, I'm not convinced that the increase in accuracy obtained is worth the significant development/computational effort, especially when the mesh is in anyway dense, as would typically be the case for a section properties calculation. So my suggestion from here would be to refactor the code so that the physical to natural coordinates is abstracted to a single function. For now this could be just left at a linear interpolation, and we can then later implement the nonlinear solver.