# Attribute error while calculating L2 error for non-matching meshes

Hello, I am facing an issue while calculating L2 error between non-matching mesh interpolation data. It is an attribute error, it was working fine previously. I don’t know whether there is an update to this or there is something wrong in my code.
I would be grateful if you can provide any insight. Please see the full code and error shown below…

``````import numpy as np
import matplotlib.pyplot as plt

from dolfinx.fem import (Constant, Function, functionspace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx import fem, io, mesh
from dolfinx.mesh import locate_entities_boundary
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad, inner
from petsc4py.PETSc import ScalarType

def solve_adv_diff(N=8, eps= 0.00025, degree=1):

# define mesh and function space
msh = create_unit_square(MPI.COMM_WORLD,N,N)
V = functionspace(msh, ("Lagrange", degree))
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(msh)

# define boundary conditions

# left facet
facet_2 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.isclose(x[0], 0))
dof_2 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_2)
bc_2 = dirichletbc(ScalarType(1), dof_2,V)

# bottom left_facet
facet_1 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_and(np.isclose(x[1],0),x[0]<0.5 -(eps/2)))
dof_1 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_1)
bc_1 = dirichletbc(ScalarType(1), dof_1,V)

# bottom right_facet
facet_3 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_and(np.isclose(x[1],0),x[0]>0.5 + (eps/2)))
dof_3 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_3)
bc_3 = dirichletbc(ScalarType(0), dof_3,V)

# linear region with length eps
a = ScalarType(-1/eps)
b = ScalarType((0.5/eps) + 0.5)
uD = Function(V)
uD.interpolate(lambda x : a*x[0] + b )

facet_lin = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_and(np.isclose(x[1],0),np.logical_and(x[0]>= 0.5 - (eps/2),x[0]<= 0.5 + (eps/2))))
dof_4 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_lin)
bc_4 = dirichletbc(uD, dof_4)

# collecting boundary conditions in a list
bcs = [bc_1,bc_2, bc_3,bc_4]

# define variational form
vec_f = Constant(msh,ScalarType((np.cos(np.pi/3),np.sin(np.pi/3))))
f = Constant(msh, ScalarType(0))
L = f * v * dx

# solving problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
return problem.solve()

def error_L2(uh, u_ex, degree_raise = 3):
# Create higher order function space
tol = 1e-9
degree = u_ex.function_space.ufl_element().degree
family = "CG"
msh1 = u_ex.function_space.mesh
W = functionspace(msh1, (family, degree + degree_raise))
u_fine = Function(W)
u_fine.interpolate(u_ex)
u_coarse = Function(W)
# Interpolate approx. solution in fine solution function space

if uh.function_space.mesh != msh1:
from dolfinx.fem import create_nonmatching_meshes_interpolation_data
interpolation_data = create_nonmatching_meshes_interpolation_data(
u_coarse.function_space.mesh._cpp_object,
u_coarse.function_space.element,
uh.function_space.mesh._cpp_object,tol)
u_coarse.interpolate(uh, nmm_interpolation_data=interpolation_data)
else:
u_coarse.interpolate(uh)

# Compute the error in fine solution function space
e_W = Function(W)
e_W.x.array[:] = u_fine.x.array - u_coarse.x.array

# Integrate the error
error = form(ufl.inner(e_W, e_W) * ufl.dx)
error_local = assemble_scalar(error)
error_global = msh1.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)

num_elements_list = [8,16,32,64]
eps_list = [1,0.1,0.01]

print("Compare to finely resolved approximation")

for eps in eps_list:
for num_elements in num_elements_list:
print(f"{num_elements}: {eps}: {error_L2(u_approx, u_fine):.3e}")
``````

See:

For the relevant changes. Note that it just means that you can remove `._cpp_object` from the two mesh arguments.

In general, if you are following the main branch, I would follow:
https://github.com/FEniCS/dolfinx/pulls?q=is%3Apr+is%3Aclosed+is%3Amerged
as it gives you the list of changes.

1 Like

Thank you Dr. @dokken for your reply. I did the changes mentioned but it is still throwing the following error. Please see the updated piece of code and error below. I am only including the updated piece of code to remove the redundancy.

Ah, i see, there is a typo in the Else clause:

