differential equations - Poisson PDE over a irregular region with FDM

The problem modeled by the Poisson PDE is related to the torsion of prismatic beams and I use the Finite Differences Method (FDM). I've managed to solve the equation over a rectangle region with quadratic mesh, here's my code (the dimensions are 0.2 m × 0.2 m):

a = 0.2;
b = 0.2;
n = 4;
m = 4;

deltax = a/n;
deltay = b/m;
G = 10^6;
theta = 0.01745;
phi[0, j_] = 0;
phi[n, j_] = 0;
phi[i_, 0] = 0;
phi[i_, m] = 0;
fi[i_, j_] = -2 G theta;
vars = Flatten[Table[phi[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];

eqns = Flatten[
Table[(phi[i + 1, j] - 2 phi[i, j] + phi[i - 1, j])/
   deltax^2 + (phi[i, j + 1] - 2 phi[i, j] + phi[i, j - 1])/
   deltay^2 == fi[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];
sol = Solve[eqns, vars][[1]]

But I also need to solve for a circular region with radius 0.1m, and I frankly don't have an idea what to do:

The corresponding difference equation is:

$\frac{2}{\Delta x^{2}}\left(\frac{1}{a(1+a)}\phi_{i-1,j}+\frac{1}{1+a}\phi_{i+1,j}-\frac{1}{a}\phi_{i,j}\right)+\frac{2}{\Delta y^{2}}\left(\frac{1}{b(1+b)}\phi_{i,j-1}+\frac{1}{1+b}\phi_{i,j+1}-\frac{1}{b}\phi_{i,j}\right)=-2G\theta$

I know how to solve it with NDSolve:

R = 0.1;
G = 1000000;
θ = 0.01745329252;
Ω = Disk[{0, 0}, {R, R}];
op = Laplacian[u[x, y], {x, y}] +2 G θ;
Subscript[Γ,D] = {DirichletCondition[u[x, y] == 0, True]};
Φ = 
NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y} ∈ Ω];


Plot3D[Φ[x, y], {x, y} ∈ Ω]

I just want to solve the problem with FDM using Cartesian coordinates.

Answer

As mentioned in the comment above, handling irregular region with FDM is cumbersome and frustrating in my view, actually that's where I stopped my self-learning of FDM and turned to finite element method (FEM), which is more suitable for this task, and finally write this. (Have a look at this post for more information. ) Nevertheless, there seems to be no implementation of FDM on irregular region with Mathematica hitherto, so I'd like to have a try.

Well, I believe that as a centuried technique, there must be some legacy code (probably in other programming language) implementing FDM on irregular region, but strangely I never found a piece. (My skill of googling is too bad?) So the following code is completely self-made and quite Mathematica-style (I think).

The method I chose to treat the irregular boundary is the one OP mentioned in his question: extrapolation. Yeah, that graph and the difference equation illustrates extrapolation. If you're unfamiliar with it, a relatively good learning material is Section 3.4 from p71 and Section 6.4 from p199 of Numerical Solution of Partial Differential Equations - An Introduction (Morton K., Mayers D)… OK, let me explain it briefly.

Graph 1

Red points for the unknowns, gray points for those outside the domain, values on green curve is known.

FDM is a technique solving unknown grids with known grids, but as shown in Graph 1, when implementing FDM on an irregular region, grid points are usually not exactly on the boundary where the function value is known, so what to do? The simplest solution is to treat the grids near the boundary as grids on the boundary, which is called zero-order interpolation and discouraged. (Those material I found doesn't mention the reason, I guess it's probably because the error caused by zero-order interpolation is large and not quantitatively controllable.) To achieve the quantitative control on the error caused by FDM, extrapolation is needed. For example, as illustrated in Graph 2 below, we can use the values on B, C, D to fit a parabola to estimate the value on A:

Graph 2

An equivalent explanation for the extrapolation is weighted difference formula, which is also covered in the book linked above. In the rest part of this answer I'll concentrate on programming. Here's the code:

G = 1000000;
θ = Rationalize[0.01745329252, 0];
m = n = 24;
R = 1/10;
dx = 2 R/m; dy = 2 R/n;


r1 = Solve[(x + R)/dx == i, x][[1]];
r2 = Solve[(y + R)/dy == j, y][[1]];
ruleInner = (x^2 + y^2 < R^2 /. r1 /. r2);
ruleOuter = (x^2 + y^2 > R^2 /. r1 /. r2);
areaInner = Table[ruleInner, {i, 0, m}, {j, 0, n}];
areaOuter = Table[ruleOuter, {i, 0, m}, {j, 0, n}];

ArrayPlot /@ Boole[{areaInner, areaOuter}]

I rationalized all the parameters for convenience. areaInner and areaOuter are masks for extracting the desired grid points:

{var, varNot} = 
  Flatten@Pick[Table[f[i, j], {i, 0, m}, {j, 0, n}], #] & /@ {areaInner, areaOuter};

var are grid points inside the circle i.e. the red points in Graph 1, varNot are the grays.

df[dx_, pos_] := (MapAt[# - 1 &, #, pos] - 2 # + MapAt[# + 1 &, #, pos])/dx^2 &
{d2fx, d2fy} = {df[dx, 1] /@ var, df[dy, 2] /@ var};

d2fx and d2fy are $\frac{\partial ^2\phi }{\partial x^2}$ and $\frac{\partial ^2\phi }{\partial y^2}$.

{d2fxOuter, d2fyOuter} = 
  Intersection[varNot, Cases[#, f[__], Infinity]] & /@ {d2fx, d2fy};
{d2fxBorder, d2fyBorder} = 
  Thread[Complement[Cases[#, f[__], Infinity], var, varNot] -> 0] & /@ {d2fx, d2fy};

d2fxOuter and d2fyOuter are grids outside the region but appearing in the equations i.e. As in Graph 2, d2fxBorder and d2fyBorder are those exactly on the boundaries, where As and Bs superpose. (It is relatively rare, actually there're only 4 in this case). Then we use extrapolation to approximate those outside the region:

{pointBorderd2fx, pointBorderd2fy} = 
  With[{ex = #[[1]], pos = #[[2]]}, 
     With[{pos2 = 2 - pos + 1}, {#, 0} & /@ 
       Flatten[Nearest[{i, j}[[pos]] /. 

             Solve[x^2 + y^2 == R^2 /. r1 /. r2, {i, j}[[pos]]] /. {i,
                j}[[pos2]] -> #[[pos2]], #[[pos]]] & /@ ex]]] & /@ 
{{d2fxOuter, 1}, {d2fyOuter, 2}};

pointBorderd2fx and pointBorderd2fy are points on the boundary that will be used in extrapolation i.e. Bs in Graph 2.

{pointInnerd2fx, pointInnerd2fy} = 
  With[{ex = #[[1]], pos = #[[2]]}, 
     Map[{#[[pos]], #} &, 
      Intersection[#, var] & /@ 
       Thread /@ 

        MapAt[{-2, -1, 1, 2} + # &, ex, {All, pos}], {2}]] & /@ 
{{d2fxOuter, 1}, {d2fyOuter, 2}};

pointInnerd2fx and pointInnerd2fy are points inside the domain that will be used in extrapolation i.e. Cs and Ds in Graph 2.

extra[{x4_, y4_}, {x1_, y1_}, {{x2_, y2_}, {x3_, y3_}}] = 
  y4 -> InterpolatingPolynomial[{{x1, y1}, {x2, y2}, {x3, y3}}, x4];

extra provides the extrapolating formula. The last step is to form the equation set with difference equations, eliminating As with Bs, Cs and Ds and then solve it:

sol = First@
    NSolve[(d2fx /. d2fxBorder /. MapThread[extra, {{#[[1]], #} & /@ d2fxOuter, 

                                                   pointBorderd2fx, 
                                                   pointInnerd2fx}]) + 
           (d2fy /. d2fyBorder /. MapThread[extra, {{#[[2]], #} & /@ d2fyOuter, 
                                                   pointBorderd2fy, 
                                                   pointInnerd2fy}]) == 
           -2 G θ // Thread, var]; // AbsoluteTiming

Finally let's draw a picture for the solution. Notice that we need to add points on the boundary because they are almost not included in sol:

bound = Flatten /@ Transpose@{#, ConstantArray[0., Length@#]} &@
   Cases[Normal@

      ContourPlot[
       x^2 + y^2 == R^2 /. r1 /. r2 // Evaluate, {i, 0, m}, {j, 0, n}], 
         Line[a__] -> a, ∞][[1]];

ListPlot3D[Join[bound, sol /. (f[a_, b_] -> c_) :> {a, b, c}]]