Trouble finding the hessian matrix on a surface

1

Hello,

I’m trying to compute a Hessian matrix of a function on the surface of a sphere. I’ve followed Rognes et al on how to set up differential operators on a 2D surface embedded in 3D space. The gradient operator converges how I would expect, so that’s good. But when I try to compute the Hessian matrix as described here and here, it converges to the wrong thing (in the case of the function included below, it converges to the cartesian hessian matrix, but I’ve found it’s not always the cartesian hessian for other functions).

If anyone has any suggestions, I would be extremely grateful! Thanks!

MWE:

from __future__ import print_function

# Environment setup
import os

# FEniCS imports
from fenics import *
from dolfin import *
import numpy

# Read the mesh and get boundaries

mesh = Mesh()
with XDMFFile("meshes/128/mesh.xdmf") as infile:
    infile.read(mesh)
# Define a mesh value collection of dim = 3, for the cell data
# read boundaries
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("meshes/128/mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

surface = BoundaryMesh(mesh, "exterior")

# Define Function Spaces
Vt = TensorFunctionSpace(surface, 'DG',0)
Vs = FunctionSpace(surface,'CG',2)

w = TrialFunction(Vt)
v = TestFunction(Vt)

# Function to differentiate
el = Vs.ufl_element()
uex = Expression("sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])",element=el)
u = interpolate(uex,Vs)

# Define problem
a = inner(w,v)*dx
L = inner(grad(grad(u)),v)*dx

# Solve Problem
hess_u = Function(Vt)
solve(a == L, hess_u)

hess_u.rename("hessian", "hess_u")
file = File('hess_u.pvd')
file << hess_u
2

My guess is that this has something to do with the fact that the geometry is still faceted, even if you use a higher-order function space to define u. Maybe you could try solving a variational problem for a weak definition of the Hessian, e.g.: Find a tensor-valued function \mathbf{H} such that

\int_\Omega\mathbf{H}:\mathbf{G}\,d\Omega = -\int_\Omega(\operatorname{grad} u)\cdot(\operatorname{div}\mathbf{G})\,d\Omega~~~~~\forall\mathbf{G}\text{ ,}

assuming a closed manifold with no boundary (so there is no boundary term).

3

With the same finite element spaces (CG2 and DG0 for the scalar and tensor spaces, respectively) using this form returns zero for whatever function I input even if the surface hessian is nonzero. For reference, two other functions I’ve been testing with are u(\vec x)=x and u(\vec x)=x^2.

If I change both spaces to P1, I do get nonzero results for these functions, but a) it’s very noisy and b) even if I do add a term to smooth out the noise it’s converging to something that isn’t the analytic solution.

I also tried splitting the equation as such:
\begin{equation} \int_\Omega \vec h \cdot \vec g = \int_\Omega \nabla u \cdot \vec g \end{equation}
\begin{equation} \int_\Omega H : G = \int_\Omega \nabla \vec h : G \end{equation}
and both solving these simultaneously and then also tried solving them separately to no avail.