Error of functions defined on two different measures

mesh
1

Hi!

I have a sequence of meshings of the sphere, and wish to calculate the error between two consecutive solutions. Find below my attempt:

from dolfin import * 
from mshr import * 
import numpy as np                                                                                                                                                                  
x_h=SpatialCoordinate(coarse_mesh)
r_h=sqrt(inner(x_h,x_h))
n=x_h/r_h
N=CellNormal(coarse_mesh)
J=(1.0/r_h)**2*abs(inner(N,n))
x=x_h/sqrt(dot(x_h,x_h))
f_fine=project(solution_fine,coarse_functionspace)
f_cors=project(solution_course,coarse_functionspace)
f_cors.set_allow_extrapolation(True)
f_fine.set_allow_extrapolation(True)
e_h=f_fine-f_cors
return(sqrt(assemble(e_h**2*J*dx(domain=coarse),form_compiler_parameters = {"quadrature_degree": 5})))

here coarse_mesh is the coarser mesh, coarse_functionspace is a functionspace over the mesh of the same type as fine_functionspace, and solution_fine and solution_course are two solutions from each mesh. It works for the coarsest two meshes but fails for the finer… I get the error:

Error:   Unable to evaluate function at point.
The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)"  on this Function to allow extrapolation.

Changing around which domain I integrate over and such yields other errors, for instance changing the integration domain to the finer mesh and extrapolating the coarse solution gets me

Error: Unable to create mesh entity.
Reason:  Mesh entity index -1 out of range [0, 320] for entity of dimension 2.

I also get a lot of “Applying str() to a metadata value of type set, don’t know if this is safe.”, how do I get rid of that?

So, any suggestions?

2

The problem comes from the call to project(). It can be avoided in your example as follows:

from dolfin import * 
from mshr import * 
import numpy as np

N_coarse = 10
N_fine = 20
r = 1.0
coarse_mesh = BoundaryMesh(generate_mesh(Sphere(Point(0,0,0),r),N_coarse),
                           "exterior")
fine_mesh = BoundaryMesh(generate_mesh(Sphere(Point(0,0,0),r),N_fine),
                         "exterior")
coarse_functionspace = FunctionSpace(coarse_mesh,"CG",1)
fine_functionspace = FunctionSpace(fine_mesh,"CG",1)

#######

# Define some placeholder functions
solution_fine = Function(fine_functionspace)
solution_coarse = Function(coarse_functionspace)

# Because the projection is done onto the coarse function space, anything
# defined on the fine mesh needs to be able to be extrapolated.
solution_fine.set_allow_extrapolation(True)

#######

x_h=SpatialCoordinate(coarse_mesh)
r_h=sqrt(inner(x_h,x_h))
n=x_h/r_h
N=CellNormal(coarse_mesh)
J=(1.0/r_h)**2*abs(inner(N,n))
x=x_h/sqrt(dot(x_h,x_h))
f_fine=project(solution_fine,coarse_functionspace)
f_cors=project(solution_coarse,coarse_functionspace)
#f_cors.set_allow_extrapolation(True) # (Harmless but unnecessary)
f_fine.set_allow_extrapolation(True)
e_h=f_fine-f_cors

print(sqrt(assemble(e_h**2*J*dx(domain=coarse_mesh),
                    form_compiler_parameters = {"quadrature_degree": 5})))
2 Likes
3

Thank you, once more, for your help, but your solution does not work. It still yields the error

Error: Unable to create mesh entity.
Mesh entity index -1 out of range [0, 320] for entity of dimension 2.
Where:   This error was encountered inside MeshEntity.cpp.