Python

The same procedure can be written in Python[24] in the following:

# test code for libstokes
# Copyright (C) 2006-2008 Kengo Ichiki <kichiki@users.sourceforge.net>
# $Id: manual.tex,v 1.5 2008/10/12 20:16:53 kichiki Exp $

import stokes

sys = stokes.stokes_init()

np = 8
nm = 8
stokes.stokes_set_np(sys, np, nm)

sys.periodic = 1 # periodic boundary condition
lx = 10.0
ly = 10.0
lz = 10.0
stokes.stokes_set_l(sys, lx, ly, lz)

ewald_tr = 60.25
xi = stokes.xi_by_tratio(sys, ewald_tr)
ewald_eps = 1.0e-12
stokes.stokes_set_xi(sys, xi, ewald_eps)

print 'xi =', xi

sys.lubmin2 = 4.0000000001
sys.lubmax = 4.0
stokes.stokes_set_iter(sys, "gmres", 2000, 20, 1.0e-6,
                       1, stokes.get_stdout())

pos = stokes.darray(np*3)
u   = stokes.darray(np*3)
f   = stokes.darray(np*3)

pos[ 0] = 0.0 # x component
pos[ 1] = 0.0 # y component
pos[ 2] = 0.0 # z component

pos[ 3] = 5.0
pos[ 4] = 0.0
pos[ 5] = 0.0

pos[ 6] = 0.0
pos[ 7] = 5.0
pos[ 8] = 0.0

pos[ 9] = 0.0
pos[10] = 0.0
pos[11] = 5.0

pos[12] = 5.0
pos[13] = 5.0
pos[14] = 0.0

pos[15] = 0.0
pos[16] = 5.0
pos[17] = 5.0

pos[18] = 5.0
pos[19] = 0.0
pos[20] = 5.0

pos[21] = 5.0
pos[22] = 5.0
pos[23] = 5.0

for i in range(np*3):
    u[i] = 1.0

print 'pos:'
for i in range(np):
    print i, pos[i*3], pos[i*3+1], pos[i*3+2]

print 'u:'
for i in range(np):
    print i, u[i*3], u[i*3+1], u[i*3+2]

stokes.stokes_set_pos(sys, pos)
stokes.solve_res_3f(sys, u, f)

nc_f = stokes.stokes_nc_init("test-stokes.res-3f.nc",
                             np,
                             0, # nf
                             0, # version
                             0, # flag_poly
                             0, # flag_Q
                             0, # flag_it (time-dependent imposed flow)
                             0, # shear_mode
                             0) # flag_BD
# f0, x, u are active
stokes.stokes_nc_set_f0(nc_f, f)
stokes.stokes_nc_set_time(nc_f, 0, 0.0)
stokes.stokes_nc_set_x(nc_f, 0, pos)
stokes.stokes_nc_set_u(nc_f, 0, u)

stokes.stokes_nc_free(nc_f)


print 'f:'
for i in range(np):
    print i, f[i*3], f[i*3+1], f[i*3+2]