Von Mises stress at quadrature points

[gittsoy]

Hi,

in Example 4 we compute the von Mises stress and obtain 1D numpy.ndarray with size [number of elements * number of quadrature points per element].

Is there a way to convert it to 1D array with size [number of nodal dofs] collecting stress values at global nodal dofs?

[kinnala]

Let's say you use piecewise linear element for the displacement. Then its derivative (i.e. stress) is elementwise constant. So you have a different stress value in each element. Thus, in order to have a unique value at a node, you must perform some sort of projection.

You could try using skfem.utils.project to project vonmises1 from ElementTriDG(ElementTriP1()) to ElementTriP1() and see what happens. I can come up with an example later if you want.

[kinnala]

Let me add that another option is to first calculate the derivatives of the displacement at nodes and then multiply and add those accordingly to arrive at the von Mises stress. I'm not sure which one is a better approach.

[gittsoy]

Let's say you use piecewise linear element for the displacement. Then its derivative (i.e. stress) is elementwise constant. So you have a different stress value in each element.

So if we use P2 element, the displacement derivative will be an elementwise linear function. But as I understand the values at the coinciding nodes of the two elements may differ anyway.

You could try using skfem.utils.project to project vonmises1 from ElementTriDG(ElementTriP1()) to ElementTriP1() and see what happens. I can come up with an example later if you want.

I used skfem.utils.project and got the following code

from skfem import *

def calculate_misses(x, Lambda, mu, mesh, gb):    
    C = linear_stress(Lambda, mu)
    # post processing
    s = {}
    e_dg = ElementTriDG(ElementTriP1())

    for itr in range(2):
        for jtr in range(2):

            @BilinearForm
            def proj_cauchy(u, v, w):
                return C(sym_grad(u))[itr, jtr] * v

            @BilinearForm
            def mass(u, v, w):
                return u * v

            ib_dg = InteriorBasis(mesh, e_dg, intorder=4)

            s[itr, jtr] = solve(asm(mass, ib_dg),
                                asm(proj_cauchy, gb, ib_dg) @ x)

            s[2, 2] = 0 # General plane stress

    vonmises = np.sqrt(.5 * ((s[0, 0] - s[1, 1]) ** 2 +
                          (s[1, 1] - s[2, 2]) ** 2 +
                          (s[2, 2] - s[0, 0]) ** 2 +
                          6. * s[0, 1]**2))
    return vonmises

mesh = MeshTri()
basis = InteriorBasis(mesh, ElementVectorH1(ElementTriP2()), intorder=4)
clamped = basis.get_dofs(lambda x: x[0] == 0.0).all()
free = basis.complement_dofs(clamped)

@LinearForm
def loading(v, w):
    return  -1. * v.value[1]

E = 21.19e4
nu = 0.277
Lambda, mu = lame_parameters(E, nu)

# stiffness matrix
K = asm(linear_elasticity(Lambda, mu), basis)
f = asm(loading, basis)
u = solve(*condense(K, f, I=free)) 

m_deformed = MeshTri(mesh.p + 1.*u[basis.nodal_dofs], mesh.t)
vonmises = calculate_misses(u, Lambda, mu, m_deformed, basis) 
ax = plot(InteriorBasis(m_deformed, ElementTriDG(ElementTriP1())), vonmises, shading='gouraud') 
draw(m_deformed, ax=ax ,node_numbering=True)

vonmises_proj=project(vonmises, basis_to=InteriorBasis(m_deformed, ElementTriP1()), basis_from=InteriorBasis(m_deformed, ElementTriDG(ElementTriP1())))

print(np.reshape(vonmises, (-1, 3)))
print(vonmises_proj)

Seems to work.

[kinnala]

Great!