How can i derive interpolated value at nodal points?

[HoseongJeong]

Dear user of scikit-fem!

I am trying to perform solid FEA using scikit-FEM, and very thankfully, overcoming my challenges through this community!!

If it is possible, could you let me know how to drive average values (like stress) at the nodes?
For example, I have stress array as an input like bellows:

Stress=
[Stress at quadrature point 1, ..., stress at quadrature point n] <- element 1
[Stress at quadrature point 1, ..., stress at quadrature point n] <- element 2

And I would like to get the interpolated values at the nodal points, and average these to obtain nodal stress!

Thanks for all the comments and help!

[kinnala]

I will describe a conceptual idea for postprocessing of the stress and later some code on how to solve it.

Let $V_h$ be the finite element space for the displacement $\underline{u}_h$. This is a vector-valued finite element space,
for example, ElementVector(ElementTriP1()). If you have performed solid FEA you should already have the DOF vector corresponding to $\underline{u}_h$.

In the continuous setting, the stress $\underline{\underline{S}}$ is defined by $\underline{\underline{S}}(\underline{u}) = 2 \mu \underline{\underline{E}}(\underline{u}) + \lambda \mathrm{tr} \underline{\underline{E}}(\underline{u}) \underline{\underline{I}}$ where $\underline{\underline{E}}(\underline{u}) = \frac 12(\nabla \underline{u} + \nabla \underline{u}^T)$.

Let $W_h$ be the finite element space for the discrete stress $\underline{\underline{S}}_h$. This is a tensor-valued finite element space, for example, ElementVector(ElementVector(ElementTriP0())).

To obtain the discrete stress, we can solve the variational problem: find $\underline{\underline{S}}_h \in W_h$ such that

$$\int_{\Omega} \underline{\underline{S}}_h : \underline{\underline{T}} dx = \int_{\Omega} \underline{\underline{S}}(\underline{u}_h) : \underline{\underline{T}} dx$$

for every $\underline{\underline{T}} \in W_h$.

Discretizing the right hand side will involve interpolating the stress at the quadrature points. If you instead choose $W_h$ as ElementVector(ElementVector(ElementTriP1())) then an averaging of the stress is done automatically by the linear solve.

[HoseongJeong]

Thank you, kinala!!! I will try your recommendation!! And I deeply appreciate your code! It is amazingly useful!!

[kinnala]

Here is a simple solid FEA solver to get the displacement (scikit-fem==10.0.1):

from skfem import *
from skfem.helpers import *


m = MeshQuad().refined(3).with_defaults()
e = ElementQuad1()
basis = Basis(m, ElementVector(e))

# material params, some testing values
mu, lam = 1.0, 0.5


def E(u):
    # expression for the strain
    return 1 / 2 * (grad(u) + transpose(grad(u)))


def S(E):
    # expression for the stress
    return 2 * mu * E + lam * eye(trace(E), 2)


@BilinearForm
def work_int(u, v, _):
    # internal virtual work
    return ddot(S(E(u)), E(v))


K = work_int.assemble(basis)


# stretch the square
u = basis.zeros()
u[basis.get_dofs('right').all('u^1')] = 0.2
u = solve(*condense(K, 0 * u, x=u, D=basis.get_dofs({'left', 'right'})))
mdefo = m.translated(u[basis.nodal_dofs])

# visualize deformation
mdefo.draw().show()

[kinnala]

Here is the postprocessing part I described above:

# postprocessing

# basis for stress tensor
Sbasis = Basis(m, ElementVector(ElementVector(e)))

@BilinearForm
def lhs(s, t, _):
    return ddot(s, t)


@LinearForm
def rhs(t, w):
    return ddot(S(E(w['u'])), t)


A = lhs.assemble(Sbasis)
B = rhs.assemble(Sbasis, u=basis.interpolate(u))


# solve the postprocessing "projection problem"
# no boundary conditions required here
stress = solve(A, B)


# visualize different components of stress
pbasis = Basis(mdefo, e)
for itr in range(4):
    pbasis.plot(stress[Sbasis.nodal_dofs[itr]],
                shading='gouraud',
                nrefs=2).show()

[kinnala]

If you already have stress at the quadrature points and want to go from there, then you need to modify this part:

...
@LinearForm
def rhs(t, _):
    return ddot(your_stress, t)  # put here stress at quadrature points, in a numpy array of correct shape

# you can check the output of Sbasis.interpolate(Sbasis.zeros()).value
# to try to figure the correct shape of the numpy array
...

Edit: Fixed the advice