Well what a gain in performance! I knew there had to be a much faster way. 10 seconds instead of 1 minute.
The way I solved it now:
from dolfin import *
from numpy import double
from InterpolateCosmo7Wind import *
class U0(Expression):
def eval(self, value, x):
wind = getWindAt(x[0],x[1],x[2])
value[0] = wind[0]
value[1] = wind[1]
value[2] = wind[2]
def value_shape(self):
return(3,)
mesh = Mesh('../cosmo7.xml')
Vv = VectorFunctionSpace(mesh,'CG',1)
V = FunctionSpace(mesh,'CG',1)
u0_dofs = Function(Vv)
u0s_dofs = Function(V)
# stopped timing after got the feeling that it takes > 1 minute
def load_u0_expr():
u0_expr = U0()
print 'writing file...'
File('u0_expr.pvd') << project(u0_expr, Vv)
# runs in under 10 seconds
def load_u0_dofs():
# assume that the DOFs are ordered!
i = 0
# iterate over V.dofmap() even if we set a vectorfunction! len(Vv.dofmap().dofs()) == 3 * len(V.dofmap().dofs())
for [x,y,z] in V.dofmap().tabulate_all_coordinates(mesh).reshape((-1,3)):
wind = getWindAt(x,y,z)
u0_dofs.vector()[i] = np.double(wind[0])
u0_dofs.vector()[i+1] = np.double(wind[1])
u0_dofs.vector()[i+2] = np.double(wind[2])
i += 3
print 'writing file...'
File('u0_dofs.pvd') << u0_dofs
IPython:
In [19]: %run u0.py
In [20]: %timeit load_u0_dofs()
writing file...
writing file...
writing file...
writing file...
1 loops, best of 3: 8.3 s per loop
In [22]: %run u0.py
In [23]: %timeit load_u0_expr()
writing file...
Solving linear system of size 136458 x 136458 (PETSc Krylov solver).
writing file...
Solving linear system of size 136458 x 136458 (PETSc Krylov solver).
writing file...
Solving linear system of size 136458 x 136458 (PETSc Krylov solver).
writing file...
Solving linear system of size 136458 x 136458 (PETSc Krylov solver).
1 loops, best of 3: 58.8 s per loop