Mesh points outside the domain

mesh
1

Hi everyone,

I am implementing an algorithm for 2D mesh adaptation. The coordinates x,y of the new mesh are defined through the variational formulation \nabla_{\xi} (w(\xi)x (\xi)) = 0 and \nabla_{\xi} (w(\xi)y(\xi)) = 0, where \xi are the 2D coordinates on a fixed uniform computational mesh. w is a monitor function defined as \sqrt{1 + |\nabla{u(x,y)}|^{2}} and u(x,y) is a scalar field that satisfies a hyperbolic PDE.

SET-UP:

mesh_comp = UnitSquareMesh(N,N)   # Computational mesh (fixed)
mesh = Mesh(mesh_comp)                  # Physical mesh 

V = FunctionSpace(mesh_comp,'CG',1)  # Function space for x,y  
V2 = FunctionSpace(mesh_comp,'CG',2) # Function space for monitor function w

U = FunctionSpace(mesh,'CG',3)  # Function space for u

u0 = Function(U)
x = Function(V)
y = Function(V)
w = Function(V2)

Variational Formulation:

x_D = Expression('x[0]',element = V.ufl_element())   
y_D = Expression('x[1]',element = V.ufl_element())

bc_x = DirichletBC(V,x_D,'on_boundary')
bc_y = DirichletBC(V,y_D,'on_boundary')

a = - inner(grad(x_test),w*grad(x_trial))*dx
L = Constant(0.0)*x_test*dx

solve(a==L,x,bc_x)

a = - inner(grad(y_test),w*grad(y_trial))*dx
L = Constant(0.0)*y_test*dx

solve(a==L,y,bc_y)

Once I have obtained x,y, I want to re-iterate the variational procedure by first interpolating u into the new mesh coordinate, update u according to a PDE (I skip the passage here) and interpolate the monitor function w into the computational mesh:

mesh_ = Mesh(mesh)

mesh_.coordinates()[:,0] = x.compute_vertex_values() 
mesh_.coordinates()[:,1] = y.compute_vertex_values()

U_interp = FunctionSpace(mesh_,'CG',3)

u_interp = Function(U_proj)
u_interp.interpolate(u0)

mesh.coordinates()[:,0] = x.compute_vertex_values()
mesh.coordinates()[:,1] = y.compute_vertex_values()

u0.assign(u_interp)

w_proj = project(sqrt(1 + inner(grad(u0),grad(u0))),U)

w.interpolate(w_proj)

In the last step I obtain the following error:

Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling “Function::set_allow_extrapolation(true)” on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
I have already tried to enforce different boundary conditions and it seems that there is no coordinates outside the squared boundary.

I have already tried to set allow_extrapolation(True), impose Dirichlet/Neumann Boundary conditions for the variational formulation, enforce the coordinates x,y inside the domain with

for index in range(n_vertices):

if near(x.vector()[index],0.0,tol) or x.vector()[index] < 0.0:

    x.vector()[index] = 0.0    

elif near(x.vector()[index],1.0,tol) or x.vector()[index] > 1.0:
    
    x.vector()[index] = 1.0


if near(y.vector()[index],0.0,tol) or y.vector()[index] < 0.0:

    y.vector()[index] = 0.0    

elif near(y.vector()[index],1.0,tol) or y.vector()[index] > 1.0:
    
    y.vector()[index] = 1.0

but I still cannot solve the issue

Thank you very much for your time and help.