Approximate $\|\Delta f\|^2_{L^2(\Omega)}$ in finite element context

Question

I have minimization problem of the form $$ G(f) + \|\Delta f\|^2_{L^2(\Omega)} \to \min $$ over all $f\in C^2(\Omega)$, $\Omega$ being closed and bounded.

Let us forgot about $G$; I'm interested in how to discretize the $L^2$-norm. In the context of finite volumes, one has $$ \begin{multline} \int_\Omega (\Delta f)^2 = \sum_i \int_{V_i} \left(\Delta f\right)^2 \approx \sum_i |V_i| (\Delta f)(x_i)^2\\ \approx \sum_i |V_i| \left(|V_i|^{-1} \int_{V_i} \Delta f\right)^2 = \sum_i |V_i|^{-1} \left(\int_{V_i} \Delta f\right)^2, \end{multline} $$ so $$ \|\Delta f\|_{L^2(\Omega)} \approx \|\Delta_h P_h(f)\|_{M^{-1}} $$ where $M$ is the mass matrix and $P_h$ the projection into the discretized space.

It seems from numerical experiments (see below) that this is also true for finite-element-type functions (withe the respective mass matrix; see the code below).

Does anyone know why that is?

---

import sympy
from dolfin import (
    Expression,
    FacetNormal,
    Function,
    FunctionSpace,
    TestFunction,
    TrialFunction,
    UnitSquareMesh,
    assemble,
    dot,
    ds,
    dx,
    grad,
    project,
    solve,
)

mesh = UnitSquareMesh(500, 500)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)

n = FacetNormal(mesh)
A = assemble(dot(grad(u), grad(v)) * dx - dot(n, grad(u)) * v * ds)
M = assemble(u * v * dx)

f = Expression("sin(pi * x[0]) * sin(pi * x[1])", element=V.ufl_element())
x = project(f, V)

Ax = A * x.vector()
Minv_Ax = Function(V).vector()
solve(M, Minv_Ax, Ax)
val = Ax.inner(Minv_Ax)

print(val)


# Exact value
x = sympy.Symbol("x")
y = sympy.Symbol("y")
f = sympy.sin(sympy.pi * x) * sympy.sin(sympy.pi * y)
f2 = -sympy.diff(f, x, x) - sympy.diff(f, y, y)
val2 = sympy.integrate(sympy.integrate(f2 ** 2, (x, 0, 1)), (y, 0, 1))
print(sympy.N(val2))

Output:

97.75031146783857
97.4090910340024

Answer 1

The problem, of course, is that the typical finite element shape functions are not twice differentiable and so $\Delta f$ doesn't exist as a function of which you can easily take the norm. To be more precise, the typical shape functions are continuous, and so the derivative is discontinuous across cell interfaces, and the second derivative therefore has an easily integrable part in the cell interiors, and a distribution-valued part that lives on the interfaces. This part can not be squared and integrated.

So the objective function you want to minimize is not well defined on typical finite element spaces. You will either have to restate the problem in terms of a norm that is well defined and somehow will have to involve the jump in the gradient of the solution, or to approximate the gradient of the solution by a continuous function whose gradient you can square and integrate.