Numbering of DG0 Elements in 1.0.0 vs 1.1.0

Asked by Thomas Fraunholz on 2013-02-18

Hallo,

when I did the upgrade on 1.1.0 I had to realize that one of my code isn't working anymore. It seems that the numbering of the DG0 Elements has changed. Is this a known issue and if yes, what is the idea to do instead, since this method is proposed in a similar way in the fenics book on page 69.

I added a commented minimum running example.

Thanks for your reply in advance,
Thomas

-----------------------------------

from dolfin import *

# Create a small mesh, if we stay here everything goes fine.
mesh = UnitSquare(4,4)

# We are doing some Refinement to produce the failure.
V = FunctionSpace(mesh, "CG", 1)
while V.dim() < 150:
    dummy_markers = MeshFunction("bool", mesh, mesh.topology().dim())
    for c in cells(mesh):
        dummy_markers[c] = True
    mesh = refine(mesh, dummy_markers)
    V = FunctionSpace(mesh, "CG", 1)

# We creat a simple marker. It is a small 2x2 chessboard.
f = Expression("(x[0]-0.5)*(x[1]-0.5)")
LocalSpace = FunctionSpace(mesh, "DG", 0)
w = TestFunction(LocalSpace)
form_estimator = f*w*dx
estimator = assemble(form_estimator, mesh=mesh)

# Create a meshfunction and set it to false
markers = MeshFunction("bool", mesh, mesh.topology().dim())
for c in cells(mesh):
    markers[c] = False

# Try to mark regions with positive values of the estimator,
# this was fine with version 1.0.0, but not with version 1.1.0.
for c in cells(mesh):
    if (estimator[c.index()] > 0.):
        markers[c] = True

# Here we can see the difference.
# Normally four square regions should be marked.
# But with version 1.1.0 everything seems to get mixed up.
plot(markers)
interactive()

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Thomas Fraunholz
Solved:
2013-02-20
Last query:
2013-02-20
Last reply:
2013-02-18
Jan Blechta (blechta) said : #1

In general indexing of dofs does not correspond to indexing of any entity. For resolution see question 221862.

Thomas Fraunholz (sleater) said : #2

Well, this is what I said. But the procedure above is not unknown to the developers. (at least Marie knows what I am talking about). So since this simple functionalty is removed, normally there should be an easy fix.

Jan Blechta (blechta) said : #3

> Well, this is what I said. But the procedure above is not unknown to
> the developers. (at least Marie knows what I am talking about). So
> since this simple functionalty is removed, normally there should be
> an easy fix.
>

There's an easy fix - you have to use dofmap to access right dofs.
Example is in Johan's answer in question 221862.

Jan

Thomas Fraunholz (sleater) said : #4

OK. I should threads to the very end. I will try Johan's hint and post my solution. I think this is interesting for more users our there.

Thomas

Thomas Fraunholz (sleater) said : #5

for c in cells(mesh):
    dof_ind = LocalSpace.dofmap().cell_dofs(c.index())[0]
    if (estimator[dof_ind] > 0.):
        markers[c] = True