It should be
` return PointOwnershipData( *_create_nonmatching_meshes_interpolation_data( mesh_to._cpp_object, element, mesh_from._cpp_object, cells, padding )`

Feel free to make a pr fixing this.
I am not at a computer for 6-8 hours, so i cant make the PR myself.

Btw you sould set `tol` as a kwarg, i.e.

``````interpolation_data = create_nonmatching_meshes_interpolation_data(
u_coarse.function_space.mesh,
u_coarse.function_space.element,
uh.function_space.mesh, tol=tol)
``````

as the current function has two kwargs, `cells` and `tol`, and currently `cells` is set to `tol`.

I’ve made a PR fixing the first part of your issue (Missing C++ accessor to mesh object in create nonmatching interpolation data by jorgensd · Pull Request #3078 · FEniCS/dolfinx · GitHub), but you still need to use he signature above after the PR is merged.

1 Like

Thank you Dr. @dokken, I appreciate your positive response. I will use the correct signature after the PR is merged.

Dr. @dokken I tried to implement the updated code but somehow it’s not working. I would appreciate if you can provide any insight. I tried to put tol as kw arg as well but still it didnt work. Please see the code and error below:

``````u_ex = solve_adv_diff(128,0.00025)
tol = 1e-9
degree = 1
N1 = 8
# define mesh and function space
msh1 = create_unit_square(MPI.COMM_WORLD,N1,N1,cell_type=mesh.CellType.triangle)
V = functionspace(msh1, ("Lagrange", degree))
u_fine = Function(V)
# Interpolate approx. solution in coarse function space
from dolfinx.fem import create_nonmatching_meshes_interpolation_data
interpolation_data = create_nonmatching_meshes_interpolation_data(
u_fine.function_space.mesh,
u_fine.function_space.element,
u_ex.function_space.mesh,msh1.topology.cell_indices,tol)
u_fine.interpolate(u_ex, nmm_interpolation_data=interpolation_data)
``````

Not sure where you got this from, as this is not a function of the topology.

What are you trying to achieve by supplying it?

I was trying to access cell indices for the argument before tol. I think I am doing it wrong. Can you please provide help with it.

Why do you need it. If you do not supply it, it sends in all cells on `msh1`.
If you really want to send something in there, it would be:

``````fine_mesh_cell_map =msh1.topology.index_map(msh1.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
``````

and pass `cells` into `create_nonmatching_meshes_interpolation_data`.

1 Like

I tried without sending it in and used tol as kw arg but it wasn’t working that’s why I thought to send in the cells of the destination mesh. I didn’t know how to access the cells of the the mesh where I have to interpolate. I appreciate your help.

It is still throwing error.

``````u_ex = solve_adv_diff(128,0.00025)
tol = 1e-9
degree = 1
N1 = 8
# define mesh and function space
msh1 = create_unit_square(MPI.COMM_WORLD,N1,N1,cell_type=mesh.CellType.triangle)
V = functionspace(msh1, ("Lagrange", degree))
fine_mesh_cell_map =msh1.topology.index_map(msh1.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
u_fine = Function(V)
# Interpolate approx. solution in coarse function space
meshtags = dolfinx.mesh.MeshTags(msh1)
from dolfinx.fem import create_nonmatching_meshes_interpolation_data
interpolation_data = create_nonmatching_meshes_interpolation_data(
u_fine.function_space.mesh,
u_fine.function_space.element,
u_ex.function_space.mesh,cells,tol)
u_fine.interpolate(u_ex, nmm_interpolation_data=interpolation_data)
``````

Please make a minimal reproducible example.
The current code that you are using isn’t reproducible, as you are missing all kinds of definitions and imports, have introduced unused variables etc.

Dr. @dokken Please see the minimal code example below:

``````import numpy as np
import matplotlib.pyplot as plt

from dolfinx.fem import (Constant, Function, functionspace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical,assemble_matrix)
from dolfinx.fem.petsc import LinearProblem
from dolfinx import fem, io, mesh
from dolfinx.mesh import locate_entities_boundary
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad, inner
from petsc4py.PETSc import ScalarType
def solve_adv_diff(N=8, eps= 0.00025, degree=1):

# define mesh and function space
msh = create_unit_square(MPI.COMM_WORLD,N,N)
V = functionspace(msh, ("Lagrange", degree))
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(msh)

# define boundary conditions

# left facet
facet_2 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.isclose(x[0], 0))
dof_2 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_2)
bc_2 = dirichletbc(ScalarType(1), dof_2,V)

# bottom left_facet
facet_1 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_and(np.isclose(x[1],0),x[0]<0.5 -(eps/2)))
dof_1 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_1)
bc_1 = dirichletbc(ScalarType(1), dof_1,V)

