Hi I am trying to reformulate the mixed poisson equation. The original formulation imposed \sigma \cdot n = g as an essential boundary condition. I wanted to do the opposite, imposing u = u_0 on \Gamma_D as the essential boundary condition. But I couldn’t get the same result. I didn’t know what went wrong. Can you help me on this?

I use DOLFINx 0.6.0 and grabbed some pieces from tutorial to set the boundary conditions.

I am also curious why in the original formulation, we have

```
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1).
```

The order in Q space should always be higher than that in P space. Otherwise no result returned, say not working for k=2 for both.

Here is my code for my formulation:

```
import numpy as np
import pyvista
from ufl import FiniteElement, MixedElement
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import(Measure, SpatialCoordinate, FacetNormal, TestFunctions, TrialFunctions,
div, grad, exp, inner, sin)
from mpi4py import MPI
from petsc4py import PETSc
domain = mesh.create_unit_square(MPI.COMM_WORLD, 12, 12, mesh.CellType.triangle)
k = 3
Q_el = FiniteElement("BDM", domain.ufl_cell(), k-1)
P_el = FiniteElement("DG", domain.ufl_cell(), k)
V_el = MixedElement([Q_el, P_el])
V = fem.FunctionSpace(domain, V_el)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = sin(5 * x[0])
dx = Measure("dx", domain)
boundaries = [(1, lambda x: np.isclose(x[1], 0.0)),
(2, lambda x: np.isclose(x[0], 1.0)),
(3, lambda x: np.isclose(x[1], 1.0)),
(4, lambda x: np.isclose(x[0], 0.0))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
a = (inner(sigma, tau) - inner(grad(u), tau) - inner(sigma, grad(v))) * dx
L = - f * v * dx - v * g * ds(1) - v * g * ds(3)
def on_left_or_right(x):
return np.isclose(x[0], 0) | np.isclose(x[0], 1)
P, _ = V.sub(1).collapse()
u_zero = fem.Function(P)
u_zero.x.array[:] = 0.0
boundary_facets = mesh.locate_entities_boundary(domain, domain.topology.dim-1,
on_left_or_right)
boundary_dofs = fem.locate_dofs_topological((V.sub(1), P), domain.topology.dim-1,
boundary_facets)
bc = fem.dirichletbc(u_zero, boundary_dofs, V.sub(1))
problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "umfpack"})
w_h = problem.solve()
w_h.x.array
```