Import function into FEniCS

Question

Hi.

I create my own mesh and know a vector function u0 on every grid point. I need its divergence in the linear form in order to solve the PDE.

How do I properly import this function together with the mesh?
The problem is that I don't know the dof_index of the vertices when I create the mesh. I create a dolfin XML file which loads into FEniCS. The mesh is static.

Thanks a lot :)!

Best,
Roger

Answer 1

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

Answer 2

I think this question has a method in the comments that could work for you.

Comments

Thanks a lot for your comment, you set me up in the right direction. I was able to solve it.