# 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]
