Laplace equation in a gapped rectangular domain with finite difference method

Question

I have a situation that is shown by this picture:

For this situation I have this code.

Clear["Global`*"]
nx = 40; ny = 20;
Do[h[i, -1] = h[i, 1], {i, 0, nx}]
Do[h[nx + 1, j] = h[nx - 1, j], {j, 0, ny}]
Do[h[i, ny] = 10, {i, 0, nx/2}]
Do[h[i, ny] = 1, {i, nx/2 + 1, nx}]
Do[h[-1, j] = h[1, j], {j, 0, ny}]
Do[eq[i, j] = h[i - 1, j] + h[i + 1, j] + h[i, j - 1] + h[i, j + 1] - 4h[i, j] == 0., {i, 0, nx}, {j, 0, ny}]
sol = Solve[Flatten[Table[eq[i, j], {i, 0, nx}, {j, 0, ny}]]];
Do[h[i, j] = h[i, j] /. sol[[1]], {i, 0, nx}, {j, 0, ny}];
h1 = Interpolation[Flatten[Table[{{x, y}, h[x, y]}, {x, 0, nx}, {y, 0, ny}], 1]];
Plot3D[h1[x, y], {x, 0, nx}, {y, 0, ny}]

hx = D[h1[x, y], x];
hy = D[h1[x, y], y];
StreamPlot[{-hx, -hy}, {x, 0, nx}, {y, 0, ny}]

However, now I have different situation as in this picture

I need to modify the code I just wrote to describe this situation, so on that box ( 6 m by 3 m ) there is no flow

I would appreciate your help.

Answer 1

Here's my FDM-based solution for your problem:

{{xl, xr}, {yl, yr}} = {{0, 33}, {0, 15}};
sf = 2;
nx = sf xr; ny = sf yr;
dx = (xr - xl)/nx; dy = (yr - yl)/ny;
xmidl = 15; xmidr = 18; ymid = 9;
h1 = 14; h2 = 2;
formula = Select[
   Flatten@Table[
      h[x - dx, y] + h[x + dx, y] + h[x, y - dy] + h[x, y + dy] - 4 h[x, y] == 0, {x, xl,
        xr, dx}, {y, yl, yr, dy}][[2 ;; -2, 2 ;; -2]], 
   FreeQ[#, h[x_, y_] /; xmidl < x < xmidr && y > ymid] &];
oneSideD1[most__, "left"] := oneSideD1[most, -1]
oneSideD1[most__, "right"] := oneSideD1[most, 1]
oneSideD1[h_, x_, step_, direction : 1 | -1] := 
 direction ((3 h@x)/(2 step) - (2 h[x - direction step])/step + h[x - 2 direction step]/
    (2 step))
bcxl = Table[oneSideD1[h[#, y] &, xl, dx, "left"] == 0, {y, yl, yr, dy}][[2 ;; -2]];
bcxr = Table[oneSideD1[h[#, y] &, xr, dx, "right"] == 0, {y, yl, yr, dy}][[2 ;; -2]];
bcyl = Table[oneSideD1[h[x, #] &, yl, dy, "left"] == 0, {x, xl, xr, dx}];
bcyr@1 = Table[h[x, yr] == h1, {x, xl, xmidl, dx}];
bcyr@5 = Table[h[x, yr] == h2, {x, xmidr, xr, dx}];
bcyr@3 = Table[
    oneSideD1[h[x, #] &, ymid, dy, "right"] == 0, {x, xmidl, xmidr, dx}][[2 ;; -2]];
bcyr@2 = Table[
    oneSideD1[h[#, y] &, xmidl, dx, "right"] == 0, {y, ymid, yr, dy}][[2 ;; -2]];
bcyr@4 = Table[
    oneSideD1[h[#, y] &, xmidr, dx, "left"] == 0, {y, ymid, yr, dy}][[2 ;; -2]];
set = Flatten@{formula, bcxl, bcxr, bcyl, bcyr /@ Range@5};
var = Union@Cases[set, h[a_, b_], ∞];
{b, mat} = CoefficientArrays[set, var];
sol = LinearSolve[mat, -N@b];
coord = List @@@ var;
ListPointPlot3D[Flatten /@ ({coord, sol}\[Transpose]), PlotRange -> All]

Remark

1. I've used one-sided difference formula i.e. oneSideD1 to discretize the Neumann boundary condtions. For your simple Neumann b.c., it's not a bad idea to handle them with reflection of course, but do notice reflection is hard to extend to more general cases. (If you want to learn more about one-sided formula, start from page 6 of this book. )

2. sf should be an Integer.

3. You can also use

   SetAttributes[h, NHoldAll]    
   sol2 = Solve[N@set, var]; // AbsoluteTiming
   

   to solve the equation set, but this approach is slower. (The speed difference isn't obvious in this simple case though.)