3 Field Hu-Washizu Elasticity

Question

I have tried to implement a 3 field variational form of nearly incompressible elasticity with the following code

from dolfin import *
mesh = UnitCubeMesh(3,3,3)

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, Q, Q])

clampBC = DirichletBC(W.sub(0), (0,0,0), "near(x[0], 0)")
pushBC = DirichletBC(W.sub(0), (0, 0.1, 0), "near(x[0], 1.0)")

w = Function(W)
u,p,d = split(w)

F = grad(u) + Identity(3)
J = det(F)
C = J**(-2.0/3.0)*F.T*F

C1 = Constant(2.5)
lmbda = Constant(30.0)

psi = C1*(tr(C) - 3) + p*(J - d) + lmbda*(ln(d)**2)

R = derivative(psi*dx, w, TestFunction(W))
J_form = derivative(R, w, TrialFunction(W))

solve(R == 0, w, [clampBC, pushBC], form_compiler_parameters = {"keep_diagonal": True})
plot(split(w)[0], interactive = True, mode = "displacement")

When I run the code I get

Solving nonlinear variational problem. Newton iteration 0: r (abs) =
-nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)

Any help is appreciated.

Answer 1

Magnus Nordaas pointed out that the ln(d) term blows down if d starts at 0. This version of the script works now.

mesh = UnitCubeMesh(3,3,3)

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, Q, Q])

clampBC = DirichletBC(W.sub(0), (0,0,0), "near(x[0], 0)")
pushBC = DirichletBC(W.sub(0), (0, 0.1, 0), "near(x[0], 1.0)")

w = Function(W)
u,p,d = w.split()

fa = FunctionAssigner(W.sub(2), Q)
fa.assign(d, interpolate(Constant(1.0), Q))

u,p,d = split(w)

F = grad(u) + Identity(3)
J = det(F)
C = J**(-2.0/3.0)*F.T*F

C1 = Constant(2.5)
lmbda = Constant(30.0)

psi = C1*(tr(C) - 3) + p*(J - d) + lmbda*(ln(d)**2)

R = derivative(psi*dx, w, TestFunction(W))
J_form = derivative(R, w, TrialFunction(W))

solve(R == 0, w, [clampBC, pushBC], form_compiler_parameters = {"keep_diagonal": True})
plot(split(w)[0], interactive = True, mode = "displacement")