Create a dolfin.Function from the dolfin.Vector containing its coefficient in the finite element basis

Question

Hi all,

I have a dolfin.Vector whose entries represent the coefficients of a Function in the finite element basis expansion.

What is the safest way to create a dolfin.Function using such dolfin.Vector as degree of freedom.

Initially I was simply using

x_as_a_func = dolfin.Function( V, x_as_a_vector)

However this does not seem to work correctly in parallel (at least in FEniCS 1.6).
So I came up with the following walkaround which seems to work both in serial and in parallel:

x_as_a_func = dolfin.Function(V)
x_as_a_func.vector().zero()
x_as_a_func.vector().axpy(1., x_as_a_vector)

Is this second approach safe?

Also this brings up my second question: How shared degree of freedom are handled by dolfin.Function during the assembly procedure?

Thanks in advance!

Comments

Here there is a code to test my issue:

import dolfin as dl

def vector2Function_1(fe_vector, Vh):
    fun = dl.Function(Vh)
    fun.vector().zero()
    fun.vector().axpy(1., fe_vector)
    return fun

def vector2Function_2(fe_vector, Vh):
    return dl.Function(Vh, fe_vector)

def vector2Function(fe_vector, Vh):
    if False:
        return vector2Function_1(fe_vector, Vh)
    else:
        return vector2Function_2(fe_vector, Vh)

if __name__ == "__main__":
    n = 10
    mesh = dl.UnitSquareMesh(n,n)
    Vh = dl.FunctionSpace(mesh, "CG",1)

    mpi_rank = mesh.mpi_comm().rank

    test = dl.TestFunction(Vh)
    trial = dl.TrialFunction(Vh)

    aform = trial*test*dl.dx
    A = dl.assemble(aform)

    x0 = dl.interpolate(dl.Expression("x[0]"), Vh)

    xv = dl.Vector()
    A.init_vector(xv,1)
    xv.axpy(1., x0.vector())

    x1 = vector2Function(xv,Vh)

    r_x0 = dl.assemble(x0*test*dl.dx)
    r_x1 = dl.assemble(x1*test*dl.dx )
    rv = A*xv

    diff_0 = r_x0 - rv
    dnorm_0 = diff_0.norm("l2")

    diff_1 = r_x1 - rv
    dnorm_1 = diff_1.norm("l2")

    if mpi_rank == 0:
        print "|| r_x0 - rv|| = {0:5e}, || r_x1 - rv || = {1:5e}".format(dnorm_0, dnorm_1)

And here the wrong outputs (using vector2Function_2):

$mpirun -n 1 python test.py
|| r_x0 - rv|| = 7.679109e-18, || r_x1 - rv || = 7.679109e-18

mpirun -n 8 python test.py 
Process 0: Number of global vertices: 121
Process 0: Number of global cells: 200
|| r_x0 - rv|| = 6.033116e-18, || r_x1 - rv || =   inf

And here the correct code outputs (using vector2Function_1):

$ mpirun -n 8 python test.py 
Process 0: Number of global vertices: 121
Process 0: Number of global cells: 200
|| r_x0 - rv|| = 6.033116e-18, || r_x1 - rv || = 5.970442e-18