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))))
a = eps*inner(grad(u), grad(v))*dx + v*inner(vec_f,grad(u))*dx
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:
u_fine = solve_adv_diff(128,eps)
for num_elements in num_elements_list:
u_approx = solve_adv_diff(num_elements,eps)
print(f"{num_elements}: {eps}: {error_L2(u_approx, u_fine):.3e}")
