difficulty with grad(div u) problem

Question

Hello all,

I am trying to solve the following problem:

PDE: Laplacian(u) + grad(div u) = f subject to Homogeneous Dirichlet boundary condition for u over the domain given by a unit square.

Code:

from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(101,101)
V= FunctionSpace(mesh, 'P', 1)
def boundaryX0(x):
    return x[0] <DOLFIN_EPS and x[0]> 1.0-DOLFIN_EPS and x[1]<DOLFIN_EPS and x[1]> 1-DOLFIN_EPS
u0= Constant(0.0)
bc= DirichletBC(V, u0, boundaryX0) 
f= Constant(1.0)
v=TestFunction(V)
u= TrialFunction(V)
a= inner(curl(u), curl(v))*dx - dot(grad(u),grad(v))*dx - dot(grad(u),grad(v))*dx
L=-inner(f,div(v))*dx
u=Function(V)
solve(a==L,u,bc)

I am getting the following error:

Calling FFC just-in-time (JIT) compiler, this may take some time.
Component and shape length don't match.
Traceback (most recent call last):

RuntimeError:
*** -------------------------------------------------------------------------
*** Error: Unable to perform just-in-time compilation of form.
*** Reason: ffc.jit failed with message:
UFLException: Component and shape length don't match.
.
*** Where: This error was encountered inside jit.py.
*** Process: 0

*** DOLFIN version: 2016.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Please help me out in finding the problem.

Thank you.

Ratnesh

Comments

Hi, please try to paste your code into a formatted block like the following. Some of it has been lost in the conversion.

 from dolfin import *
 import numpy as np
 # more...

Also, you seem to define a problem with a vector valued solution. Search for VectorFunctionSpace in the docs and this forum. There seems to be some mistake in the definition of the bilinear form too (duplicated code).

---

Thanks, I will keep that in mind for next time.

Can you please point out what the problem with the bilinear form is

---

I meant the following snippet in the bilinear form

  - dot(grad(u),grad(v))*dx - dot(grad(u),grad(v))*dx

which I guess is not what you intended since you could just multiply by a Constant(2.0).

---

Yes, I see that now.

Thanks.

Answer 1

In your code, you write

inner(curl(u), curl(v))dx - dot(grad(u),grad(v))dx - dot(grad(u),grad(v))dx

However, u is defined as a scalar function instead of a vector function, so the function curl makes no sense.

Comments

Yes, it should be VectorFunctionSpace. Thanks