unable to define nonlinear variational problem

Question

I'm solving a mixed problem and I define the energy function and using derivative() to obtain the weak form according the hyperelasticity example(https://fenics.readthedocs.io/projects/dolfin/en/latest/demos/demo_hyperelasticity.py.html#demo-hyperelasticity). But I'm not able to get it right somehow. Below is the minimal working code snippet. Any help is greatly appreciated!

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(32,32)

U = VectorElement("Lagrange",mesh.ufl_cell(), 1, dim=2) # displacement
M = FiniteElement("Lagrange",mesh.ufl_cell(), 1)       # mu
W = VectorElement("Lagrange",mesh.ufl_cell(), 1, dim=2) # Lagrangian multipliers
S = FunctionSpace(mesh, MixedElement([U,M,W]))

v, q, w = TestFunctions(S)
u, mu, lamd = TrialFunctions(S)
L = inner(u,u)*dx + inner(grad(mu),grad(mu))*dx + \
inner(mu*(grad(lamd)+grad(lamd).T), grad(u)+grad(u).T)*dx + inner(div(lamd),div(u))*dx

ds = TestFunction(S)
du = TrialFunction(S)
s = Function(S)
F = derivative(L, s, ds)
JF = derivative(F, s, du)

bc = DirichletBC(S.sub(0), Constant([0.0,0.0]), "on_boundary")

problem = NonlinearVariationalProblem(F,s,bc,JF)
solver = NonlinearVariationalSolver(problem)
solver.solve()
u_s, mu_s, lamd_s = s.split()

The error message is
*** Error: Unable to define nonlinear variational problem F(u;v)=0 for all v.
*** Reason: Expecting the residual F to be a linear form (not rank 2).

By the way, I tried this on both fenics 2016.1.0 and 2016.2.0 and got the same error.

Answer 1

Hello,

I've done something similar to this type of problem before, and ran into some similar issues. I'm not on my work computer though so bear with me. From what I can remember I think there are a few issues with your code. Firstly, your variable you solve for, s, isn't tied in any way to your weak formulation. Secondly, the nonlinear solver doesn't use a Trial Function, only a Test Function (I think this is the cause of the error you got). Saying something to the effect of,

s = Function(S)
u, mu, lamd = s.split()

right before you define your weak form (instead of defining the Test Functions) may help. I can double check with some old code of mine tomorrow.

edit: I just tried it with my suggestion, and the nonlinear solver seems to work now. It doesn't converge now to an actual value but I think this is a different problem (or it may have something to do with the fact you only posted just a snippet of your actual code, not sure).

Comments

Hi,

Thanks a lot for your reply. Yeah, you are right. Your suggestion works. It diverges now maybe because my zero initial guess. I'm working on that now.

---

After fixing the problem according to your suggestion, now it's able to finish the first iteration but then it throws out a petsc error I'm not able to comprehend. I think that might be another question. So I asked another one here. I'm not sure if you also have experience on this one. Thanks in advance!