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):
        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
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)

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)
​
Ah yes, this is the same as grad(u) - dot(grad(u), n) * n. Thanks for the hint!
written 4 months ago by Nico Schlömer  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »