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','x','x'),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*x-x*x)/(16*sqrt(pi))",degree=2) #error Y22=Expression("sqrt(15)*(x*x-x*x)/(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 sizes[i]=v np.savetxt('ver.txt',(np.log(errors),np.log(sizes))) plt.plot(np.log(sizes),np.log(errors)) plt.savefig("ver.png")