Calculating the traction vector field for a surface

[cumberworth]

I would like to be able to calculate and visualize the traction vector field for a given surface in a problem with hyperelasticity and Dirichlet boundary conditions. Below is sample code for an example problem in which a displacement is enforced on one side of a cube while the other is held in place.

import dolfin as dlf
import matplotlib.pyplot as plt
import numpy as np


def main():

    # Optimization options for the form compiler
    dlf.parameters['form_compiler']['cpp_optimize'] = True
    dlf.parameters['form_compiler']['representation'] = 'uflacs'

    # Define system
    E = 10.0
    nu = 0.3
    mu = dlf.Constant(E/(2*(1 + nu)))
    lmbda = dlf.Constant(E*nu/((1 + nu)*(1 - 2*nu)))

    L = 0.1
    W = 0.1
    mesh = dlf.BoxMesh(dlf.Point(0, 0, 0), dlf.Point(L, W, W), 2, 2, 2)
    V = dlf.VectorFunctionSpace(mesh, 'CG', 1)

    # Define boundaries
    facets = dlf.MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
    facets.set_all(0)
    left = dlf.CompiledSubDomain('near(x[0], side) && on_boundary', side=0.0)
    left.mark(facets, 1)
    right = dlf.CompiledSubDomain('near(x[0], side) && on_boundary', side=L)
    right.mark(facets, 2)

    # Create boundary conditions
    cl = dlf.Constant((0.0, 0.0, 0.0))
    cr = dlf.Expression(('gamma', '0.0', '0.0'), gamma=0.1, degree=1)
    bcl = dlf.DirichletBC(V, cl, left)
    bcr = dlf.DirichletBC(V, cr, right)
    bcs = [bcl, bcr]


    # Solve for displacement vector field
    u = dlf.Function(V)
    du = dlf.TrialFunction(V)
    v = dlf.TestFunction(V)
    dim = len(u)
    F = dlf.Identity(dim) + dlf.grad(u)
    body_force = dlf.Constant((0.0, 0.0, 0.0))
    traction_force = dlf.Constant((0.0, 0.0, 0.0))
    E = green_lagrange_strain(u)
    energy_density = lmbda/2*dlf.tr(E)**2 + mu*dlf.tr(dlf.dot(E, E))
    pot = (energy_density*dlf.dx - dlf.dot(body_force, u)*dlf.dx -
           dlf.dot(traction_force, u)*dlf.ds)
    F = dlf.derivative(pot, u, v)
    J = dlf.derivative(F, u, du)
    dlf.solve(F == 0, u, bcs, J=J)

    # Save solution in VTK format
    out_file = dlf.File('outs/saint_venant-kirchhoff_displacement.pvd')
    out_file << u

    # Calculate second Piola-Kirchhoff pseudo traction vector field
    N = dlf.FacetNormal(mesh)
    T = dlf.dot(second_piola_kirchhoff_stress(u, lmbda, mu), N)


def green_lagrange_strain(u):
    linear_part = dlf.grad(u) + dlf.grad(u).T
    nonlinear_part = dlf.dot(dlf.grad(u), dlf.grad(u).T)

    return 1/2*(linear_part + nonlinear_part)


def second_piola_kirchhoff_stress(u, lmbda, mu):
    E = green_lagrange_strain(u)

    return 2*mu*E + lmbda*dlf.tr(E)*dlf.Identity(len(u))


if __name__ == '__main__':
    main()

I seem to be able to integrate over the surface to calculate the mean normal force vector,

Tnn_c = dlf.dot(T, N)
Tnn = Tnn_c*N
ds = dlf.ds(domain=mesh, subdomain_data=facets)
left_normal_force = dlf.assemble(Tnn_c*ds(1))
right_normal_force = dlf.assemble(Tnn_c*ds(2))

However, attempts to output the traction vector so that it can be visualized with Paraview fail. If I naively try to project the traction vector field,

W = dlf.VectorFunctionSpace(mesh, 'DG', 1)
dlf.project(T, W)

I end up with the error,

ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal.

If I instead try to manually project and integrate over surfaces only,

W = dlf.VectorFunctionSpace(mesh, 'DG', 1)
w = dlf.TrialFunction(W)
v = dlf.TestFunction(W)
a = dlf.dot(w, v)*dlf.ds
L = dlf.dot(T, v)*dlf.ds
w = dlf.Function(W)
dlf.solve(a == L, w)
out_file = dlf.File('outs/saint_venant-kirchhoff_traction.pvd')
out_file << w

