I am trying to solve the problem \nabla u = f where the solution is u has lower and upper bounds i.e 0 \leq u \leq 1 and \int u\;dx=c_0. I follow steps and example posted in fenics demos. So I used SNES solver for upper and lower bounds and Lengrange Function Space for integral constraint. Integral constraint is working fine, I am having problems with the lower and upper bounds. Any help will be much appreciated
import numpy
from dolfin import *
mesh = UnitSquareMesh(64, 64)
V = FiniteElement("CG", mesh.ufl_cell(), 1)
R = FiniteElement("R", mesh.ufl_cell(), 0)
W_elem = MixedElement([V, R])
W = FunctionSpace(mesh, W_elem)
c0 = Constant(5)
u_lower = Constant(0.0) # lower bound
u_upper = Constant(1.0) # upper bound
# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("-sin(5*x[0])", degree=2)
a = inner(grad(u), grad(v)) * dx + c * v * dx + u * d * dx
L = f * v * dx + c0 * d * dx
F = a - L
# Compute solution
w = Function(W)
bc = []
# J = derivative(F, w, TrialFunction(W))
J = derivative(L, w, TestFunction(W))
H = derivative(L, w, TrialFunction(W))
snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver": {"linear_solver": "lu",
"maximum_iterations": 20,
"report": True,
"error_on_nonconvergence": False}}
problem = NonlinearVariationalProblem(J, w, bc, H)
# handle the bounds
lower = Function(W)
upper = Function(W)
ninfty = Function(W.sub(0).collapse()); ninfty.vector()[:] = -numpy.infty
pinfty = Function(W.sub(0).collapse()); pinfty.vector()[:] = numpy.infty
fa = FunctionAssigner(W, [V])
fa.assign(lower, [interpolate(y_lower, V)])
fa.assign(upper, [interpolate(y_upper, V)])
# set bounds and solve
problem.set_bounds(lower, upper)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
info(solver.parameters, True)
(iter, converged) = solver.solve()
# w = Function(W)
# solve(a == L, w)
(u, c) = w.split()
plot(u, title='Final solution')
filename = 'demo_u.png'
plt.savefig(filename)
plt.close()
# Save solution in VTK format
file = File("output/demo_u.pvd")
file << u```