Hi!

I am using **dolfinx** version **0.6.0**. It was installed with **conda-forge**. I am also using **multiphenicsx**. My aim is to solve a complicated set of linear elasticity equations with Robin boundary conditions using **multiphenicsx** and block forms. Also, I eventually want to impose additional constrains with Lagrange multipliers.

As a starting point, I am trying to solve Poisson’s equation with Pure Neumann boundary conditions. To do so, I am trying to translate this code using **multiphenicsx** and block forms. Below is a minimal working code example that shows how I am trying to do it. The output that this code produces shows that I am doing something terribly wrong.

How should I write Poisson’s equation with Pure Neumann boundary conditions and Lagrange multipliers in block form? Why am I getting such wrong output? What should I fix in the code below to solve the problem correctly?

```
import numpy as np
from mpi4py import MPI
import dolfinx.fem
import dolfinx.fem.petsc
from dolfinx import mesh, plot
from petsc4py import PETSc
import ufl
import pyvista
import multiphenicsx.fem
import multiphenicsx.fem.petsc
def solve_poisson_demo():
# MESH
# mpi vars
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
# rectangular mesh and facets
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
cell_type=mesh.CellType.triangle,)
facets = mesh.locate_entities_boundary(domain, dim=1,
marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 2.0)))
# FUNCTION SPACES / BLOCK FUNCTION SPACE
V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2))
R = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2))
(u, c) = (ufl.TrialFunction(V), ufl.TrialFunction(R))
(v, d) = (ufl.TestFunction(V), ufl.TestFunction(R))
# IMPORTANT EXPRESSIONS
f = dolfinx.fem.Function(V)
f.interpolate(lambda x: 10 * np.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02))
g = dolfinx.fem.Function(V)
g.interpolate(lambda x: np.sin(5 * x[0]))
_zero = dolfinx.fem.Function(V)
_zero.interpolate(lambda x: np.zeros(len(x[0])))
dx = ufl.Measure("dx",domain=domain)
ds = ufl.Measure("ds",domain=domain)
# LINEAR AND BILINEAR FORMS
aa = [[ufl.inner(ufl.grad(u),ufl.grad(v))*dx,v*c*dx],
[u*d*dx,None]]
ff = [ufl.inner(f,v)*dx + ufl.inner(g,v)*ds,_zero*d*dx]
aa_cpp = dolfinx.fem.form(aa)
ff_cpp = dolfinx.fem.form(ff)
# ASSEMBLE BLOCK EQUATIONS
AA = dolfinx.fem.petsc.assemble_matrix_block(aa_cpp, bcs=[])
AA.assemble()
FF = dolfinx.fem.petsc.assemble_vector_block(ff_cpp, aa_cpp, bcs=[])
# SOLVE
uu = dolfinx.fem.petsc.create_vector_block(ff_cpp)
# solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(AA)
solver.setType("preonly")
solver.getPC().setType("lu")
solver.getPC().setFactorSolverType("superlu_dist")
# solve
solver.solve(FF,uu)
uu.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
solver.destroy()
# SPLIT
# Split the block solution in components
(uh,c) = (dolfinx.fem.Function(V), dolfinx.fem.Function(R))
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(uu, [V.dofmap, R.dofmap]) as
uu_wrapper:
for u_wrapper_local, component in zip(uu_wrapper, (uh,c)):
with component.vector.localForm() as component_local:
component_local[:] = u_wrapper_local
# POST-PROCESS
if comm.Get_size() == 1:
# plot displacement
topology, cell_types, geometry = plot.create_vtk_mesh(V)
u = uh.x.array
grid = pyvista.UnstructuredGrid(topology,cell_types,geometry)
grid.point_data["u (cm)"] = u
warped = grid.warp_by_scalar("u (cm)", factor=1.0)
plotter = pyvista.Plotter()
plotter.add_mesh(warped, show_scalar_bar=True, show_edges=False, scalars="u (cm)")
plotter.show(screenshot='poisson_demo_visualize.png')
plotter.clear()
return None
if __name__=='__main__':
solve_poisson_demo()
```