Hyperelasticity Demo (Dolfinx <-- Dolfin)

[bhaveshshrimali]

I was looking for a hyperelasticity example, but didn't find one. I am attaching a working demo (basically translating the old demo from dolfin), in case if someone comes looking up for it

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import BoxMesh, DirichletBC, Function, VectorFunctionSpace, cpp, Constant, NewtonSolver, NonlinearProblem
from dolfinx.cpp.mesh import CellType
from dolfinx.fem import (apply_lifting, assemble_matrix, assemble_vector,
                         locate_dofs_geometrical, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.log import set_log_level, LogLevel
from ufl import (
    Identity, SpatialCoordinate, TestFunction, TrialFunction, as_vector, dx, ds, grad, inner, ln, det,
    derivative)

set_log_level(LogLevel.INFO)
class Hyperelasticity(NonlinearProblem):
    def __init__(self, a, L, bcs, V):
        super().__init__()
        self.L = L
        self.a = a
        self._F, self._J = None, None
        self._bcs = bcs
        self._V = V

    def form(self, x):
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

    def F(self, x):
        if self._F is None:
            self._F = assemble_vector(self.L)

        else:
            with self._F.localForm() as f_local:
                f_local.set(0.0)
            self._F = assemble_vector(self._F, self.L)

        apply_lifting(self._F, [self.a], [self._bcs], [x], -1.0)
        self._F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        set_bc(self._F, self._bcs, x, -1.0)
        return self._F

    def J(self, x):
        if self._J is None:
            self._J = assemble_matrix(self.a, self._bcs)
        else:
            self._J.zeroEntries()
            self._J = assemble_matrix(self._J, self.a, self._bcs)
        self._J.assemble()
        return self._J

mesh = BoxMesh(
    MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
                     np.array([1.0, 1.0, 1.0])], [24, 16, 16],
    CellType.tetrahedron, cpp.mesh.GhostMode.none)

# boundary tolerance
btol = 1.e-6

def left(x):
    return x[0] < btol

def right(x):
    return np.abs(x[0]-1.) < btol

def rightBoundary(x):
    return np.stack((np.zeros(x.shape[1]),
        (0.5 + (x[1] - 0.5) * np.cos(np.pi/3) - (x[2] - 0.5) * np.sin(np.pi/3) - x[1])/2.,
        (0.5 + (x[1] - 0.5) * np.sin(np.pi/3) + (x[2] - 0.5) * np.cos(np.pi/3) - x[2])/2.))

x = SpatialCoordinate(mesh)
bodyForce = as_vector([0, -0.5, 0])
surfaceTraction = as_vector([0.1, 0, 0])
E, nu =  Constant(mesh, 10.0), Constant(mesh, 0.3)

mu = E/2./(1+nu)
lmbda = 2 * mu * nu / (1. - 2 * nu)

V = VectorFunctionSpace(mesh, ("Lagrange", 1))
u = Function(V, name="u")
v = TestFunction(V)
u_ = TrialFunction(V)
F = Identity(len(u)) + grad(u)
I1 = inner(F, F)
J = det(F)

psi = mu/2. * (I1 - 3) - mu * ln(J) + lmbda / 2. * (J-1)**2
pi = derivative(psi, u, v) * dx - inner(bodyForce, v) * dx - inner(surfaceTraction, v) * ds
Jac = derivative(pi, u, u_)

ur = Function(V)
ur.interpolate(rightBoundary)

ul = Function(V)
with ul.vector.localForm() as bc_local:
    bc_local.set(0.0)

bcs = [
    DirichletBC(ul, locate_dofs_geometrical(V, left)),
    DirichletBC(ur, locate_dofs_geometrical(V, right))
]

problem = Hyperelasticity(Jac, pi, bcs, V)
solver = NewtonSolver(MPI.COMM_WORLD)
solver.convergence_criterion = "incremental"
solver.max_it = 25
solver.rtol = 1.e-8
solver.solve(problem, u.vector)

with XDMFFile(MPI.COMM_WORLD, "hyperelasticity.xdmf", "w") as wfil:
    wfil.write_mesh(mesh)
    wfil.write_function(u)

[jhale]

Thanks for this!

It would be even more useful if you submitted it as a pull request to a python/demo/hyperelasticity directory.

[bhaveshshrimali]

Sure, I'd be happy to submit a PR if it's useful. The problem description remains the same, right?

I just saw that a cpp version and ufl already exists. So just a pointer on how to use NewtonSolver (edit: already in Cahn-Hilliard demo) instead of the old solve (feel free to suggest otherwise).