Specifying 'F' for Mixed finite element for Nearly Incompressible Elasticity

Hi

I'm trying to modify source code that Corrado Maurini provided a while back. In it, the strain energy (W) was provided and 'F' was calculated automatically using F = derivative(W,.....).

However, I would like to specify 'F' manually as I'd like to manipulate the 2PK stress tensor. It's a two-field mixed formulation with u (displacement) and p (pressure) as unknowns.

From literature, 'F' should be:

DPi(u,p).du = integral( [J(u)*p*C(u)^-1 + 2*(dW/dC) ] : dE(u) )dV + DPi_ext
DPi(u,p).dp = integral( [J(u) - 1]*dp )dV

My attempts result in either the nonlinear variational solver producing nans or PETSC ERRORs. Below, I've included the code. Could I please get some help as to what I'm doing wrong.

Many thanks

Pietro
=======

from dolfin import *

# Create mesh and define function space
mesh = UnitSquare(20,20)
#plot(mesh)

# Mark boundary subdomains
def left(x, on_boundary):
    return x[0]==0 and on_boundary

b=Constant((0.,.04))

# Material parameters
Y, nu = 1.0, 0.49999
mu, lmbda = Constant(Y/(2*(1 + nu))), Constant(Y*nu/((1 + nu)*(1 - 2*nu)))

# Defintion of function spaces
porder=2
Pu = VectorFunctionSpace(mesh, "CG", porder) # space for displacements
Pp = FunctionSpace(mesh, "CG", porder-1) # space for pressure
V = MixedFunctionSpace([Pp,Pu]) # mixed space
ndim = Pu.cell().d # number of space dimensions

# Create test and trial functions
#
vq = TestFunction(V) # Testfunction in the mixed space
dup = TrialFunction(V) # Trial function in the mixed space to define the incremental solution
du, dp = split(dup)

up = Function(V) # Function in the mixed space to store the solution
dupsol = Function(V) # Function in the mixed space to store the solution of the linearized problem
p, u = split(up) # Function in each subspace to write the functional
q, v = split(vq)

# Boundary conditions
bcL= DirichletBC(V.sub(1), Constant((0.0,0.0)), left)

# Kinematics
I = Identity(ndim) # Identity tensor
F = I + grad(u) # Deformation gradient
C = transpose(F)*F # Right Cauchy-Green tensor
E = 0.5*( C - I ) # Green-Lagrange tensor

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Stored strain energy density (compressible neo-Hookean model, mixed formulation)
W = (mu/2)*(Ic - ndim) - mu*ln(J) + p*ln(J) - 1/(2*lmbda)*p**2
# Potential energy
#Pi_bulk = W*dx - inner(b,u)*dx

# First derivatives of the bulk and surface contributions to the potential energy
#F_bulk=derivative(Pi_bulk,up,vq)

#========================================================
#Specifying F Manually
#========================================================
dE = 0.5 * ( grad(v).T*F + F.T*grad(v) )

CC = variable(C)
S = J*p*inv(C) + 2*diff(W,CC)

F_bulk = inner( S, dE )*dx + inner(b,v)*dx + (J-1)*q*dx
#========================================================

## Second derivatives of the bulk and surface contributions to the potential energy
dF_bulk=derivative(F_bulk,up,dup)

# Set-up the FEniCs build-in Newton
MyNonlinearVariationalProblem = NonlinearVariationalProblem(F_bulk, up, bcL, dF_bulk)
solver = NonlinearVariationalSolver(MyNonlinearVariationalProblem)
solver.parameters["linear_solver"]="petsc" # use PETSc builtin LU solver
# Solve the problem
solver.solve()
# Plot and hold solution
plotu=plot(u, mode = "displacement",axes=True,interactive=True)
#################################################################################