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