Laplace operator on curved surface

1

Dear FEniCS users,
I’d like to compute the mean curvature on curved surface. I referred to the presentation by Marie E. Rognes here. It seems that differential operators on curved surface such grad(), div() have been implemented in FEniCS.

Take unit sphere for example. I expect the following code will give the mean curvature of unit sphere, which is -1. But the actual results of the following code are very small numbers, like array([ 9.95379322e-12, 9.03699026e-12, -6.24288584e-12, 1.24366276e-11, ...]). Is there something wrong in the code or misunderstanding of how to use the differential operators on curved surface on my part?
Thanks in advance.

from dolfin import *
import pygalmesh
import meshio

s = pygalmesh.Ball([0, 0, 0], 1.0)
mesh = pygalmesh.generate_surface_mesh(
    s, angle_bound=5, radius_bound=0.1, distance_bound=0.1 )
meshio.write("sphere.xml", mesh)
mesh = Mesh("sphere.xml")

order = 3
e_X = Expression(("x[0]", "x[1]", "x[2]"), degree=order)
V0 = VectorFunctionSpace(mesh, "CG", order, dim=3)
V1 = FunctionSpace(mesh, "DG", order-2)
X = project(e_X, V0)
H = project(0.5*dot(div(grad(X)), X), V1)
H.compute_vertex_values()
1 Like
2

See for instance: How to compute curvature of a boundary

3

Hi, Jorgen. Thanks for your quick reply. I don’t know if I get the idea in the linked post. The mesh in my problem is on sphere surface (topological dimension 2). I tried the following code, but with no luck. It gives approximately the same results as those of previous code.

from dolfin import *
import pygalmesh
import meshio

s = pygalmesh.Ball([0, 0, 0], 1.)
mesh = pygalmesh.generate_surface_mesh(
    s, angle_bound=5, radius_bound=0.1, distance_bound=0.1 )
meshio.write("sphere.xml", mesh)
mesh = Mesh("sphere.xml")

order = 3
e_X = Expression(("x[0]", "x[1]", "x[2]"), degree=order)
V0 = VectorFunctionSpace(mesh, "CG", order, dim=3)
V1 = FunctionSpace(mesh, "DG", order-2)
X = project(e_X, V0)

u = TrialFunction(V1)
v = TestFunction(V1)
A = assemble(u*v*dx, keep_diagonal=True)
A.ident_zeros()
L = assemble(0.5*inner(div(grad(X)), X)*v*dx)
H = Function(V1)
solve(A, H.vector(), L)
H.compute_vertex_values()

Previously, I managed to implement Laplace operator vertex-wise on 2D surface according to the paper by M. Meyer et al.. However, I really want to know the reason why the “native” Laplace operator div(grad()) does not work in my code. So I can use it correctly next time.

4

If you want to use my advice, start with a 3D mesh. Use boundarymesh to reduce it to a manifold, then compute the curvature as suggested in that post (You just don’t need to transfer the values back to the volume mesh at the end).

In that way, you can get the FacetNormal for your 2D mesh from the 3D mesh.

You should also normalize your normal, see this example with built in meshes:

from dolfin import *

mesh = SphericalShellMesh.create(MPI.comm_world, 1)

xs = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, "CG", 3)

n = project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs)), V)
    
Z = VectorFunctionSpace(mesh, "DG", 1)
dgn = project(div(grad(n)), Z)

T = FunctionSpace(mesh, "DG", 1)
print(project(0.5*dot(dgn, xs), T).vector().get_local())

gives

[-1.34063378 -1.34037984 -1.34168336 -1.34069641 -1.34027472 -1.34212199
 -1.34069641 -1.34027472 -1.34212199 -1.34063378 -1.34037984 -1.34168336
 -1.34069641 -1.34027472 -1.34212199 -1.34063378 -1.34037984 -1.34168336
 -1.34063378 -1.34037984 -1.34168336 -1.34069641 -1.34027472 -1.34212199
 -1.34040043 -1.34040043 -1.34113709 -1.3404033  -1.3404033  -1.34122741
 -1.3404033  -1.3404033  -1.34122741 -1.34040043 -1.34040043 -1.34113709
 -1.34048022 -1.34048022 -1.34182124 -1.34048022 -1.34048022 -1.34182124
 -1.34047517 -1.34047517 -1.34203122 -1.34047517 -1.34047517 -1.34203122
 -1.34094085 -1.34023178 -1.34214651 -1.34094085 -1.34023178 -1.34214651
 -1.34094085 -1.34023178 -1.34214651 -1.34094085 -1.34023178 -1.34214651]

The spherical mesh is very coarse in this case (I do not have pygalmesh).

5

Thanks for the detailed explanation. I think I can follow your advice to obtain the desired result. Meanwhile, I tried the following code, which gives correct mean curvature.

from dolfin import *

r = 1
sphere = Sphere(Point(0.0, 0.0, 0.0), r)
vmesh = generate_mesh(sphere, 20)
mesh = BoundaryMesh(vmesh, "exterior")

xs = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, "CG", 3)
# x = project(as_vector((xs[0],xs[1],xs[2])), V)  # does not work
x = project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs))*r, V)  # works
n =project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs)), V)
T = FunctionSpace(mesh, "DG", 1)
print(project(0.5*dot(div(grad(x)), n), T).vector().get_local())

What puzzles me is that why

x = project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs))*r, V)

works but

x = project(as_vector((xs[0],xs[1],xs[2])), V)

does not. Why the normalization of coordinates makes all the difference?
I am using FEniCS 2019.1.0 installed via Conda.

6

This should work fine. I don’t use mshr but there is no obvious reason I see for it to not work. Can you post a snippet of the error message ?

