I am badly struck at this particular formulation of a boundary-value problem. Consider that a function's second partial derivative equals a constant throughout a well-defined circular domain. There are two boundary conditions- one on the bound of the domain, where the flux of the function vanishes and another at each small circular domain inside the big circular domain, where the function has a fixed value. i.e there are small circular regions inside the big circle containing them. $$ k_1\nabla^2\cdot f=R\ \ \ , [x,y\in\Omega=\{x^2+y^2\lt1\}]\\q\cdot\hat{n}=0, q=\nabla\cdot f\ \ \ , [x,y\in\partial\Omega]\\f=k_2\ \ \ ,[x,y\in\partial\Gamma_i]\ \ \ ,\Gamma_i=\{(x-x_i)^2+(y-y_i)^2\lt r^2\},\\ \texttt{where each small circle lies entirely inside the unit circle}\ \\0\lt x_i,y_i,r\lt1$$. How can we formulate the above as a poisson equation with well known boundary conditions like Dirichlet, Neumann etc.? Also, how can we numerically find the maximum possible value of $R$ for a given threshold $f$? Any ideas. Thanks beforehand.
The problem arose in a study modeling oxygen diffusion in blood capillaries scattered in a striated muscle fibre bundle. Note that the small circles represent the capillaries and the unit circle represents the muscle fibre bundle.