Mixed formulation and compute the error computation

General
1

I am writing code for mixed formulaton for poisson using fenicsx, but I am not getting error computation using exact solution. Here, is the minimal code.
Please suggest something

    from mpi4py import MPI
    from dolfinx import fem, mesh, plot
    import numpy as np
    import ufl
    from petsc4py.PETSc import ScalarType
    from dolfinx.fem.petsc import LinearProblem
    import pyvista as pv
    
    
    domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                     nx=2, ny=2,
                                     cell_type=mesh.CellType.triangle)
    P1 = ufl.FiniteElement(family='BDM', cell=domain.ufl_cell(), degree=1)
    P2 = ufl.FiniteElement(family='DG', cell=domain.ufl_cell(), degree=0)
    W = fem.FunctionSpace(mesh=domain, element=ufl.MixedElement([P1, P2]))
    
    # test and trial functions
    tau, v = ufl.TestFunctions(function_space=W)
    sigma, u = ufl.TrialFunctions(function_space=W)
    
    
    # function determining if a node is on the tray top
    def on_left_and_right(x):
        return np.logical_or(np.isclose(x[0], 0),
                             np.isclose(x[0], 1))
    
    W0, _ = W.sub(0).collapse()
    u_zero = fem.Function(W0)
    u_zero.x.array[:] = 0.0 
    # print(u_zero.x.array[:])
    
    
    boundary_facets = mesh.locate_entities_boundary(mesh=domain,
                                                    dim=domain.topology.dim-1,
                                                    marker=on_left_and_right)
    boundary_dofs = fem.locate_dofs_topological(V=(W.sub(0), W0),
                                                entity_dim=domain.topology.dim-1,
                                                entities=boundary_facets)
    
    # apply Dirichlet BC
    bc = fem.dirichletbc(value=u_zero,
                         dofs=boundary_dofs,
                         V=W.sub(0))  # V.sub(0): Corresponds to BDM element
    
    
    # apply Neumann BC
    x = ufl.SpatialCoordinate(domain=domain)
    f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
    g = ufl.sin(5 * x[0])
    
    # Variational form
    a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
    L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds
    
    # Solve
    problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()
    
     ###### Define Exact solution
    Q, _ = W.sub(0).collapse()
    # dofs_top = fem.locate_dofs_topological((W.sub(1), Q), fdim, facets_top)
    uD = fem.Function(Q)
    uD.interpolate(lambda x: x[0]*( 1 - x[0]) + x[1]*( 1 - x[1])   ) ###### Exact solution
    
    ##### Error_max
    # error_max = np.max(np.abs(uD.x.array-uh.x.array))
    
    L2_error = fem.form(ufl.inner(uh - uD, uh - uD) * ufl.dx)
2 
 from mpi4py import MPI
from dolfinx import fem, mesh, plot
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
from dolfinx.fem.petsc import LinearProblem
import pyvista as pv


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=2, ny=2,
                                 cell_type=mesh.CellType.triangle)
P1 = ufl.FiniteElement(family='BDM', cell=domain.ufl_cell(), degree=1)
P2 = ufl.FiniteElement(family='DG', cell=domain.ufl_cell(), degree=0)
W = fem.FunctionSpace(mesh=domain, element=ufl.MixedElement([P1, P2]))

# test and trial functions
tau, v = ufl.TestFunctions(function_space=W)
sigma, u = ufl.TrialFunctions(function_space=W)


# function determining if a node is on the tray top
def on_left_and_right(x):
    return np.logical_or(np.isclose(x[0], 0),
                         np.isclose(x[0], 1))

W0, _ = W.sub(1).collapse()
u_zero = fem.Function(W0)
u_zero.x.array[:] = 0.0
# u_zero(lambda x: 5.0  ) ###### Exact solution
# print(u_zero.x.array[:])


boundary_facets = mesh.locate_entities_boundary(mesh=domain,
                                                dim=domain.topology.dim-1,
                                                marker=on_left_and_right)
boundary_dofs = fem.locate_dofs_topological(V=(W.sub(1), W0),
                                            entity_dim=domain.topology.dim-1,
                                            entities=boundary_facets)

# apply Dirichlet BC
bc = fem.dirichletbc(value=u_zero,
                     dofs=boundary_dofs,
                     V=W.sub(1))  # V.sub(0): Corresponds to BDM element


# apply Neumann BC
x = ufl.SpatialCoordinate(domain=domain)
f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])

# Variational form
a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

# Solve
problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
ph = problem.solve()
print(ph.x.array)
# Solve
# problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# uh = problem.solve() ##### Computed solution

try:
    w_h = problem.solve()
except PETSc.Error as e:
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, uh = w_h.split()
print(uh.x.array)
print(sigma_h.x.array)

def u_exact(x):
    return np.sin(x[0])*np.sin(x[1])  ###### Random Exact solution

 ###### Define Exact solution
Q, _ = W.sub(1).collapse()
# dofs_top = fem.locate_dofs_topological((W.sub(1), Q), fdim, facets_top)
uD = fem.Function(Q)
# uD.interpolate(lambda x: x[0]*( 1 - x[0]) + x[1]*( 1 - x[1])   ) ###### Random Exact solution
uD.interpolate(u_exact ) ###### Exact solution

##### Error_max
# error_max = np.max(np.abs(uD.x.array-uh.x.array))

L2_error = fem.form(ufl.inner(uh - uD, uh - uD) * ufl.dx)

And the output is

[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf]
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf]
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf]