Calculating error on surface mesh

1

Hi again,

I am trying to verify that my method behaves as intended. I want to solve equations on the sphere, and below is my attempt to verify the method through solving the Laplace-Beltrami equation, with the 2,2-real spherical harmonic times its eigenvalue as right-hand side, which is an eigenfunction, so I know the result. Theoretically I know that the error of the method should behave as h^2, where h is the max size of the triangle. However, what comes out behaves strangley. I think what I do can be relatated to how I calculate the error. See below, errornorm gives warning and should according to old discussions, ~2016 not be used. Instead I try to compile the error, but it does not work at all. I ideally want to integrate over “sphere”, but that does not seem to work. I think I have to compare them on the same mesh, so I choose the finest, like I do below, but it does not work. Any suggestions? And when should I use ds and not dx in fenics? Best regards!

from mshr import *
from dolfin import * 
import numpy as np 
import matplotlib.pyplot as plt
#Domain
N=[8,16,32,64]
n=len(N)
print(n)
errors=np.zeros([n,1])
sizes=np.zeros([n,1])

sphere=Sphere(Point(0,0,0),1)
global_normal=Expression(('x[0]','x[1]','x[2]'),degree=2)
surfaces=[BoundaryMesh(generate_mesh(sphere,i),"exterior") for i in N]

ref_surface=surfaces[-1]ref_surface=surfaces[-1]
f = Expression("6*sqrt(15)*(x[0]*x[0]-x[1]*x[1])/(16*sqrt(pi))",degree=2) #error 
Y22=Expression("sqrt(15)*(x[0]*x[0]-x[1]*x[1])/(16*sqrt(pi))",degree=2)
file=File("./testim/ver.pvd")


def verificator(N,surface):
    surface.init_cell_orientations(global_normal) 
    P1 = FiniteElement("P", surface.ufl_cell(), 1)
    C=FiniteElement("R",surface.ufl_cell(),0)
    V = FunctionSpace(surface,P1*C)
    # Define variational problem
    (u,c) = TrialFunction(V)
    (v,d) = TestFunction(V)
    a = inner(grad(u), grad(v))*dx+(c*v+u*d)*dx
    L = f*v*dx 
    # Compute solution
    w = Function(V)
    solve(a == L, w)
    (x,y)=w.split()
    eig=6 #l*(l+1)
    file<<x
    size=surface.hmax() #have to do for now
   # return([errornorm(Y22,x,'L2',5,mesh=surface),size]) #yields a strange error
    return([assemble((Y22-x)**2*dx(ref_surface)),size])

for i in range(n):
    v=verificator(N[i],surfaces[i])
    errors[i]=v[0]
    sizes[i]=v[1]

np.savetxt('ver.txt',(np.log(errors),np.log(sizes)))
plt.plot(np.log(sizes),np.log(errors))
plt.savefig("ver.png")