FEniCS Project

Solution array (u.vector().array()) does not correspond to mesh vertices

Asked by Toby Meierbachtol

Hello,

Calling the solution as an array seems to scramble the results. For instance, printing the vertex coordinates and corresponding solutions in the poisson tutorial yields an erroneous solution. Does calling the expansion coefficients of the solution function indeed result in a scrambled array? If so, is there a work-around?

Thanks for your help.

Toby

Here is my code:

from pylab import *
from dolfin import *

#Create meash and define function space
mesh = UnitSquareMesh(6, 4)

V = FunctionSpace(mesh, 'Lagrange', 1)

#Define boundary conditions
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')

def u0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, u0_boundary)

#Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

#Compute solution
u = Function(V)
solve(a == L, u, bc, solver_parameters={'linear_solver': 'cg', 'preconditioner': 'ilu'})

#Grab all solution values 'u'
u_nodal_values = u.vector()
u_array = u_nodal_values.array()

#Get mesh information
coords = mesh.coordinates()
num_cells = mesh.num_cells()
num_vertices = mesh.num_vertices()
print mesh #Print mesh info

#print out solution and coordinates
if num_vertices == len(u_array):
    for i in range(num_vertices):
        print 'u(%8g,%8g) = %g' % (coor[i][0], coor[i][1], u_array[i])

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Best Johan Hake (johan-hake) said : 		#1

Since dolfin 1.1 there are no guaranties for dof numbering. They are not
following the vertex numberings.

Check out:

  V.dofmap().vertex_to_dof_map(mesh)

Johan

On 03/20/2013 06:21 PM, Toby Meierbachtol wrote:
> New question #224748 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/224748
>
> Hello,
>
> Calling the solution as an array seems to scramble the results. For instance, printing the vertex coordinates and corresponding solutions in the poisson tutorial yields an erroneous solution. Does calling the expansion coefficients of the solution function indeed result in a scrambled array? If so, is there a work-around?
>
> Thanks for your help.
>
> Toby
>
> Here is my code:
>
> from pylab import *
> from dolfin import *
>
>
> #Create meash and define function space
> mesh = UnitSquareMesh(6, 4)
>
> V = FunctionSpace(mesh, 'Lagrange', 1)
>
> #Define boundary conditions
> u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')
>
> def u0_boundary(x, on_boundary):
> return on_boundary
>
> bc = DirichletBC(V, u0, u0_boundary)
>
> #Define variational problem
> u = TrialFunction(V)
> v = TestFunction(V)
> f = Constant(-6.0)
> a = inner(nabla_grad(u), nabla_grad(v))*dx
> L = f*v*dx
>
> #Compute solution
> u = Function(V)
> solve(a == L, u, bc, solver_parameters={'linear_solver': 'cg', 'preconditioner': 'ilu'})
>
> #Grab all solution values 'u'
> u_nodal_values = u.vector()
> u_array = u_nodal_values.array()
>
> #Get mesh information
> coords = mesh.coordinates()
> num_cells = mesh.num_cells()
> num_vertices = mesh.num_vertices()
> print mesh #Print mesh info
>
> #print out solution and coordinates
> if num_vertices == len(u_array):
> for i in range(num_vertices):
> print 'u(%8g,%8g) = %g' % (coor[i][0], coor[i][1], u_array[i])
>
>
>

Revision history for this message
Toby Meierbachtol (toby-meierbachtol) said : 		#2

I was unable to get the vertex_to_dof_map(mesh) command to work. However, I did find another solution to the problem:

u.compute_vertex_values()

Thanks for the quick response!

Revision history for this message
Henrik Garde (1etr3k) said : 		#3

I can not ge the 'vertex_to_dof_map' working either, and after updating to dolfin 1.1.0 a lot of my code has been broken due to this change. Is there some way to revert to dolfin 1.0.0?

Revision history for this message
Johan Hake (johan-hake) said : 		#4

On 03/21/2013 09:11 AM, Henrik Garde wrote:
> Question #224748 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/224748
>
> Henrik Garde posted a new comment:
> I can not ge the 'vertex_to_dof_map' working either,

If you want us to help you, you need to tell us what is not working.

> and after updating
> to dolfin 1.1.0 a lot of my code has been broken due to this change. Is
> there some way to revert to dolfin 1.0.0?

Sure, just remove the present installation and install the older one.

Johan

Revision history for this message
Maximilian Albert (cilix) said : 		#5

Hi Henrik,

some instructions for downgrading dolfin on Ubuntu systems can be found in the thread at [1] (see e.g. comment #4). There is more information in other threads which I can dig up if you're interested.

However, a much easier solution with dolfin-1.1 is the following:

    import dolfin as df
    df.parameters.reorder_dofs_serial = False

This switches off the re-ordering of the dofs so that your code will continue to work (at least in serial).

Furthermore, as far as I'm aware, the 'vertex_to_dof_map' isn't included in any released versions yet, so I guess this is why you couldn't get it to work (am I correct to presume that the command doesn't exist when you tried it?). However, the developers told me that a new release is due any day now so your problems should be fixed very soon. (Incidentally, we had the same problem with our code and are also eagerly awaiting this new release. ;-)) For a nice, concise recipe how to use the new 'vertex_to_dof_map' functionality once it becomes available, see Johan's answer #2 in [2].

Hope this helps!

[1] https://answers.launchpad.net/fenics/+question/220534
[2] https://answers.launchpad.net/dolfin/+question/222509

Revision history for this message
Henrik Garde (1etr3k) said : 		#6

Sorry Johan I should have been more precise in my statement. However, Maximilian guessed correctly that 'vertex_to_dof_map' was not included in the version I have. The trick with 'parameters.reorder_dofs_serial = False' did solve the ordering, so thank you very much!