I get the following error

*** Error:   Unable to compute matrix factorisation.
*** Reason:  The provided data did not satisfy the prerequisites.
*** Where:   This error was encountered inside EigenLUSolver.cpp.

I wondered whether I should be trying to use the boundary mesh, so I ran

bmesh = dlf.cpp.mesh.BoundaryMesh(mesh, 'exterior')
N = dlf.FacetNormal(bmesh)
T = dlf.dot(second_piola_kirchhoff_stress(u, lmbda, mu), N)

but ended with the error

TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'

The closest posts on the forums that I could find are here, which involves a different variational form but seems to be specific to a Neumann boundary condition, and here, which has much more going on and many things that are not clear to me and don’t seem to be documented anywhere.

I am using dolfin version 2019.1.0, which I installed directly (i.e. I am not running it with docker).

Any input would be greatly appreciated.

[dokken]

You should be able to do something along the lines of: Iterate connected mesh entities in parallel

[cumberworth]

Thanks, based on that answer I was able to construct something that seems to work. I have included it below for anyone else with a similar issue.

# Calculate second Piola-Kirchhoff pseudo traction vector field
N = dlf.FacetNormal(mesh)
T = dlf.dot(second_piola_kirchhoff_stress(u, lmbda, mu), N)
dleft = dlf.ds(domain=mesh, subdomain_data=facets, subdomain_id=1)
dright = dlf.ds(domain=mesh, subdomain_data=facets, subdomain_id=2)
V = dlf.VectorFunctionSpace(mesh, 'CG', 1)
u = dlf.Function(V, name='traction-vector')
du = dlf.TrialFunction(V)
v = dlf.TestFunction(V)
L = dlf.assemble(dlf.dot(T, v)*(dleft + dright))
A = dlf.assemble(dlf.dot(du, v)*(dleft + dright), keep_diagonal=True)
A.ident_zeros()
dlf.solve(A, u.vector(), L)
out_file = dlf.File('outs/forum-post_traction-vector.pvd')
out_file.write(u)

I have to admit its not totally clear to me why it only works with passing the keep_diagonal=True option and the A.ident_zeros(). I am also not sure if it is better to use CG or DG (I have seen in other forum posts that calculate stresses, DG seems to be used, although its not so clear to me why). I am also not clear on the reasoning behind using degree 2 instead of degree 1 as was done in the answer linked to, so here I just stuck with degree 1. In testing all four, the changes in the result are minor for this example.

[dokken]

So, as you only integrate over the exterior facets, there will be several roes in your matrix filled with only zeros, making it invertible. ident_zeros identifies these rows and put a 1 on the diagonal to make it invertible.

Let’s assume you have used a CG-2 function-space to calculate deformation. Then, to get stresses, you differentiate this function. Differentiating a CG-2 function leads to a DG-1 function.

[cumberworth]

Okay, so if I want to project a quantity that involves the derivative of the deformation, I should consider the space I used to estimate the deformation.

But here, when I go back and calculate the deformation in CG-2 space and then calculate the stress vector in DG-1 space, I end up with some facets missing their vectors. Specifically, the vectors around two edges of the right surface (the surface which is being pulled by the boundary condition) are missing,

Compare this to when I instead use CG-1 space to calculate the stress vector, which includes these vectors,

Although even here there is an asymmetry present, as the vectors indicate a shear stress that should not exist. Perhaps this is to do with the way the mesh is constructed, instead of it being an artifact of the projection?

The missing vectors seem to stem from the face normal vectors. When I project the facet normal vectors onto DG-1, the same vectors are missing,

When I project onto CG-1 space, they are present.

So then it would seem better to project to CG-1 space?

[dokken]

Note that to visualize DG functions, you need to use XDMFFile.write_checkpoint, see: Loading xdmf data back in

[cumberworth]

Ah, okay, thanks. With that all the vectors are visible. However, now at least some of the vertices appear to have two slightly different traction vectors. This is not visible when only projecting the normal vectors. I use the following to output,

out_file = dlf.XDMFFile('outs/forum-post_traction-vector.xdmf')
out_file.write_checkpoint(u, 'f', encoding=dlf.XDMFFile.Encoding.ASCII)

which when viewed in Paraview looks like this

[dokken]

Of course there will be multiple tractions at the vertices. Each cell has a separate dof at this point (as it is a DG space) and the normal at this vertex is not uniquely defined.

[cumberworth]

Oh, yes that makes sense. Thanks for all your help!