### Laplace + surface Laplace boundary condition

165
views
1
4 months ago by
I would like to solve Poisson's equation on a given domain $\Omega$Ω, but instead of the ordinary Dirichlet or Neumann boundary conditions, I would like to impose  $\Delta_su=0$Δsu=0 on the boundary (where  $\Delta_s$Δs is the surface Laplacian). (The idea here is that if the domain Laplacian and surface Laplacian are 0, the function behaves linearly towards the boundary of the domain, a feature useful for a curve fitting.)

Now, the surface Laplacian problem by itself can be solved by
from dolfin import *
import meshzoo
import numpy

def create_fenics_mesh(points, cells):
# convert to fenics mesh
editor = MeshEditor()
mesh = Mesh()
# topological and geometrical dimension 2
editor.open(mesh, 'triangle', 2, 3, 1)
editor.init_vertices(len(points))
editor.init_cells(len(cells))
for k, point in enumerate(points):
for k, cell in enumerate(cells.astype(numpy.uintp)):
editor.close()
return mesh

points, cells = meshzoo.iso_sphere(ref_steps=4)
mesh = create_fenics_mesh(points, cells)

V = FunctionSpace(mesh, 'CG', 1)

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

b = assemble(Expression('x[0]*x[1]*x[2]', degree=3) * v * dx)

x = Function(V)
solve(A, x.vector(), b)

xdmf = XDMFFile('out.xdmf')
xdmf.write(x)​

See [the article by Rognes et al.](https://www.geosci-model-dev.net/6/2099/2013/gmd-6-2099-2013.pdf).

[This answer](https://www.allanswered.com/post/awlpv/get-test-functions-with-no-support-on-the-boundary/) allows to only get the "interior" equations of the Laplacian of the domain, so the only thing that remains to be done is to combine that with the surface Laplacian.

It's unclear to me however how to do the above if the Functions u and v are really inside the domain. There's of course ds,

dot(grad(u), grad(v)) * ds

but those grads take into account the non-surface component of the gradient as well.

Any hints?

Community: FEniCS Project

4
4 months ago by
You could try subtracting out the normal part of the gradient, i.e., using the projection $P_\Gamma = \mathbf{I} - \mathbf{n}\otimes\mathbf{n}$ from Section 2.2 of

https://arxiv.org/pdf/1509.08597.pdf

to get the tangential gradient $\nabla_\Gamma = P_\Gamma\nabla$.  I'm not sure if that's what you're looking for, but the following at least runs without error for me:

from dolfin import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh,"Lagrange",1)
n = FacetNormal(mesh)
I = Identity(2)

i,j = indices(2)
return as_tensor((I[i,j]-n[i]*n[j])*Dx(u,j),(i,))

# (Some made-up problem to test compilation)
x = SpatialCoordinate(mesh)
f = sin(x[0])*cos(x[1])*x[0]

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