Mixed formulation of Poisson problem in dolfinx

dolfinx
4

Hi @dokken,
I could resolve the above mentioned error by following your suggestion like this.
I am trying to reproduce the solution to Poisson equation using the mixed formulation as done using legacy FEniCS.
However, the solution is erroneous as it contains only inf values.
Attached is the new code.
Any help is highly appreciated.

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=32, ny=32,
                                 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(V=W0)
u_zero.x.array[:] = 0.0  # ----------------------------------------------------   Changed here
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()