# Solve Poisson equation without integrating by parts

Hello,
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)
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

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.