How to restrict solutions or stresses to subdomains?

I/O
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