# Weak formulation for square of gradient

Hello

I’m a beginner in FEniCS and have been following the “Solving PDEs in Python” tutorial (which is excellent - thank you!).

Considering e.g. the simple Poisson example problem: -nabla^2 u(x) = f(x). How do you deal with terms like (nabla u(x))^2?

E.g. -nabla^2 u(x) + (nabla u(x))^2 = f(x), by which I mean in 1d -u''(x) + (u'(x))^2 = f(x). Is this possible? How do you write the weak form - do I need to integrate by parts?

My attempt: starting from the provided example which works:

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)

# Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)


If i don’t integrate by parts, I should just modify it as follows:

a = dot(grad(u), grad(v)) * dx + dot(grad(u), grad(u)) * v * dx


This gives the error:

NotImplementedError: Cannot take length of non-vector expression.


Where am I going wrong? Or is it a conceptual mistake?

Thanks

1 Like

Welcome robert. First: You can use LaTeX for formulas using \$ signs.

The weak form is as usual obtained through multiplying by test function v and integration (assuming homogenous Neumann boundary conditions)

You are on the right track in your approach, but the syntax solve(a==L, ...) is used for linear problems and the one you want to solve is nonlinear. But we can simply define a nonlinear variational form as

such that your solution has to satisfy F(u,v)=0. FEniCS will automatically compute the Jacobian and iteratively solve your problem when solve is called in the right way (cf. my answer here: Function() v.s. TrialFunction())

Consider the following code:

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)

# Define variational problem

u = Function(V)
v = TestFunction(V)
f = Constant(-6.0)