Laplace equation in a gapped rectangular domain with finite difference method

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

• Well, I think this is problem is quite similar to the former, what trouble do you have? – xzczd Jan 1 '17 at 13:07
• I dont know what the necessary modifaction do I need to do on the code to get the new one. I would apprecuate someone's help – Khs Jan 1 '17 at 14:08
• Just figure out the meaning of the selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &]; line in this answer with the help of document, then you'll have an idea about what to do. – xzczd Jan 1 '17 at 14:48
• @xzczd I am sorry to bother you and and I apprecaite your help. Actually I did not get the code you sent . Until now I have not figured a way to do it :( – Khs Jan 1 '17 at 15:18
• Just press F1 and check the document of every function in that code line. If you still don't know how to use document, check this post. I can give an answer, but I'm afraid you won't understand it if you can't even figure out the meaning of the code line above. – xzczd Jan 1 '17 at 15:27

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.)

• (+1) It works. Since the OP is just learning the ropes, let me just point out the valuable tutorial in the documentation on The Numerical Method of Lines. In particular, that page describes FiniteDifferenceDerivative which is a good way to obtain the matrices automatically. – Jens Jan 1 '17 at 18:10
• @Jens Yeah, NDSolveFiniteDifferenceDerivative is powerful. (I dug out the one-sided difference formula by observing NDSolveFiniteDifferenceDerivative[1, #, f /@ #, DifferenceOrder -> 2] &@(d Range@5) actually. ) Sadly there seems to be no simple way to extend its usage to irregular regions. – xzczd Jan 1 '17 at 18:21
• I think there's something wrong with oneSideD1 because the solution looks like it obeys Dirichlet conditions instead of Neumann. – Jens Jan 1 '17 at 18:45
• @Jens You're right. I missed a 2` in the definition. Edited. Thx for pointing out. – xzczd Jan 1 '17 at 18:50
• Yes, now it looks right! – Jens Jan 1 '17 at 18:56