Dear FEniCS Community,

I am a new user and this is my first time posting a question here: I apologize in advance for any mistakes or if my questions are too basic.

I am trying to use FEniCS to solve a time-dependent hyperelasticity problem, using a mixed formulation in which the unknowns are the displacement (a vector field **u**) and the pressure (a scalar field *p*). Starting from the hyperelasticity demo, I attempted to write my own code, but I have an issue: the program seems to correctly solve the problem at the first time step, but then nothing else happens, and the solution is not updated in the following iterations of the time loop.

Here is a simplified version of the code:

```
import dolfin
from dolfin import *
dolfin.parameters["form_compiler"]["cpp_optimize"] = True
dolfin.parameters["form_compiler"]["representation"] = "uflacs"
from ufl import grad as ufl_grad
def Grad(v):
return ufl_grad(v)
mesh = UnitCubeMesh(6, 6, 6)
#Elements and Function Spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2) #displacement u
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1) #pressure p
TH = MixedElement([P2, P1]) #u + p
V = FunctionSpace(mesh, TH) #(u,p)
# Mark boundary subdomains
upp = CompiledSubDomain("near(x[2], side) && on_boundary", side = 1.0)
down = CompiledSubDomain("near(x[2], side) && on_boundary", side = 0.0)
#Define parameters
u0 = Constant((0.0, 0.0, 0.0))
pp = Constant(0.0)
T = 1e-08
num_steps = 100
dt = T / num_steps
TT = Constant((0.0, 0.0, -0.1)) # Traction force on the boundary - current configuration
k = Constant(5e-9)
K0 = as_matrix([ [k, 0, 0], [0, k, 0], [0, 0, k] ])
#Boundary Conditions: zero displacement at z = 0, zero pressure at z = 1
bcd_z = DirichletBC(V.sub(0), u0, down)
bcp = DirichletBC(V.sub(1), pp, upp)
bcs = [bcd_z, bcp]
# Define functions
dup = TrialFunction(V) # Incremental displacement and pressure
du, dp = split(dup)
u_, p_ = TestFunctions(V)
up = Function(V) # Unknown Displacement and Pressure (n+1)
u, p = split(up)
up_prev = Function(V) # Previous iteration
u_prev, p_prev = split(up_prev)
up_init = Expression( ("0.0", "0.0", "0.0", "0.0"), degree=1)
up_prev = interpolate(up_init, V) #Initial condition for u and p
# Kinematics
d = len(u)
I = Identity(d)
F = I + Grad(u) #Deformation gradient
#Jacobian of F
def Js(u):
return det(F)
#First Piola-Kirchhoff stress tensor - Mooney Rivlin
def P(u, p):
C = (F.T)*F
Ic = variable(tr(C))
IIc = variable(0.5*(tr(C)**2 - tr(C*C)))
#Mooney-Rivlin strain energy
mu1 = Constant(0.4)
mu2 = Constant(0.25)
psi = (mu1/2)*(Ic - 3) + (mu2/2)*(IIc - 3)
gamma1 = diff(psi, Ic) + Ic*diff(psi, IIc)
gamma2 = -diff(psi, IIc)
P_s = 2*F*(gamma1*I + gamma2*C)
return P_s - Js(u)*p*inv(F).T
#Permeability tensor
def K_star(K0, u):
return Js(u)*inv(F)*K0*inv(F).T
#Traction in material configuration
def T_star(TT, u):
return Js(u)*inv(F).T*TT
dolfin.parameters["form_compiler"]["quadrature_degree"] = 6
#Loop for time stepping
t = 0
n = 1
while (t <= T):
t += dt
print("Iteration", n, "-th", "Time", t)
#Variational form
L = Js(u)*p_*dx + dt*inner(Grad(p_), K_star(K0, u)*Grad(p))*dx - inner(P(u, p), Grad(u_))*dx - Js(u_prev)*p_*dx + inner(T_star(TT, u), u_)*ds
j = derivative(L, up, dup)
problem = NonlinearVariationalProblem(L, up, bcs, J=j)
solver = NonlinearVariationalSolver(problem)
solver.solve()
assigner = FunctionAssigner(V, V)
assigner.assign(up_prev, up)
n = n+1
```

It looks as if the nonlinear variational form L (or its jacobian) is not updated properly after the first time step: the first iteration is fine and the solution behaves as expected, but then it remains constant. I tried to make use of “FunctionAssigner” instead of the simple “assign” command as suggested in a similar topic I found here, but unfortunately it did not solve my problem. Defining everything inside the time loop does not seem to fix the issue as well.

Any help or suggestion would be greatly appreciated. Thank you all in advance!