I can run your MWE with a mesh generated using pygmsh instead of mshr, i.e. using Gmsh

import pygmsh as pm
geom = pm.opencascade.Geometry()
geom.add_ball([0., 0., 0.], 1)
msh = pm.generate_mesh(geom, dim=3)
import meshio
meshio.xdmf.write("trialMesh.xdmf", meshio.Mesh(
    msh.points, cells = {"tetra": np.vstack((
        cell.data for cell in msh.cells if cell.type == "tetra"))
        }
    ))

from dolfin import *
mesh = Mesh()
XDMFFile("trialMesh.xdmf").read(mesh)
7

Thanks for your response.

from dolfin import *

r = 1.
sphere = Sphere(Point(0.0, 0.0, 0.0), r)
vmesh = generate_mesh(sphere, 20)
mesh = BoundaryMesh(vmesh, "exterior")

xs = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, "CG", 3)
# x = project(as_vector((xs[0],xs[1],xs[2])), V)  # does not work
x = project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs))*r, V)  # works
n =project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs)), V)
T = FunctionSpace(mesh, "DG", 1)
print(project(0.5*dot(div(grad(x)), n), T).vector().get_local())

gives correct mean curvature of sphere with radius 1

[-1.0156043  -1.01669694 -1.01459874 ... -1.01318209 -1.01374275
 -1.01894586]

However,

from dolfin import *

r = 1.
sphere = Sphere(Point(0.0, 0.0, 0.0), r)
vmesh = generate_mesh(sphere, 20)
mesh = BoundaryMesh(vmesh, "exterior")

xs = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, "CG", 3)
x = project(as_vector((xs[0],xs[1],xs[2])), V)  # does not work
# x = project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs))*r, V)  # works
n =project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs)), V)
T = FunctionSpace(mesh, "DG", 1)
print(project(0.5*dot(div(grad(x)), n), T).vector().get_local())

gives

[-3.09199752e-11 -1.15894689e-10 -5.32450640e-11 ... -1.02792785e-10
 -6.85204378e-11  3.19003427e-12]

which are far too small.

8

div(grad(x)) should be identically zero, right?

1 Like
9

x is in the space of degree 3. If div() and grad() are indeed defined on 2D manifold, i.e., spherical surface here, I cannot see why div(grad(x)) has to be identically zero.

1 Like
10

I will take a look over the paper a bit more carefully later. You are right in that it isn’t obvious to have grad(div(x)) to be identically zero over an arbitrary manifold in general (unless I am missing anything). For this case

import pygmsh as pm, numpy as np, meshio
geom = pm.opencascade.Geometry()
geom.add_ball([0,0,0], 1)
msh = pm.generate_mesh(geom, dim=3)
points, cells = msh.points, np.vstack((cell.data for cell in msh.cells if cell.type == "tetra"))
meshio.xdmf.write("s.xdmf", meshio.Mesh(points, cells = {"tetra": cells}))

from dolfin import *
mesh = Mesh()
XDMFFile("s.xdmf").read(mesh)
bmesh = BoundaryMesh(mesh, "exterior")
V = VectorFunctionSpace(bmesh, "CG", 3)
Vd = VectorFunctionSpace(bmesh, "DG", 1)
x = SpatialCoordinate(bmesh)
m1 = project(div(grad(x)), Vd)
m2.vector()[:]   
# array([0., 0., 0., ..., 0., 0., 0.])
11

I think that if you normalize x it stays on the sphere, so you get the right curvature. Instead if you don’t normalize, x lies on the facets of the mesh which are flat so you get locally 0

12

Thanks for the response.
If it is the case, grad() should also give 0. But

from dolfin import *
from mshr import *

r = 1.
sphere = Sphere(Point(0.0, 0.0, 0.0), r)
vmesh = generate_mesh(sphere, 20)
mesh = BoundaryMesh(vmesh, "exterior")
mesh.init(0,1)

xs = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh, "CG", 3)
x =project(as_vector((xs[0],xs[1],xs[2])), V)
x1 = project(as_vector((xs[0],xs[1],xs[2]))/sqrt(dot(xs,xs)), V)
W = TensorFunctionSpace(mesh, "DG", 2)
gradx = project(grad(x), W)
gradx1 = project(grad(x1), W)
print(gradx(Vertex(mesh, 0).point()))
print(gradx1(Vertex(mesh, 0).point()))

gives

[ 0.98866153  0.07453345 -0.0751976   0.07453345  0.51005447  0.49431131
 -0.0751976   0.49431131  0.501284  ]
[ 0.99550091  0.0783516  -0.07244443  0.05990564  0.48328798  0.46998683
 -0.06197578  0.52349369  0.52821502]

Since grad(x) and grad(x1) give approximately the same result, it is likely that div() causes the difference between div(grad(x) and div(grad(x1)). But I don’t know how.

13

I think grad(x) cannot be zero, because x is not constant and lives on the facets of the mesh. But as the normal \mathbf{n} of the facets is locally constant, because the mesh is polyhedral, hence x\cdot \mathbf{n}=const on each facet so that \nabla x\cdot \mathbf{n}=0. Therefore when you compute div (grad(x))\cdot \mathbf{n} you get 0 (the curvature of the polyhedral approximation of the sphere is zero indeed)
whereas \nabla x_1\cdot \mathbf{n}\neq 0 because x_1 is exactly on the sphere, being normalized. It should be easy to check \nabla x\cdot \mathbf{n}=0 and \nabla x_1\cdot \mathbf{n}\neq 0 with Fenics but one must be careful because the FacetNormal(mesh) gives the normal on the boundary, but since the mesh is the boundary of the initial mesh, it is not straightforward to compute these quantities, at least for me, who is not an expert in fenics.
Hope this helps.