I have to solve a Poisson Equation using FEM in the FEniCS environment. The equation is$$\nabla^2 F=R\ \ \ in\ \ \omega\\\nabla F\cdot \hat{n}=0\ \ \ in\ \ \partial\omega$$(where $\omega$ is a disk and $R$ is a constant) with Dirichlet boundary conditions at circular subdomains( boundary of small circles lying entirely inside $\omega$) within $\omega$. In addition, I must obtain the maximum possible value of $R$ for different arrangement of the circular subdomains. Note that $F$ is bounded on the above by a threshold value, say $F_{\max}$ in $\omega$ and must be positive. How do I write the code in Python to accomplish this. Should I use step size increment for $R$ and then use for-loop?
Maximum value of load in Poisson Equation
finite element;python;fenics