Solve Poisson equation without integrating by parts

I am dealing with a problem where I want to avoid by-part integration in the variational formulation. To make a simple test, I want to solve the Poisson equation with Dirichlet boundary conditions
\nabla^2 u = -f \textrm{ in } \Omega
u = u_D \textrm{ on } \partial \Omega

in the unit square \Omega without integrating by parts. In other words, in the variational formulation I leave the term

\int_\Omega (\nabla^2 u) v

as it is. Here is a minimal example:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np

xdmffile_u = XDMFFile("u.xdmf")

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 8)

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)
a = div(grad(u))* v*dx
L = -f*v*dx

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

xdmffile_u.write(u, 0)

However, the solution is strange and incorrect:

What am I doing wrong? Thank you :smiley:

There’ll be other other issues to consider, I expect, but at a basic level you’re proposing to take the second derivative of linear functions (P1 elements),which of course is 0. So not surprising that it doesn’t work. Try it with higher order elements.

Hello, the same issue arises if I increase the order of the polynomials, see the revised original post.

Dumb question / suggestion:

Are you not supposed to use different variable names for trial function, u, and the solution function, u? I realize that you define bilinear and linear forms prior and I don’t have FEniCSx access at work so I couldn’t try it.

1 Like

Please note that from an analysis point of view, what you are doing here is not well defined, as:
\nabla^2u is not well defined over \Omega for any H^1 finite element, which is what we support in FEniCS. This is because any function from a higher order Lagrange space is only part of H^1, while the equation you are considering is only valid in H^2. The whole benefit of the FEM is to reduce the regularity of the solution from H^2 to H^1.