Assume I have a problem on a manifold. Let this for the time being be the sphere. I can solve a model problem with a known-right hand side on the manifold, but would like to calculate the error in a particular way. If Y_22 is my exact solution, and u_n is my approximate solution defined on a mesh. How do I calculate the error over my sphere?
I would like to “lift” the function u_n to the sphere, by defining u_n^l on the sphere by evaluating in the point on the mesh that is the euclidian projection of the point on the sphere that we would like to evaulate in, and then calculate the L2-norm. How do I go about to do this, or is there a different way that does not include “integrating” on a mesh?
For the specific case of a sphere, you can project a point on the mesh to the sphere by calculating its polar angular coordinates from its Cartesian spatial coordinates. If your exact solution is expressed in terms of these angles (or some other coordinate system on the sphere that you can convert to in closed form), you could do something like the following:
from dolfin import *
from mshr import *
N = 20
mesh = BoundaryMesh(generate_mesh(Sphere(Point(0,0,0),1),N),"exterior")
x = SpatialCoordinate(mesh)
r = sqrt(inner(x,x))
phi = atan(x[1]/x[0])
theta = acos(x[2]/r)
u_exact = cos(phi)*sin(theta)
V = FunctionSpace(mesh,"CG",1)
uh = project(u_exact,V)
e_u = uh - u_exact
import math
print(math.sqrt(assemble(inner(e_u,e_u)*dx)))
Ah, I see, thank you for your answer.
There is no way of performing the integral on the sphere, and not some mesh? The thing is, my theoretical results build on being able to compute the error “on” the sphere, and not a discretization of it.
Given the closed-form mapping \phi:x\mapsto rx/\Vert x\Vert, where x is a point in the mesh and \phi(x) is the corresponding point on the sphere, you can treat the mesh as a parameterization of the sphere. The only thing you need to add to integrate “over the sphere” instead of “over the mesh” is the Jacobian determinant associated with this mapping, to appropriately scale the integration measure (in the same way that you ordinarily integrate “over the physical domain” by including a Jacobian determinant in quadrature rules that are defined to integrate “over the reference element”).
Oh, I got rid of it by upping the number of quadrature nodes. But, Still, running something, like trying to verify convergence of order h on Poissons equation on the sphere yields strange results. Do you see any issues with how I calc the error?
def error_calcer(approx,surface):
x=SpatialCoordinate(surface)
u_exact=sqrt(15)*(x[0]**2-x[1]**2)/sqrt(16*pi) "Y_22 of real spherical harmonics
e_h=u_exact-approx
jacobdet=1-(x[0]+x[1]+x[2])/sqrt(x[0]**2+x[1]**2+x[2]**2)
#jacobdet=1.0 #does no difference!
return(sqrt(assemble(e_h**2/jacobdet*dx,form_compiler_parameters = {"quadrature_degree":10})))
The Jacobian determinant for the surface deformation may be a little more involved than you’re thinking. I didn’t work it out completely myself, but, since the formula for the deformation is given in terms of a 3D ambient space, you’d end up using Nanson’s formula to get the deformed area element. A technical problem here is that the formula for \phi, if considered a volumetric mapping, will have a singular Jacobian, because it maps to a 2D surface; what you really want is some constant extension of \phi along the radial direction, in some neighborhood of the faceted surface. For Nanson’s formula, though, all you need is the Jacobian of this extension evaluated at the surface, which would be D\phi + n\otimes n, where n = x/\Vert x\Vert.
EDIT: Also, for specifying the exact solution in terms of x on the sphere, while treating the mesh as a parameterization of the sphere, you probably meant to compose the mesh’s “x” with the mapping from the mesh to the sphere, e.g.,
x_h = SpatialCoordinate(surface)
x = x_h/sqrt(dot(x_h,x_h)) # Assuming unit sphere
u_exact = ...
EDIT2: Take a look at the following:
from dolfin import *
from mshr import *
N = 20
r = 1.0 # The exact radius of the sphere
mesh = BoundaryMesh(generate_mesh(Sphere(Point(0,0,0),r),N),"exterior")
x_h = SpatialCoordinate(mesh)
r_h = sqrt(inner(x_h,x_h))
n = x_h/r_h # Deformed normal
N = CellNormal(mesh) # Un-deformed normal
# Jacobian of mapping from mesh to exact sphere, guessed from geometrical
# heuristics:
J = ((r/r_h)**2)*abs(inner(N,n))
# Check area integral:
import math
print("Surface area of mesh = "+str(assemble(dot(N,N)*dx)))
print("Exact area of sphere = "+str(4.0*math.pi*(r**2)))
print("Integrated area of sphere = "+str(assemble(J*dx)))
Maybe you can help me with something similar. Imagine I do not have an exact solution, but wish to integrate over two consecutive solutions from two meshes with max size h_1 and h_2. How would I go about to do that? I have tried to run
The following runs, and the error decreases with refinement, although I’m not 100% sure what the interpretation of “extrapolation” is for manifold meshes, and there will be a significant quadrature error, since the quadrature over mesh_fine doesn’t in general conform to mesh_coarse.
from dolfin import *
from mshr import *
# Parameters:
N_coarse = 10
N_fine = 20
r = 1.0 # The exact radius of the sphere
# Coarse and fine meshes, function spaces, and functions:
mesh_coarse = BoundaryMesh(generate_mesh(Sphere(Point(0,0,0),r),N_coarse),
"exterior")
mesh_fine = BoundaryMesh(generate_mesh(Sphere(Point(0,0,0),r),N_fine),
"exterior")
V_coarse = FunctionSpace(mesh_coarse,"CG",1)
V_fine = FunctionSpace(mesh_fine,"CG",1)
# Populate functions with some nontrivial data:
x_coarse = SpatialCoordinate(mesh_coarse)
x_fine = SpatialCoordinate(mesh_fine)
f_coarse = project(x_coarse[0]**2,V_coarse)
f_fine = project(x_fine[0]**2,V_fine)
# Allow f_coarse to be extrapolated to points that do not lie exactly on
# mesh_coarse:
f_coarse.set_allow_extrapolation(True)
# Integrate difference on mesh_fine:
e = f_fine - f_coarse
print(assemble(e*dx(domain=mesh_fine)))
# Sanity check: Value of error integral does go down as N_coarse and
# N_fine are doubled.
I’m guessing a bit, without seeing what’s being passed as arguments, but the error is probably due to one of the calls to project() rather than the integral for e_h. It looks like either f or fc is not defined on the mesh corresponding to the function space cors[1]. (Keep in mind that project() is shorthand for assembling and solving an L^2 projection.) So, if f and/or fc are of type Function, they would need to be set to allow extrapolation.