w = Function(W)
solve(a == L, w, bcs)
u, p = w.split()
def symgrad(u, i):
return sym(grad(u))
K = np.zeros((3,3))
for i in range(3):
for j in range(3):
K[i,j] = assemble(2*mu*inner(symgrad(u, i), symgrad(u, j))*dx)
print("K_ij =", K[i,j])
This is the minimal code I think, what i would to do is to use the output u from the solver and then integrate over dx the formula above. You think this way is appropriate for the forum?