How to restrict solutions or stresses to subdomains?

I/O
1

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! :slightly_smiling_face:

2

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
3

Thanks, dokken, for your quick and helpful answer! You enabled me to implement your suggestion and it works like a charm. :+1:

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.

4

There are examples in the dolfin bitbucket repo: Bitbucket

5

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?

1 Like
6

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.