# bottom right_facet
facet_3 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_and(np.isclose(x[1],0),x[0]>0.5 + (eps/2)))
dof_3 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_3)
bc_3 = dirichletbc(ScalarType(0), dof_3,V)

# linear region with length eps
a = ScalarType(-1/eps)
b = ScalarType((0.5/eps) + 0.5)
uD = Function(V)
uD.interpolate(lambda x : a*x[0] + b )

facet_lin = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
marker=lambda x: np.logical_and(np.isclose(x[1],0),np.logical_and(x[0]>= 0.5 - (eps/2),x[0]<= 0.5 + (eps/2))))
dof_4 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_lin)
bc_4 = dirichletbc(uD, dof_4)

# collecting boundary conditions in a list
bcs = [bc_1,bc_2, bc_3,bc_4]

# define variational form
vec_f = Constant(msh,ScalarType((np.cos(np.pi/3),np.sin(np.pi/3))))
f = Constant(msh, ScalarType(0))
L = f * v * dx

# solving problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
return problem.solve()
# interpolating
tol = 1e-9
degree = 1
N1 = 8
# define mesh and function space
msh1 = create_unit_square(MPI.COMM_WORLD,N1,N1,cell_type=mesh.CellType.triangle)
V = functionspace(msh1, ("Lagrange", degree))
fine_mesh_cell_map =msh1.topology.index_map(msh1.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
u_fine = Function(V)
# Interpolate approx. solution in coarse function space
from dolfinx.fem import create_nonmatching_meshes_interpolation_data
interpolation_data = create_nonmatching_meshes_interpolation_data(
u_fine.function_space.mesh,
u_fine.function_space.element,
u_ex.function_space.mesh,cells,tol)
u_fine.interpolate(u_ex, nmm_interpolation_data=interpolation_data)
``````

Suggestion for the (immediate) future: stop posting screenshots, just copy their content in “```”. In taking the screenshot you missed the final line of the error, which may have been helpful in understanding what is going on. Furthermore, screenshots are not easily searchable, which may mean (for you) that people are less likely to help you and (for the rest of the community) that people may not be able to find your error message using the search button.

Please note that your code is far from minimal.
Here is a minimal code that reproduces your error:

``````from dolfinx.fem import create_nonmatching_meshes_interpolation_data
import numpy as np

from dolfinx.fem import Function, functionspace

from dolfinx import mesh
from dolfinx.mesh import create_unit_square
from mpi4py import MPI

N = 128
degree = 1
msh = create_unit_square(MPI.COMM_WORLD, N, N)
V = functionspace(msh, ("Lagrange", degree))
u_ex = Function(V)

tol = 1e-9
N1 = 8
# define mesh and function space
msh1 = create_unit_square(MPI.COMM_WORLD, N1, N1,
cell_type=mesh.CellType.triangle)
V = functionspace(msh1, ("Lagrange", degree))
fine_mesh_cell_map = msh1.topology.index_map(msh1.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
u_fine = Function(V)
# Interpolate approx. solution in coarse function space
interpolation_data = create_nonmatching_meshes_interpolation_data(
u_fine.function_space.mesh,
u_fine.function_space.element,
u_ex.function_space.mesh, cells, tol)
u_fine.interpolate(u_ex, nmm_interpolation_data=interpolation_data)
``````

Turns out my last fix wasn’t a fix at all.
Will be addressed in Revert change to Python interface of non-matching interpolation data by jorgensd · Pull Request #3088 · FEniCS/dolfinx · GitHub
with the fix:

``````interpolation_data = create_nonmatching_meshes_interpolation_data(
u_fine.function_space.mesh.geometry,
u_fine.function_space.element,
u_ex.function_space.mesh, cells, tol)
``````

Dr. @francesco-ballarin Thank you for your suggestion. I will make sure not to post any screenshots in future.

Thank you Dr. @dokken for your kind help. I am truly grateful for your resposne. I hope it will be fixed soon.