Thank you for the link!

It might be helpful to somebody (or even my future self) to have a **minimal working example** (mwe) on how to do the hole thing, so here is what I came up with. I do not claim that this is the best way to do it and am very open to constructive criticism.

The following mwe consists of three bodies, a box and a cube connected by a cylinder. The “bottom” side of the box is clamped and a constant traction acts on the “top” surface of the cube.

First create a dummy mesh with gmsh (I used version 4.6.0):

```
""
Generate a dummy mesh with three bodies: a box and a cube box connected by a cylinder.
"""
import subprocess
import gmsh
def main():
gmsh.initialize()
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
name2tag = dict()
tag = gmsh.model.occ.addBox(-1, -1, 0, 2, 2, -0.2)
box = (3, tag)
name2tag['box'] = box
tag = gmsh.model.occ.addBox(-0.5, -0.5, -2, 1, 1, -1)
cube = (3, tag)
name2tag['cube'] = cube
tag = gmsh.model.occ.addCylinder(0, 0, -0.1, 0, 0, -2, 0.25)
cyl = (3, tag)
cutcyl = gmsh.model.occ.cut([cyl], [cube, box], removeTool=False)[0][0]
gmsh.model.occ.fragment([cutcyl], [cube, box]) # makes mesh conformal
name2tag['cylinder'] = cyl
gmsh.model.occ.synchronize()
name2tag['clamped'] = (2, [27])
name2tag['forced'] = (2, [20])
# add physical groups
grouptags = dict()
for name, (dim, tags) in name2tag.items():
if isinstance(tags, int):
tags = [tags]
grouptag = gmsh.model.addPhysicalGroup(dim, tags)
gmsh.model.setPhysicalName(dim, grouptag, name)
grouptags[name] = (dim, grouptag)
print(name, grouptags[name])
gmsh.model.occ.synchronize()
# generate mesh (output only for PhysicalGroup elements)
elementsize = 0.10
gmsh.model.mesh.setSize(gmsh.model.getEntities(dim=0), elementsize)
gmsh.model.mesh.generate()
meshfilebasename = 'mesh_mwe'
gmsh.write('{}.msh'.format(meshfilebasename))
gmsh.fltk.run()
gmsh.finalize()
# convert msh to xml
command = 'dolfin-convert {f}.msh {f}.xml'.format(f=meshfilebasename)
subprocess.check_call(command.split())
if __name__ == '__main__':
main()
```

Then read the xml and setup and run the FEniCS calculation:

```
"""
Linear elastic problem (no body force -> r.h.s. is zero):
-div(sigma(u)) = 0
Restrict von Mises output to selected bodies (zero otherwise).
"""
from __future__ import print_function
from fenics import *
from ufl import nabla_div
import matplotlib.pyplot as plt
def main():
elementorder = 2
# Create mesh and define function space
mesh = Mesh('mesh_mwe.xml')
markers = MeshFunction('size_t', mesh, 'mesh_mwe_physical_region.xml')
boundaries = MeshFunction('size_t', mesh, 'mesh_mwe_facet_region.xml')
dx = Measure('dx', domain=mesh, subdomain_data=markers)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
V = VectorFunctionSpace(mesh, 'P', elementorder)
# Define boundary condition
bc = [DirichletBC(V, Constant((0, 0, 0)), boundaries, 4)] # clamp surface 4
# Define strain and stress
def epsilon(u):
return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
E = [0, 194000, 194000, 194000] # 0 -> dummy entry this remains unused
nu = [0, 0.3, 0.3, 0.3] # 0 -> dummy entry this remains unused
def Sigma(E, nu):
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
def sigma(u):
return lmbda * nabla_div(u) * Identity(d) + 2 * mu * epsilon(u)
return sigma
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
a = sum(inner(Sigma(E[k], nu[k])(u), epsilon(v))*dx(k) for k in [1, 2, 3])
L = dot(Constant((0, 0, -70)), v)*ds(5) # constant traction on surface 5
# Compute solution
u = Function(V)
solve(a == L, u, bc)
print('plot and write von Mises stresses for selected bodies only:')
tags = [1, 2] # tags of the bodies I am interested in
W = FunctionSpace(mesh, 'P', elementorder - 1) # derivatives in von_Mises reduce elementorder by one
v = TestFunction(W)
L = 0 # dummy 0, just to initialize L
for tag in tags:
sigma = Sigma(E[tag], nu[tag])
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) # u most NOT be a TrialFunction already
von_Mises = sqrt(3./2*inner(s, s))
L += inner(von_Mises, v)*dx(tag) # v must be a TestFunction here
u = TrialFunction(W)
a = inner(u, v)*dx
stress = Function(W)
solve(a==L, stress)
plot(stress)
plt.show()
File('elasticity/von_mises_tag{}.pvd'.format(''.join(map(str, tags)))) << stress
if __name__ == '__main__':
main()
```

Your example, dokker, worked perfectly well as long as the bodies didn’t have different material properties (and sigma didn’t depend on the body). It took me a while to realize I have to split this line in two to make it work in that case:

u, v = TrialFunction(VDG), TestFunction(VDG)

So thanks again for your help! I will look into MeshViews, too, though. Using submeshes seems like a good idea and maybe I can get rid of the zeros and get stresses *only* in the bodies I am interested in.

By the way, the element order of the FunctionSpace for von Mises stresses are reduced by one compared to that of the displacments, because of the derivative in sigma and epsilon. Did I get that right?