49 lines
1.3 KiB
Text
49 lines
1.3 KiB
Text
code real RlRes=46, RlOut=48;
|
|
def S = 10;
|
|
|
|
proc SetBoundary(MV, MF);
|
|
real MV; int MF;
|
|
[MF(1,1):= 1; MV(1,1):= 1.0;
|
|
MF(6,7):= -1; MV(6,7):= -1.0;
|
|
];
|
|
|
|
func real CalcDiff(MV, MF, DV, W, H);
|
|
real MV; int MF; real DV; int W, H;
|
|
int I, J, N; real V, Total;
|
|
[Total:= 0.0;
|
|
for I:= 0 to H-1 do
|
|
for J:= 0 to W-1 do
|
|
[V:= 0.0; N:= 0;
|
|
if I then [V:= V + MV(I-1,J); N:= N+1];
|
|
if J then [V:= V + MV(I,J-1); N:= N+1];
|
|
if I+1 < H then [V:= V + MV(I+1,J); N:= N+1];
|
|
if J+1 < W then [V:= V + MV(I,J+1); N:= N+1];
|
|
V:= MV(I,J) - V/float(N); DV(I,J):= V;
|
|
if MF(I,J) = 0 then Total:= Total + V*V;
|
|
];
|
|
return Total;
|
|
];
|
|
|
|
func real Iter(MV, MF, W, H);
|
|
real MV; int MF, W, H;
|
|
real DV, Diff, Cur; int I, J;
|
|
[DV:= RlRes(W); for I:= 0 to W-1 do DV(I):= RlRes(H);
|
|
Diff:= 1E10;
|
|
Cur:= [0.0, 0.0, 0.0];
|
|
while Diff > 1E-24 do
|
|
[SetBoundary(MV, MF);
|
|
Diff:= CalcDiff(MV, MF, DV, W, H);
|
|
for I:= 0 to H-1 do
|
|
for J:= 0 to W-1 do
|
|
MV(I,J):= MV(I,J) - DV(I,J);
|
|
];
|
|
for I:= 0 to H-1 do
|
|
for J:= 0 to W-1 do
|
|
Cur(MF(I,J)+1):= Cur(MF(I,J)+1) +
|
|
DV(I,J) * float(-(I>0) - (J>0) - (I<H-1) - (J<W-1));
|
|
\middle=4; side=3; corner=2
|
|
return (Cur(2)-Cur(0))/2.0;
|
|
];
|
|
|
|
real MeshV(S,S); int MeshF(S,S);
|
|
RlOut(0, 2.0 / Iter(MeshV, MeshF, S, S))
|