### Laplace + surface Laplace boundary condition

165

views

1

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$Δ

Now, the surface Laplacian problem by itself can be solved by

_{s}`u`=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):
editor.add_vertex(k, point)
for k, cell in enumerate(cells.astype(numpy.uintp)):
editor.add_cell(k, cell)
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)
A = assemble(dot(grad(u), grad(v)) * dx)
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 `grad`s take into account the non-surface component of the gradient as well.

Any hints?

Community: FEniCS Project

### 1 Answer

4

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:

Ah yes, this is the same as grad(u) - dot(grad(u), n) * n. Thanks for the hint!

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)
def grad_Gamma(u):
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)
res = inner(grad(u),grad(v))*dx + inner(u,v)*dx \
+ inner(grad_Gamma(u),grad_Gamma(v))*ds - inner(f,v)*dx
u = Function(V)
solve(lhs(res)==rhs(res),u)
```

written
4 months ago by
Nico Schlömer

Please login to add an answer/comment or follow this question.