 # How to restrict solutions or stresses to subdomains?

Hi,

I have a test setup for a simple linear elastic problem: I apply an external force and want to compute displacements and von Mises stresses. My setup is basically that of this Tutorial and consists of two bodies that are connected by a third one.

Departing from the tutorial example, each body has its own material properties. Consequently, the stress tensor sigma depends on the subdomain, since the Lamé parameters change between bodies.

I implemented different material properties by adding parameters to the definition of sigma:

``````def sigma(u, E, nu):
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
return lmbda * nabla_div(u) * Identity(d) + 2 * mu * epsilon(u)
``````

For context: Since my mesh already comes with marked subdomains, I use

``````mesh = Mesh('mesh_mwe.xml')
markers = MeshFunction('size_t', mesh, 'mesh_mwe_physical_region.xml') # three 3D bodies are marked 1, 2, 3
boundaries = MeshFunction('size_t', mesh, 'mesh_mwe_facet_region.xml') # two 2D boundaries are marked 4, 5
dx = Measure('dx', domain=mesh, subdomain_data=markers)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
``````

to load the mesh and update the measures dx and ds. The variational problem then becomes

``````u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
a = sum(inner(sigma(u, E[k], nu[k]), epsilon(v))*dx(k) for k in [1, 2, 3])
L = sum(dot(traction[k], v)*ds(k) for k in [4, 5])
``````

with predefined lists for E, nu, and traction (which in my case is simply constant on each boundary).

I can compute the solution and displacements, but I am stuck with von Mises stresses. The tutorial uses,

``````s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1) # Why is this line necessary?
von_Mises = project(von_Mises, V)
``````

but my sigma depends on the subdomain via the definition above.

How can I solve this?

• Can I restrict u to a subdomain (or subdomains) and compute stresses only for those subdomains? If so, how?

This is basically the same question as posted here (t/stress-for-one-subdomain-in-multiple-subdomain-problem/919), which remained unanswered. Also, I could not find any information on this topic in this Tutorial, which explains how to setup subdomains but not how to adapt the post-processing.

Any help is appreciated! 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

1 Like

Thanks, dokken, for your quick and helpful answer! You enabled me to implement your suggestion and it works like a charm. Out of curiosity: Could you provide a link to any documentation on MeshView, please? I could not find more than API reference but could not understand how to apply it without any examples.

There are examples in the dolfin bitbucket repo: Bitbucket

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)
gmsh.model.occ.fragment([cutcyl], [cube, box]) # makes mesh conformal
name2tag['cylinder'] = cyl

gmsh.model.occ.synchronize()

name2tag['clamped'] = (2, )
name2tag['forced'] = (2, )

grouptags = dict()
for name, (dim, tags) in name2tag.items():
if isinstance(tags, int):
tags = [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):

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?

Note that you can use a DG0 function for the Lame parameters, such that you only need one definition og sigma, see: How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

That is correct.