calculate Von-Mises Stress element by element in a 3D Case

[Preetzl]

Hey kinnala,
in exp. 4 you already implemented the Von Mises Stress and viusalized it.
My plan was to calculate the von mises stress tetrahedron by tetrahedron and compare whether the yield stress is exceeded or not. I was hoping that I could bring it to work for this Code:
`import numpy as np
from skfem import *
from skfem.helpers import dot, grad, div
from skfem.models.poisson import *
from scipy.sparse import spmatrix
from scipy.sparse.linalg import LinearOperator
from skfem.models.elasticity import linear_elasticity, lame_parameters
import meshio
from skfem.visuals.matplotlib import draw
m = MeshTet()
m = m.load("C:/Users/AlexP/Desktop/BA/Meshes/Cube.vtk")

e1 = ElementTetP1() # ElementHex1() # ElementTetP1()
e = ElementVector(e1)
ib = Basis(m, e, MappingIsoparametric(m, e1), 3)

K = asm(linear_elasticity(*lame_parameters(3e9, 0.33)), ib) # young's modulus given in GPa, Poisson ratio for PLA

dofs = {
'left': ib.get_dofs(lambda x: np.isclose(x[0], x[0].min())),
'right': ib.get_dofs(lambda x: np.isclose(x[0], x[0].max())),
}

u = ib.zeros()

force applied to knots

Kraft = 5000000
F = np.zeros(u.shape)
left_dofs = dofs['left'].nodal['u^1']
F[left_dofs] = Kraft

fixed_dofs = np.hstack([dofs['right'].all()])
u[fixed_dofs] = 0

#set dofs and solve DGL
free_dofs = np.setdiff1d(np.arange(K.shape[0]), fixed_dofs)
K_free = K[free_dofs][:, free_dofs]
F_free = F[free_dofs]
u_free = solve(K_free, F_free)
u[free_dofs] = u_free

calculating upper shift only for the hard coded example

verschiebung_links = u[left_dofs]

implement shifts

u[left_dofs] = verschiebung_links

I = ib.complement_dofs(dofs)

u = solve(*condense(K, x=u, I=I))

print(u.max(), u.min())

sf = 2.0
m = m.translated(sf * u[ib.nodal_dofs])
draw(m)
#m.save("output_hex" + '.vtk', {"affected": u[ib.nodal_dofs][0]})`

Do I have to work with a faceted base to get the stress? Because at the moment I'm only using knots that don't have an area.

Would it be possible to iterate through the tetrahedrons one by one?

Thank you for your nice library.

Best Regards

Preetzl

[kinnala]

Here are some thoughts:

• What do you mean by a knot?
• It is possible to loop over elements but it is against the philosophy of scikit-fem...
• We try to avoid looping over elements, since it is computationally inefficient to do loops in Python, this is why numpy was invented
• It is perhaps a better idea to calculate von Mises stress in all elements in one numpy call and then find which of them will yield

So you can then just slightly modify example 4:

s = {}
dgb = ib.with_element(ElementTetP0())
up = ib.interpolate(u)

for i in [0, 1]:
    for j in [0, 1]:
        s[i, j] = dgb.project(C(sym_grad(up))[i, j])

s[2, 2] = poisson_ratio * (s[0, 0] + s[1, 1])

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))

Finally, the variable vonmises should have one value per element and you can check which of them yield.