Please note that you should not use CG
1 function spaces for Von Mises stresses, as the stresses are a DG 0 function (if the displacement is CG 1. I.e.
VDG = FunctionSpace(mesh, "DG", 0)
von_Mises = project(von_Mises, V)
You can compute the stresses on a subdomain by projecting on a subdomain, as follows:
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = sqrt(3./2*inner(s, s))
u, v = TrialFunction(VDG), TestFunction(VDG)
a = inner(u, v)*dx
L = inner(von_Mises, v)*dx(1)
stress = Function(VDG)
solve(a==L, stress)
would only compute the stresses on domain 1, and be zero everywhere else.
This can be optimized in a similar fashion to How to extract values in a surface of 3D mesh to another 2D mesh? - #4 by dokken
where you would use dx(1)
in the definition of a
, then assemble it into A
, and use ident_zeros
.
Another solution would be to use MeshViews