73 lines
1.4 KiB
Nim
73 lines
1.4 KiB
Nim
const S = 10
|
|
|
|
type
|
|
|
|
NodeKind = enum nodeFree, nodeA, nodeB
|
|
|
|
Node = object
|
|
v: float
|
|
fixed: NodeKind
|
|
|
|
Mesh[H, W: static int] = array[H, array[W, Node]]
|
|
|
|
|
|
func setBoundary(m: var Mesh) =
|
|
m[1][1].v = 1.0
|
|
m[1][1].fixed = nodeA
|
|
m[6][7].v = -1.0
|
|
m[6][7].fixed = nodeB
|
|
|
|
|
|
func calcDiff[H, W: static int](m,: Mesh[H, W]; d: var Mesh[H, W]): float =
|
|
for i in 0..<H:
|
|
for j in 0..<W:
|
|
var v = 0.0
|
|
var n = 0
|
|
if i > 0:
|
|
v += m[i - 1][j].v
|
|
inc n
|
|
if j > 0:
|
|
v += m[i][j - 1].v
|
|
inc n
|
|
if i + 1 < m.H:
|
|
v += m[i + 1][j].v
|
|
inc n
|
|
if j + 1 < m.W:
|
|
v += m[i][j + 1].v
|
|
inc n
|
|
v = m[i][j].v - v / n.toFloat
|
|
d[i][j].v = v
|
|
if m[i][j].fixed == nodeFree:
|
|
result += v * v
|
|
|
|
|
|
func iter[H, W: static int](m: var Mesh[H, W]): float =
|
|
var
|
|
d: Mesh[H, W]
|
|
cur: array[NodeKind, float]
|
|
diff = 1e10
|
|
|
|
while diff > 1e-24:
|
|
m.setBoundary()
|
|
diff = calcDiff(m, d)
|
|
for i in 0..<H:
|
|
for j in 0..<W:
|
|
m[i][j].v -= d[i][j].v
|
|
|
|
for i in 0..<H:
|
|
for j in 0..<W:
|
|
var k = 0
|
|
if i != 0: inc k
|
|
if j != 0: inc k
|
|
if i < m.H - 1: inc k
|
|
if j < m.W - 1: inc k
|
|
cur[m[i][j].fixed] += d[i][j].v * k.toFloat
|
|
|
|
result = (cur[nodeA] - cur[nodeB]) / 2
|
|
|
|
|
|
when isMainModule:
|
|
|
|
var mesh: Mesh[S, S]
|
|
let r = 2 / mesh.iter()
|
|
echo "R = ", r
|