Mixed formulation problem

Question

Hi, I was trying to solve a simple incompressible neo-hookean hyperelasticity problem using mixed formulation but I somehow ran into the following error.

problem = NonlinearVariationalProblem(F1, w, bcs, Jac)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 145, in init
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 53, in init
"Expected a ufl.Form.")
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/common.py", line 2611, in dolfin_error
return _common.dolfin_error(location, task, reason)
RuntimeError:

---

DOLFIN encountered an error. If you are not able to resolve this issue
using the information listed below, you can ask for help at

fenics@fenicsproject.org

Remember to include the error message listed below and, if possible,
include a minimal running example to reproduce the error.

---

Error: Unable to creating dolfin.Form.
Reason: Expected a ufl.Form..
Where: This error was encountered inside form.py.
Process: 0

DOLFIN version: 1.6.0
Git changeset: unknown

---

Been spending my time looking at the code but couldn't find where exactly I did wrong. Any help would be greatly appreciated. Thanks a bunch for your help in advance!! My code are given below

from dolfin import *

def Fmat(u):
    d = u.geometric_dimension()
    I = Identity(d)             # Identity tensor
    F = I + grad(u)             # Deformation gradient
    return F

# Define boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[2] - 0.0) < DOLFIN_EPS and on_boundary

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[2] - 1.0) < DOLFIN_EPS and on_boundary

class Lower(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1] - 0.0) < DOLFIN_EPS and on_boundary

class Upper(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1] - 1.0) < DOLFIN_EPS and on_boundary

parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


nx = 1  # divisions in r direction
ny = 1  # divisions in theta direction
nz = 1  # divisions in theta direction
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), nx, ny, nz)

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG",1)
W = V*Q
# Initialize sub-domain instances
left = Left()
right = Right()
bottom = Bottom()
top = Top()
lower = Lower()
upper = Upper()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
bottom.mark(boundaries, 3)
top.mark(boundaries, 4)
lower.mark(boundaries, 5)
upper.mark(boundaries, 6)

ds = Measure("ds")[boundaries]
n = FacetNormal(mesh)

# Set up boundary condition at left end
c = Constant((0.0))
bcleft0 = DirichletBC(W.sub(0).sub(0), c, boundaries, 1)
bcleft1 = DirichletBC(W.sub(0).sub(1), c, boundaries, 1)
bcleft2 = DirichletBC(W.sub(0).sub(2), c, boundaries, 1)
bclower0 = DirichletBC(W.sub(0).sub(0), c, boundaries, 5)
bclower1 = DirichletBC(W.sub(0).sub(1), c, boundaries, 5)
bclower2 = DirichletBC(W.sub(0).sub(2), c, boundaries, 5)
bcbottom2 = DirichletBC(W.sub(0).sub(2), c, boundaries, 3)

uspecified = Expression("u", u = 0.0)
bc_right0 = DirichletBC(W.sub(0).sub(0), uspecified, boundaries, 2)
bc_right1 = DirichletBC(W.sub(0).sub(1), c, boundaries, 2)
bc_right2 = DirichletBC(W.sub(0).sub(2), c, boundaries, 2)

# Set up boundary conditions
bcs = [bcleft0, bcleft1, bcleft2, bc_right0, bc_right1, bc_right2]

## Define variational problem
w = Function(W)
(u,p) = w.split();

C = Fmat(u).T*Fmat(u)                   # Right Cauchy-Green tensor
Ic = tr(C)
J = det(Fmat(u))
d = u.geometric_dimension()
S = Identity(d)
Pmat = Fmat(u)*S
psi = (Ic - 3) - p*(J - 1)
F1 = derivative(psi, w, TestFunction(W))

# Compute Jacobian of F
Jac = derivative(F1, w, TrialFunction(W))

problem = NonlinearVariationalProblem(F1, w, bcs, Jac)

Answer 1

The first argument of the NonlinearVariational problem must be a ufl.Form object but in your code F1 is a ufl.differentiation.CoefficientDerivative object.

The solution to your probem is to use

F1 = derivative(psi*dx, w, TestFunction(W))

instead of

F1 = derivative(psi, w, TestFunction(W))

The above suggestion works because you want to minimize the total potential energy psi*dx and in your problem psi is the elastic stored energy density (see this demo as a reference).

Comments

Thanks a lot!! I should have notice that ..........
By the way, I was wondering why is it that it is much slower and require more quadrature points if I started off with the functional F1 i.e.,

F1 = inner(grad(v),Pmat - J*p*invF.T)*dx + inner(q,J - 1)*dx

I can get the same answer but starting off with F1 is much much slower rather than having

F1 = derivative(psi*dx, w, TestFunction(W))