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 \Omegawithout 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)
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.
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.
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.
