Generate Neo-Hookean force & force Jacobian functions using FEniCS

neohookean.py

from ufl import *
import ffc.compiler

# Define cell type we solve on
cell = tetrahedron
d = cell.geometric_dimension()

# The finite element function space
element = VectorElement("Lagrange", cell, 1, dim=d)

# Displacement that is used as linearization point
u = Coefficient(element)
# Increment that is solved for in each iteration
du = TrialFunction(element)
# Test functions
v = TestFunction(element)

# Material parameters
rho = Constant(cell)
mu = Constant(cell)
lmbda = Constant(cell)

I = Identity(d)
# Deformation gradient
F = I + grad(u)
F = variable(F)
# Right Cauchy-Green strain tensor
C = F.T * F

# Hyperelastic invariants
Ic = tr(C)
J = det(F)

# Neo-Hookean strain energy density
psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
# First Piola-Kirchhoff stress tensor
P = diff(psi, F)

# Elastic forces
L = inner(P, grad(v)) * dx
# Jacobian of the elastic forces
a = derivative(L, u, du)

# Mass matrix
M = rho*inner(du, v) * dx

# Total kinetic energy
E_kin = (rho/2) * inner(u, u) * dx
# Total strain energy
E_strain = psi*dx

parameters = ffc.default_parameters()
parameters["representation"] = "uflacs"
parameters["optimize"] = True

code_c, code_h = ffc.compiler.compile_form([L], prefix="neohookean", parameters=parameters)

with open("neohookean.c", mode="w") as f:
    f.write(code_c)