At right-hand-side of my Poisson equation I have expression parametrized with two parameters: `A`

and `B`

. For a given `A`

I need to select such `B`

that my equation is compatible with Neumann BC.

I think the simplest way to do so is to integrate the expression over the spatial domain with `A, B = 0, 1`

and `A, B = 1, 0`

and calculate ratio of the results. How may I perform the integration with FEniCS?

A toy example of my code is below (in the real example I have a 3D mesh of complex geometry):

```
from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('A * exp(-0.5 * ((x[0] - X)*(x[0] - X) + (x[1] - Y)*(x[1] - Y) / S2)) - B', degree=1, X=0.5, Y=0.5, S2=0.04, A=1, B=0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solver = KrylovSolver("cg", "ilu")
solver.parameters["maximum_iterations"] = 20000
solver.parameters["absolute_tolerance"] = 1E-8
solver.parameters["monitor_convergence"] = True
solver.solve(assemble(a), u.vector(), assemble(L))
```