DOF ordering of P1 functions on a structured mesh

Asked by Christian Waluga on 2013-02-21

Dear all,

I have a small question concerning the DOF-numbering. I thought that in a structured mesh, the nodes are numbered first along the x-axis, then along the y-axis and so on. A simple test yields that this is indeed the case (for the mesh):

mesh = UnitSquareMesh(3,3)
print mesh.coordinates()[:]

the numbering is like:
12 13 14 15
 8 9 10 11
 4 5 6 7
 0 1 2 3

However, when I create a linear Lagrangian space, the numbering seems to have changed. Here the DOFs seem to start in the lower-right corner and (if I recall correctly) the numbering continues somehow like:
15 13 10 6
14 11 7 3
12 8 4 1
 9 5 2 0

this has the effect that the examples dealing with the visualization of structured mesh data from the tutorial fail to work properly:

http://fenicsproject.org/documentation/tutorial/fundamentals.html#visualization-of-structured-mesh-data

A minimal example:

from dolfin import *

n = 100
mesh = RectangleMesh(-1.0, -1.0, 1.0, 1.0, n, n)
V = FunctionSpace(mesh, 'Lagrange', 1)
u0expr = Expression('((x[1]-0.25)*(x[1]-0.25)+(x[0]-0.25)*(x[0]-0.25)) < 0.025')
u = interpolate(u0expr, V)
plot(u, interactive=True)

import scitools.BoxField
import scitools.easyviz as ev
u_box = scitools.BoxField.dolfin_function2BoxField(u, mesh, (n, n), uniform_mesh=True)
ev.contour(u_box.grid.coorv[0], u_box.grid.coorv[1], u_box.values)

(my version of scitools is 0.9.0)

Does anybody happen to know how to fix this? I mean, I can easily come up with workarounds (manual renumbering, etc.), but it would be nice to stick with the builtin functionality. Perhaps I am missing something here?

Best, Christian

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
2013-02-21
Last query:
2013-02-21
Last reply:
2013-02-21
Christian Waluga (chris-1) said : #1

I just figured that this is a known issue:
https://code.google.com/p/scitools/issues/detail?id=24

Best Johan Hake (johan-hake) said : #2

Since dolfin 1.1.0 the dofs are not following the mesh vertices. To map
dof data to vertex based data you need to use a remapping:

  map = V.dofmap().vertex_to_dof_map()
  data = u.vector().array()
  data[map] = u.vector().array()

The values in data will now follow the vertex ordering.

You need to convince Hans Petter to update scitools so the built in
functionalities works as expected. He is _very_ well aware of the issue.

Johan

On 02/21/2013 06:56 PM, Christian Waluga wrote:
> New question #222509 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/222509
>
> Dear all,
>
> I have a small question concerning the DOF-numbering. I thought that in a structured mesh, the nodes are numbered first along the x-axis, then along the y-axis and so on. A simple test yields that this is indeed the case (for the mesh):
>
> mesh = UnitSquareMesh(3,3)
> print mesh.coordinates()[:]
>
> the numbering is like:
> 12 13 14 15
> 8 9 10 11
> 4 5 6 7
> 0 1 2 3
>
> However, when I create a linear Lagrangian space, the numbering seems to have changed. Here the DOFs seem to start in the lower-right corner and (if I recall correctly) the numbering continues somehow like:
> 15 13 10 6
> 14 11 7 3
> 12 8 4 1
> 9 5 2 0
>
> this has the effect that the examples dealing with the visualization of structured mesh data from the tutorial fail to work properly:
>
> http://fenicsproject.org/documentation/tutorial/fundamentals.html#visualization-of-structured-mesh-data
>
> A minimal example:
>
> from dolfin import *
>
> n = 100
> mesh = RectangleMesh(-1.0, -1.0, 1.0, 1.0, n, n)
> V = FunctionSpace(mesh, 'Lagrange', 1)
> u0expr = Expression('((x[1]-0.25)*(x[1]-0.25)+(x[0]-0.25)*(x[0]-0.25)) < 0.025')
> u = interpolate(u0expr, V)
> plot(u, interactive=True)
>
> import scitools.BoxField
> import scitools.easyviz as ev
> u_box = scitools.BoxField.dolfin_function2BoxField(u, mesh, (n, n), uniform_mesh=True)
> ev.contour(u_box.grid.coorv[0], u_box.grid.coorv[1], u_box.values)
>
> (my version of scitools is 0.9.0)
>
> Does anybody happen to know how to fix this? I mean, I can easily come up with workarounds (manual renumbering, etc.), but it would be nice to stick with the builtin functionality. Perhaps I am missing something here?
>
> Best, Christian
>
>

Christian Waluga (chris-1) said : #3

Thanks Johan. That explains the problem. Now I have to figure out how to get the unstable Version to run on OS X 10.8. I am still using the 1.1.0 binary released in late January where vertex_to_dof_map seems not to be implemented yet.

Christian Waluga (chris-1) said : #4

Thanks Johan Hake, that solved my question.

Johannes Ring (johannr) said : #5

This is fixed now in the development version of SciTools.