# Issues Solving Pure Neumann Problems in dolfinx **URL:** \ **Category:** dolfinx\ **Created:** [February 17, 2023, 9:32pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468 "2023-02-17T21:32:50Z")\ **Posts on this page:** 10\ **Page:** 1
**Author:** ![jkrokowski](https://avatars.discourse-cdn.com/v4/letter/j/0ea827/32.png) [@jkrokowski](https://fenicsproject.discourse.group/u/jkrokowski)\ **Post date:** [February 17, 2023, 9:32pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/1 "2023-02-17T21:32:50Z")
Is there any currently recommended way to solve pure neumann problems in Dolfinx? This seems to be a major limitation for me using dolfinx. There are many cases were we would need to enforce that some average condition for the solution (contraining rigid body modes in mechanics, solving the pressure-poisson equation with incompressible fluid flows, etc), but currently I don’t know of a recommended way to solve these types of problems in dolfinx. The normal way (involving the use of LaGrange multipliers) works fine in legacy dolfin: [https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo\_neumann-poisson.py.html](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html) but it is well know that this is not possible in dolfinx because of the Real Type element does not exist: > > > Builds on previous work of @mscroggs. > > Currently brings in tensor related bran…ches from @garth-wells so would need to wait before merging. Does anyone have a work around or a best practice until real type elements are fully implemented? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [February 17, 2023, 11:09pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/2 "2023-02-17T23:09:01Z")
You can remove the constant nullspace with PETSc. See for instance > > ---
**Author:** ![jkrokowski](https://avatars.discourse-cdn.com/v4/letter/j/0ea827/32.png) [@jkrokowski](https://fenicsproject.discourse.group/u/jkrokowski)\ **Post date:** [March 16, 2023, 10:17pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/3 "2023-03-16T22:17:36Z")
Thanks @dokken. I was able to replicate the Singular Poisson problem in dolfinx using your help. The dolfin example is [here.](https://fenicsproject.org/olddocs/dolfin/2018.1.0/python/demos/singular-poisson/demo_singular-poisson.py.html) It might be nice to include this example in the dolfinx tutorial. Either way, I’ll just post the code I was able to get working here (apologies if this already exists elsewhere, it just took me a bit to put the pieces together): ```auto from mpi4py import MPI from dolfinx import mesh,plot from dolfinx.fem import FunctionSpace,Function,form from dolfinx.fem.petsc import create_vector,assemble_matrix,assemble_vector from ufl import inner,TrialFunction,TestFunction,ds,dx,grad,exp,sin,SpatialCoordinate from petsc4py import PETSc import pyvista # Create mesh and define function space domain = mesh.create_unit_square(MPI.COMM_WORLD, 20,20, mesh.CellType.quadrilateral) V = FunctionSpace(domain, ("CG", 1)) u = TrialFunction(V) v = TestFunction(V) x = SpatialCoordinate(domain) f = 10*exp(-(((x[0] - 0.5)**2) + (x[1] - 0.5)**2) / 0.02) g = -sin(5*x[0]) a = inner(grad(u), grad(v))*dx L = f*v*dx + g*v*ds # Assemble system A = assemble_matrix(form(a)) A.assemble() b=create_vector(form(L)) with b.localForm() as b_loc: b_loc.set(0) assemble_vector(b,form(L)) # Solution Function uh = Function(V) # Create Krylov solver solver = PETSc.KSP().create(A.getComm()) solver.setOperators(A) # Create vector that spans the null space nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD) A.setNullSpace(nullspace) # orthogonalize b with respect to the nullspace ensures that # b does not contain any component in the nullspace nullspace.remove(b) # Finally we are able to solve our linear system :: solver.solve(b,uh.vector) #plotting stuff tdim = domain.topology.dim topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim) grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) plotter = pyvista.Plotter() plotter.add_mesh(grid, show_edges=True,opacity=0.25) u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V) u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry) u_grid.point_data["u"] = uh.x.array.real u_grid.set_active_scalars("u") plotter.add_mesh(u_grid.warp_by_scalar("u",factor=.5), show_edges=True) plotter.view_vector((-0.25,-1,0.5)) if not pyvista.OFF_SCREEN: plotter.show() ``` ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 16, 2023, 10:42pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/4 "2023-03-16T22:42:13Z")
If you want to add it to the tutorial, feel free to make a Pull request at: [https://github.com/jorgensd/dolfinx-tutorial/pulls](https://github.com/jorgensd/dolfinx-tutorial/pulls) and I will review it. ---
**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\ **Post date:** [April 30, 2024, 8:48pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/5 "2024-04-30T20:48:17Z")
If I want to change the above code to make the average of u equal to 7 instead of 0, how would I do it? This would help me to understand what is in the code. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 1, 2024, 11:36am UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/6 "2024-05-01T11:36:08Z")
You can simply add the constant that makes the average 7 after the fact, as a Poisson problem with pure Neumann conditions are determined only up to a constant. ---
**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\ **Post date:** [May 1, 2024, 12:49pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/7 "2024-05-01T12:49:21Z")
That’s cheating. This is a baby problem, where we already know the answer. The motivation for this is a somewhat more complicated problem, where we don’t know the answer. The version of this I am trying to solve (another baby problem…) is here > [@Implementing a potential problem with a current condition](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491): > As I understand it, the old fenics had the capability to introduce additional constraints into the PDE, which fenicsx apparently doesn’t allow. Is there some workaround for this which doesn’t rely on already knowing what the solution in the null space looks like? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [May 1, 2024, 4:07pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/8 "2024-05-01T16:07:48Z")
> [@Daniel.Baker](#): > > As I understand it, the old fenics had the capability to introduce additional constraints into the PDE, which fenicsx apparently doesn’t allow. I guess you are talking about a Lagrange multiplier from the real space? There are ways of implementing this in DOLFINx, > > > \### Describe new/missing feature > > It is known that the Poisson equation has infi…nite solutions when the boundary conditions like pure neuman are considered. In the Legacy DOLFIN, we use the Lagrange Multiplier to add extral constraint, which make the system has a unique solution, whose average is 0, see more details in \[demo\](https://fenicsproject.org/olddocs/dolfin/1.4.0/python/demo/documented/neumann-poisson/python/documentation.htmll). > > However, the implementation of Lagrange Multiplier involves the \*\*Real number\*\* function space, which is not avaliable in DOLFINx, so it has been a problem. Now I would like to propose a method to implement the Lagrange Multiplier, this method does not require any new feature, and just use all the functional in current DOLFINX. > > Let us consider the following Variational form of Poisson equation: > !\[Screenshot from 2024-04-09 11-42-13\](https://github.com/FEniCS/dolfinx/assets/98745087/103d33cb-4c6c-423a-94e7-c4a045b452ff) > The LHS of the above VF results the following matrix: > !\[Screenshot from 2024-04-09 11-49-57\](https://github.com/FEniCS/dolfinx/assets/98745087/88ffb021-cf9c-459b-8a9a-11278c52e040) > we can see that, when we apply the Lagrange Multiplier, the resulting matrix is just the matrix given by regular Poisson solver, with one more column and row (marked as red), because the basis function of Real number space is just 1. We can easily obtain the red vector, by defining the following linear form in DOLFINX: > > L = v \* dx > where v is test function in the function space, and then if we assemble the vector based on this linear form, we can obtain the vector that is going to added to the (m+1)th row and (m+1)th column, if the size of the original matrix in m by m. > > I have a demo to show the implementation and performance of the idea. Basically, it is adoped from the Poisson demo in DOLFINX library, but I just remove all the Dirichlet BC, which makes the syetem ill-posed. I can obtain a unique solution with the method mentioned above, and this method can be applied with DOLFINX\_MPC, there is no conflict. All the codes are written in DOLFINX v0.7.3. > > \*\*I believe the best way to implement the above idea is to use the nested operator, since in that case, the solution vector is blocked from the addition unknow value, and we can easily split the solution. I also have some codes using nested matrix and vector, and I can share it as well. But I can not use direct solver for nested operator, and the iteration solvers sometimes produce bad solution (not stable performance). I'm not that familiar with PETSc, but I guess this issue can be resolved by using proper PCtype and KSPtype. Also, the current code just works in series.\*\* > > Any comments or suggests are welcomed. > > > > > \### Suggested user interface > > \`\`\`python3 > \# This demo is modified from this demo https://github.com/FEniCS/dolfinx/blob/v0.7.3/python/demo/demo\_poisson.py, > \# but remove the Dirichlet BC > > \# ## Equation and problem definition > \# > \# For a domain $\\Omega \\subset \\mathbb{R}^n$ with boundary $\\partial > \# \\Omega = \\Gamma\_{D} \\cup \\Gamma\_{N}$, the Poisson equation with > \# particular boundary conditions reads: > \# > \# $$ > \# - \\nabla^{2} u &= f \\quad {\\rm in} \\ \\Omega, \\\\ > \# u &= 0 \\quad {\\rm on} \\ \\Gamma\_{D}, \\\\ > \# \\nabla u \\cdot n &= g \\quad {\\rm on} \\ \\Gamma\_{N}. \\\\ > \# $$ > \# > \# where $f$ and $g$ are input data and $n$ denotes the outward directed > \# boundary normal. The variational problem reads: find $u \\in V$ such > \# that > \# > \# $$ > \# a(u, v) = L(v) \\quad \\forall \\ v \\in V, > \# $$ > \# > \# where $V$ is a suitable function space and > \# > \# $$ > \# a(u, v) &:= \\int\_{\\Omega} \\nabla u \\cdot \\nabla v \\, {\\rm d} x, \\\\ > \# L(v) &:= \\int\_{\\Omega} f v \\, {\\rm d} x + \\int\_{\\Gamma\_{N}} g v \\, {\\rm d} s. > \# $$ > \# > \# The expression $a(u, v)$ is the bilinear form and $L(v)$ > \# is the linear form. > > > from mpi4py import MPI > from petsc4py import PETSc > > > \# + > import numpy as np > > import ufl > from dolfinx import la > from dolfinx import fem, io, mesh, plot > from dolfinx.fem.petsc import create\_matrix, create\_vector, assemble\_vector, assemble\_matrix\_mat, apply\_lifting > from dolfinx.fem.forms import form as \_create\_form > from dolfinx.fem.function import Function > from ufl import ds, dx, grad, inner > > \# - > > \# Note that it is important to first \`from mpi4py import MPI\` to > \# ensure that MPI is correctly initialised. > > \# We create a rectangular {py:class}\`Mesh \\` using > \# {py:func}\`create\_rectangle \\`, and > \# create a finite element {py:class}\`function space > \# \\` $V$ on the mesh. > > \# + > msh = mesh.create\_rectangle(comm=MPI.COMM\_WORLD, > points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16), > cell\_type=mesh.CellType.triangle) > V = fem.functionspace(msh, ("Lagrange", 1)) > \# - > > \# The second argument to {py:func}\`functionspace > \# \\` is a tuple \`(family, degree)\`, where > \# \`family\` is the finite element family, and \`degree\` specifies the > \# polynomial degree. In this case \`V\` is a space of continuous Lagrange > \# finite elements of degree 1. > \# > \# To apply the Dirichlet boundary conditions, we find the mesh facets > \# (entities of topological co-dimension 1) that lie on the boundary > \# $\\Gamma\_D$ using {py:func}\`locate\_entities\_boundary > \# \\`. The function is provided > \# with a 'marker' function that returns \`True\` for points \`x\` on the > \# boundary and \`False\` otherwise. > > > \# + > u = ufl.TrialFunction(V) > v = ufl.TestFunction(V) > x = ufl.SpatialCoordinate(msh) > f = 10 \* ufl.exp(-((x\[0\] - 0.5) \*\* 2 + (x\[1\] - 0.5) \*\* 2) / 0.02) > g = ufl.sin(5 \* x\[0\]) > a = inner(grad(u), grad(v)) \* dx > L = inner(f, v) \* dx + inner(g, v) \* ds > > > \# Create the bilinear and linear forms for the Poisson equation > \# and assemlbe the matrix and vector > a\_form = \_create\_form(a) > A = create\_matrix(a\_form) > L\_form = \_create\_form(L) > b = create\_vector(L\_form) > assemble\_matrix\_mat(A, a\_form) > A.assemble() > b.set(0) > assemble\_vector(b, L\_form) > > > \# Create the linear form vdx, and assemble the vector, this vector will be added to the stiffness matrix of Poisson equation > Lt = \_create\_form(v \* dx) > bt = create\_vector(Lt) > bt.setFromOptions() > bt.set(0) > assemble\_vector(bt, Lt) > > > > \# Create the nested matrix based on the stiffness matrix > \# nested\_A = \[A bt > \# bt' 0\] > m, n = A.getSize() > nested\_A = PETSc.Mat().create() > nested\_A.setSizes(\[m+1, n+1\]) > nested\_A.setFromOptions() > nested\_A.setUp() > idxm = list(range(m)) > idxn = list(range(n)) > nested\_A.setValues(idxm, idxn, A.getValues(idxm, idxn)) > nested\_A.setValues(idxm, \[n\], bt.getValues(list(range(bt.getSize())))) > nested\_A.setValues(\[m\], idxn, bt.getValues(list(range(bt.getSize())))) > nested\_A.setValue(m, n, 0.0) > nested\_A.assemblyBegin() > nested\_A.assemblyEnd() > > > \# Create the nested vector based on the RHS > \# nested\_b = \[b > \# 0\] > nested\_b = PETSc.Vec().create() > nested\_b.setSizes(m+1) > nested\_b.setUp() > nested\_b.setValues(idxm, b.getValues(list(range(b.getSize())))) > nested\_b.setValue(m, 0.0) > nested\_b.assemblyBegin() > nested\_b.assemblyEnd() > nested\_b.setUp() > > > \# Solve the system and otain the solution > uh = Function(V) > solver = PETSc.KSP().create(msh.comm) > solver.setOperators(nested\_A) > solver.setFromOptions() > solver.solve(nested\_b, nested\_b) > > x = la.create\_petsc\_vector\_wrap(uh.x) > x.setValues(idxm, b.getValues(list(range(nested\_b.getSize()-1)))) > > > \# + > with io.XDMFFile(msh.comm, "out\_poisson/poisson.xdmf", "w") as file: > file.write\_mesh(msh) > file.write\_function(uh) > \# - > > \# and displayed using \[pyvista\](https://docs.pyvista.org/). > > \# + > try: > import pyvista > cells, types, x = plot.vtk\_mesh(V) > grid = pyvista.UnstructuredGrid(cells, types, x) > grid.point\_data\["u"\] = uh.x.array.real > grid.set\_active\_scalars("u") > plotter = pyvista.Plotter() > plotter.add\_mesh(grid, show\_edges=True) > warped = grid.warp\_by\_scalar() > plotter.add\_mesh(warped) > if pyvista.OFF\_SCREEN: > pyvista.start\_xvfb(wait=0.1) > plotter.screenshot("uh\_poisson.png") > else: > plotter.show() > except ModuleNotFoundError: > print("'pyvista' is required to visualise the solution") > print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'") > \# - > \`\`\` The problem you are trying to solve seem related to: > > > Thanks a lot for the work on this library, which is really useful. > > I am tryi…ng to use dolfinx\_mpc to impose a boundary to have the same displacement (something like \[this\](https://fenicsproject.discourse.group/t/rigid-plate-boundary-condition/2056 )) > > In particular for a linear elasticity problem with a vector-valued displacement field I want only to have "rigidity" on one of the components. > > I managed to obtain the results with the following method. However, this is very slow and surely not the good way of doing it. > Could you suggest me how to proceed? > > \`\`\` > def create\_rigid\_constraint(mpc, dofs, subspace=None): > dof\_coordinates = mpc.V.tabulate\_dof\_coordinates()\[dofs, :\] > for i in range(len(dof\_coordinates) - 1): > mpc.create\_general\_constraint( > {dof\_coordinates\[i + 1, :\].tobytes(): {dof\_coordinates\[0, :\].tobytes(): 1}}, > subspace, > subspace, > ) > right\_facets = locate\_entities\_boundary(mesh, 1, right\_boundary) > rigid\_dofs = fem.locate\_dofs\_topological(V, 1, right\_facets) > mpc = dolfinx\_mpc.MultiPointConstraint(V) > create\_rigid\_constraint(mpc, rigid\_dofs, 0) > \`\`\` > > Besides, when considering periodic boundary conditions could you consider adding an option to impose it only on one of the components of a vector-valued field. This would be also a possible solution to the problem above, because "rigidity" can easily be written with a particular periodicity map. > > I think that this "rigidity" constraint is very useful in many applications as shown also by the several questions on the forums about that. Maybe it could be added between the built-in methods off \`MultiPointConstraint\`? > > For completeness, I give below also a full code for a simple example to play with > > \`\`\` > import dolfinx.fem as fem > import dolfinx\_mpc.utils > import numpy as np > > from dolfinx.common import Timer, list\_timings, TimingType > from dolfinx.generation import UnitSquareMesh > from dolfinx.io import XDMFFile > from dolfinx.mesh import locate\_entities\_boundary > from mpi4py import MPI > from petsc4py import PETSc > from ufl import (Identity, > SpatialCoordinate, > TestFunction, > TrialFunction, > as\_vector, > dx, > grad, > inner, > sym, > tr, > ) > > > mesh = UnitSquareMesh(MPI.COMM\_WORLD, 50, 50) > V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) > > \# Generate Dirichlet BC on lower boundary (Fixed) > u\_bc = fem.Function(V) > > with u\_bc.vector.localForm() as u\_local: > u\_local.set(0.0) > > > def left\_boundary(x): > return np.isclose(x\[0\], np.finfo(float).eps) > > > left\_facets = locate\_entities\_boundary(mesh, 1, left\_boundary) > topological\_dofs = fem.locate\_dofs\_topological(V, 1, left\_facets) > bc = fem.DirichletBC(u\_bc, topological\_dofs) > bcs = \[bc\] > > \# Define variational problem > u = TrialFunction(V) > v = TestFunction(V) > > \# Elasticity parameters > E = PETSc.ScalarType(1.0) > nu = 0.3 > mu = fem.Constant(mesh, E / (2.0 \* (1.0 + nu))) > lmbda = fem.Constant(mesh, E \* nu / ((1.0 + nu) \* (1.0 - 2.0 \* nu))) > > \# Stress computation > def sigma(v): > return 2.0 \* mu \* sym(grad(v)) + lmbda \* tr(sym(grad(v))) \* Identity(len(v)) > > > x = SpatialCoordinate(mesh) > \# Define variational problem > u = TrialFunction(V) > v = TestFunction(V) > a\_form = inner(sigma(u), grad(v)) \* dx > L\_form = inner(as\_vector((10 \* (x\[1\] - 0.5) \*\* 2, 0)), v) \* dx > > \# Create MPC for imposing u\[1\] constant on the bottom face > def right\_boundary(x): > return np.isclose(x\[0\], 1.0) > > > right\_facets = locate\_entities\_boundary(mesh, 1, right\_boundary) > rigid\_dofs = fem.locate\_dofs\_topological(V, 1, right\_facets) > > > def create\_rigid\_constraint(mpc, dofs, subspace=None): > dof\_coordinates = mpc.V.tabulate\_dof\_coordinates()\[dofs, :\] > for i in range(len(dof\_coordinates) - 1): > mpc.create\_general\_constraint( > {dof\_coordinates\[i + 1, :\].tobytes(): {dof\_coordinates\[0, :\].tobytes(): 1}}, > subspace, > subspace, > ) > > with Timer("~RIGID: Initialize MPC"): > mpc = dolfinx\_mpc.MultiPointConstraint(V) > create\_rigid\_constraint(mpc, rigid\_dofs, 0) > mpc.finalize() > > problem = dolfinx\_mpc.LinearProblem(a\_form, L\_form, mpc, bcs=bcs) > solver = problem.solver > u\_h = problem.solve() > u\_h.name = "u" > outfile = XDMFFile(MPI.COMM\_WORLD, "rigid\_boundary.xdmf", "w") > outfile.write\_mesh(mesh) > outfile.write\_function(u\_h) > > list\_timings(MPI.COMM\_WORLD, \[TimingType.wall\]) > > \`\`\` ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [May 1, 2024, 4:18pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/9 "2024-05-01T16:18:12Z")
I am locking this to preserve the discussion as close as possible to the original question. @Daniel.Baker please answer to @dokken in the other post that you linked. ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [May 1, 2024, 4:18pm UTC](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/10 "2024-05-01T16:18:21Z")