Eigenvalues Laplacian on the 2D half sphere

1

Hi Fenics users,

It’d be great if you can help a newbie here with the following problem. I’m trying to obtain numerical eigenvalues of the Laplacian on a 2D semisphere (\mathbb{S}^2_+) with zero Dirichlet boundary conditions applied on the equatorial circle.

So far I’ve only achieved to obtain eigenvalues for the Neumann problem on the full sphere. I’m failing at describing correctly both the half sphere and the boundary conditions. Even the simplest case with zero boundary conditions on the whole equator isn’t working.

I’ve tried to use the function SubDomains to express the problem but fenics gives me the warning: “Found no facets matching domain for boundary condition.” I attach below the code for the Neumann problem in case you could point out what should I add to implement the boundary conditions. In particular, I’d like to impose that the eigenfunctions are zero on an arbitrary arc segment on the equator. Say, for instance, like the one depicted in the figure below where we have Dirichlet conditions on half circle and Neumann conditions on the other half (here the picture shows only a projection):

arc
arc540×539 15.1 KB

The code:

from dolfin import *
from mshr import *

            
# Mesh and function space
sphere = Sphere(Point(0.0, 0.0, 0.0), 1.0)
vmesh = generate_mesh(sphere, 10)
mesh = BoundaryMesh(vmesh, "exterior")
V = FunctionSpace(mesh, "CG", 1)

#This part of the code should implement the zero Dirichlet problem but doesn't work.

#def boundary(x, on_boundary):
#   return (on_boundary and near(x[0],0))

#class SemiSphere(SubDomain):
#    def inside(self, x, on_boundary):
#      return (on_boundary and x[0]>0.0)

#bc = DirichletBC(V, 0.0, boundary)

# Basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
m = u*v*dx

# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)

# Boundary conditions
#bc.apply(A)         
#bc.apply(M)

# Eigensolver
eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'

# Eigenvalues computation
print("Computing eigenvalues")
eigensolver.solve()
neigs=4

assert eigensolver.get_number_converged() > 0
for i in range(0,neigs):
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print("eigenvalue: ",i+1, r)

Thanks in advance!