How to solve Poisson on the surface of a manifold?

1

Hi. I am struggling with implementing the Poisson equation on the surface of a sphere. In essence I want to solve the equation on a 2D manifold embedded in 3D. The equation is -\bigtriangleup_\Gamma u = f for some f, where \Gamma is the manifold. The problem I have is that the results that FEniCS computes are too large to be correct and hence I get a very large L^2 norm error compared to the correct solution.
The case below in my code is for a sphere with f=2\sin^2\theta (where \theta is the azimuth in spherical coordinates). I used the contravariant Piola transformation to bound the problem to tangential space.

global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree=3)
mesh.init_cell_orientations(global_normal)

order = 3

# Defining the exact solution and approximation spaces
Ve = VectorFunctionSpace(mesh, "Lagrange", order+1)

mesh_cell = ufl.Cell(mesh.cell_name(), geometric_dimension=mesh.geometry().dim())

V = FiniteElement("RT", mesh_cell, order)
Q = FiniteElement("DG", mesh_cell, order-1)
element = MixedElement([V, Q])

W = FunctionSpace(mesh, element)
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

x = SpatialCoordinate(mesh)

# Functions for the problem
theta = acos(x[2])
phi = atan(x[1]/x[0])
f = -2*pow(sin(theta), 2) # Expression("-2*pow(sin(theta),2)", theta=theta, degree=2)

# Solutions
u_sol = sin(theta)*sin(phi)
sigma_sol = Expression((("1-x[0]*x[0]", "-x[0]*x[1]", "-x[0]*x[2]"), ("-x[1]*x[0]", "1-x[1]*x[1]", "-x[1]*x[2]"), ("-x[2]*x[0]", "-x[2]*x[1]", "1-x[2]*x[2]")), degree=3)*grad(u_sol)

# The Forms
a = (inner(tau, sigma) - div(tau)*u + div(sigma)*v)*dx
L = (f*v)*dx

# Solving the problem
uv = Function(W)
solve(a==L, uv)
(sigma, u) = uv.split()

# Computing errors
L2_error_u = errornorm(Expression(str(u_sol), degree=5), u, "L2", degree_rise=1, mesh=mesh)
print(L2_error_u)

I’m sorry for the long code but I don’t know where the error could lie. Any help will be very much appreciated.

2

One issue that stands out immediately is that there is nothing constraining the space of constant functions. If u is a solution to -\Delta u = f on a sphere with no boundaries, then u + C is as well, for arbitrary C\in\mathbb{R}. Thus, there is no unique solution without imposing some auxiliary constraint, such as a DirichletBC acting on a single degree of freedom, or a global constraint on the average value (as demonstrated here).

1 Like
3

@kamensky thank you so much for your reply. I implemented some constraints by increasing the number of trial functions and test functions to get the following form

a = (inner(tau, sigma) + div(tau)*u + div(sigma)*v + r*v + t*u)*dx

But I have read a paper that seems to imply that this should be the way I should do the calculation on any mesh with geometric dimension 3 and topological dimension 2:

a = inner(grad(u), grad(v))*dx

This on the other hand gives me very strange results. Do you know what the correct way to implement it is?

4

The second formulation, a = inner(grad(u),grad(v))*dx, also requires an additional constraint to have a unique solution. It’s also possible you may have forgotten to switch the space for u to a continuous Galerkin space ("CG") when testing the second formulation; it would be ill-posed with the DG space used for u and v in your mixed formulation, and would produce strange results, if any.