1. FEniCS Project Avatar FEniCS Project
  2. Core
  3. DOLFIN
  4. Issues

Global_index() of vertices wrong if more than one processor and shared_facet ghost cells used

Issue #493 invalid
Stephan Schmidt created an issue

Dear all, I am experiencing bugged behavior if I access the global index of vertices with shared_facet ghost cells on more than one processor.

Please see attached minimal example: I am counting how many facets are connected to each vertex by first iterating over the facets and then incrementing a counter based on the global_index() of each vertex associated with the facet I am currently processing.

I attach two pictures of the results I get: Of course the parallel partitioning is clearly visible, which is as expected. But if you study the results with ghost cells closely, one can spot a vertex that is not attached to any facet! It is located near but not on the partition boundary in the top of the picture.

To me, either the connectivity facets -> vertices is bugged or the global_index() of said vertices. I have tested with facets.normals() and those normals seem fine, so I think the global_index() routine is bugged if ghost layers + parallel happens at the same time.

The results were created with mpiexec -np 2 python IndexBug.py

Comments (37)

  3. Chris Richardson

    The following code seems to be OK for me: 

    parameters["ghost_mode"] = "shared_facet"

mesh = UnitSquareMesh(100,100) 
mesh.init(0, 1)

fmax = 0
fmin = 100
for v in vertices(mesh, "regular"):
    nfacets = v.num_entities(1)
    fmax = max(fmax, nfacets)
    fmin = min(fmin, nfacets)

print fmax, fmin

    or how about:

    parameters["ghost_mode"] = "shared_facet"
mesh = UnitSquareMesh(100,100)

dim = mesh.topology().dim()
Touched = numpy.zeros(mesh.topology().size_global(0))

for MyFace in facets(mesh, "regular"):
    for Vert in vertices(MyFace):
        VertIndex = Vert.global_index()
        Touched[VertIndex] += 1.0

minval = 10
for i in range(len(Touched)):
    sum_proc =  MPI.sum(mesh.mpi_comm(), Touched[i])
    minval = min(minval, sum_proc)

print minval

    I agree that your code produces a strange zero, but I can't quite understand what you are doing - I'm also not sure that dof_to_vertex etc. works well in parallel...

  4. Stephan Schmidt reporter

    Sorry for posting this twice, my internet connection is fairly bad right now.

    My code is a very much stripped down minimal example of what I actually need: A FEniCS function that contains the vertex normal on boundary vertices and a vector zero inside the domain, in parallel with ghost layers. Hence, I need to count boundary vertices for averaging the facet normals to boundary vertices.

    I can confirm that the FEniCS mesh.move(V) functionality also does not work for vertices for which my minimal example above gives zero answer when used with my vertex normal CG 1 VectorFunction. It appears that there are vertices, which are not reachable when iterating over facets...

    A vertex normal (rather than a facet normal) would be a great feature for doing shape optimization and deforming geometries with FEniCS: http://www.mathematik.uni-wuerzburg.de/~schmidt/femorph

  5. Stephan Schmidt reporter

    The dof_to_vertex() has never failed me in parallel so far, with or without ghost layers enabled, it always maps from local to local indices.

  6. Jan Blechta

    dof_to_vertex_map is well-defined in parallel although I have no idea whether and how it interacts with ghost_mode. Only pitfall is that dof_to_vertex_map did not map unowned DOFs at the time the function has been written (because of global DOF numbering at that time) and I'm not sure whether the function has been rewritten to contain also unowned portion of DOFs with the switch to local DOF numbering by Garth.

  7. Jan Blechta

    I'd suggest that @Epoxid tries to prove that global indices are incorrect by a more obvious way. By definition: global index is the same on all processes for each entity, rank-wise unique and < global number of entities. Is some of those requirements violated and easily observable?

  8. Stephan Schmidt reporter

    Hmm... I will try to help where I can. But can we also include the dof_to_vertex map in this?

  9. Jan Blechta

    Well, it'd would be good if you were able to prove the claim without touching any dofmap at all. But if there's no other way...

  10. Stephan Schmidt reporter

    Oh, I am not strictly claiming that it's the global index, I am just saying, that under the circumstances described above, the chain facet -> vertex -> dof -> function seems to fail occasionally...

  11. Stephan Schmidt reporter

    Returning to the two alternatives Chris Richardson has kindly provided, without checking them in all detail, they don't produce the same answer for me when more than two processors are used. The first example at some point starts to report a minimum of fmin = 3, but the second version always reports 2 for minval, independent of the number of processors...

  12. Jan Blechta

    What about

    for MyFace in facets(mesh, "regular"):
    if True:
    #if MyFace.exterior() == True:
        for Vert in vertices(MyFace):
            VertIndex = Vert.global_index()
            Touched[VertIndex] += 1.0
for MyFace in facets(mesh, "ghost"):
    if True:
    #if MyFace.exterior() == True:
        for Vert in vertices(MyFace):
            VertIndex = Vert.global_index()
            Touched[VertIndex] += 1.0

    This solves the problem. Although I'm not sure that I understand how a ghosted facet is defined. Also notice vertex iterator in inner loop. This is probably also exclusively non-ghosted.

  13. Chris Richardson

    @Epoxid - the results seem correct to me. In the first case, as you increase the process count, some processes will no longer have the corners (with two facets per vertex). In the second case, all vertices are being summed by MPI, so it is always the same.

    @blechta - using "all" should be the same as "ghost" + "regular". A ghost facet is one which only belongs to a ghost cell (and is not shared with a regular cell).

  14. Jan Blechta

    Ok, let's make it more simple and avoid dofmaps at all:

    from dolfin import *

parameters["ghost_mode"] = "shared_facet"
mesh = UnitSquareMesh(100, 100)

# This loops all the vertices
local_indices = set()
global_indices = set()
for f in facets(mesh, "all"):
    for v in vertices(f):
        local_indices.add(v.index())
        global_indices.add(v.global_index())
assert len(local_indices) == len(global_indices)
assert len(local_indices) == mesh.num_vertices()

# This loops a subset
local_indices = set()
global_indices = set()
for f in facets(mesh, "regular"):
    for v in vertices(f):
        local_indices.add(v.index())
        global_indices.add(v.global_index())
assert len(local_indices) == len(global_indices)
assert len(local_indices) < mesh.num_vertices() or MPI.size(mesh.mpi_comm()) == 1
assert len(local_indices) == mesh.topology().ghost_offset(0)

    @chris_richardson, can you confirm the second loop (vertices of regular facets) does what intended? Should it really keep some vertices unvisited?

    BTW, there should be a check for opt=="all" in MeshEntityIterator(const Mesh& mesh, std::size_t dim, std::string opt) otherwise it keeps space for uncatched typos, like for c in cells(mesh, 'regullar'):.

  16. Jan Blechta

    Then I suggest using new-style enums, see example. For instance 

    enum class EntityType : int32_t { all = 0, regular = 1, ghost = 2 };

    which is wrapped by SWIG automatically as EntityType_all, EntityType_regular, EntityType_ghost. (I'm not sure whether this is what IPython users are dying for. Can somebody test that void list_timings(TimingClear clear, std::set<TimingType> type); works with all the IPython magic?)

  17. Chris Richardson

    Yes, I guess it can be an enum.

    The following code illustrates that all vertices are visited, even when using "regular" facets - at least, it works for me. The vertices of "regular" facets provide full coverage of all vertices of the global mesh. Maybe there is some confusion about num_vertices() which I think also includes the ghosts. 

    global_indices = set()
for f in facets(mesh, "regular"):
    for v in vertices(f):
        global_indices.add(v.global_index())

rmin = 100000
for i in range(mesh.size_global(0)):
    r = MPI.sum(mesh.mpi_comm(), int(i in global_indices))
    rmin = min(r, rmin)

print rmin
  18. Jan Blechta

    You're right. Adding line assert len(local_indices) == mesh.topology().ghost_offset(0) makes it clear.

    The question is whether the vertex coverage by (regular facets vertices) is the same as ownership of respective vertex nodes in dofmap. This would explain Stephan's problem. (I'm currently to busy to investigate this.)

  19. Stephan Schmidt reporter

    Dear all, thank you very much for the support! On a side node from a user: Might I suggest to make this difference between mesh.num_vertices() and mesh.size_global(0) more clear in the documentation?

    http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/cpp/mesh/Mesh.html#dolfin.cpp.mesh.Mesh.num_vertices and http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/cpp/mesh/Mesh.html#dolfin.cpp.mesh.Mesh.size_global

    To me as a user, reading both documentations, I would expect those to do the same, which is returning what "mesh.topology().ghost_offset(0)" actually returns...

  20. Chris Richardson

    Yes, there is always room to improve documentation.

    num_vertices() is the same is size(0) which is the number of local vertices

    size_global(0) is the number of global vertices across all processes.

    topology().ghost_offset(0) is the number of local, non-ghosted vertices.

  21. Stephan Schmidt reporter

    Dear all, thank you so much for looking into this. In the meantime, I studied the behavior of the dof_to_vertex_map with ghost facets in parallel. I think the dof_to_vertex_map returns ghost vertices instead of regular vertices, which caused the original problem? Anyway... please have a look at this minimal example:

    from dolfin import *
parameters["ghost_mode"] = "shared_facet"

mesh = UnitSquareMesh(100, 100)
mesh.init()

NSpaceS = FunctionSpace(mesh, "CG", 1)
d2v = dof_to_vertex_map(NSpaceS)
N = Function(NSpaceS)
dmin, dmax = N.vector().local_range()

for dof in range(dmax-dmin):
    vert = Vertex(mesh, d2v[dof])
    index = vert.index()
    if index > mesh.topology().ghost_offset(0):
        print "I touched a ghost vertex"

    As a user, I would expect the dof_to_vertex_map to always give me regular or physical vertices and not ghost ones...

  22. Chris Richardson

    The ghost vertices will have dofs associated with them, just as shared vertices do in a non-ghosted mesh. If you don't need them, you can just ignore them, I guess.

  23. Stephan Schmidt reporter

    Ah! Great! Thanks for the answer! I can very much understand having these dofs in the DG context. But in the CG context, somehow all of these different local dofs need to be synchronised into one value.

    I suppose this is done in the N.vector().apply("") step. Maybe the problem is rooted there? unfortunately, I could not find the apply function in the online documentation...

  24. Stephan Schmidt reporter

    Thank you again for the great discussion. Please have a look at this little program. To me, this demonstrates that the local vertex associated with the local dof number zero is never touched when looping through the local mesh...

    from dolfin import *
parameters["ghost_mode"] = "shared_facet"

mesh = UnitSquareMesh(100, 100)
mesh.init()

NSpaceS = FunctionSpace(mesh, "CG", 1)
d2v = dof_to_vertex_map(NSpaceS)
VertexOfDofZero = d2v[0]

print "Vertex of dof zero", VertexOfDofZero

Touched = False

for MyFace in facets(mesh, "regular"):
    if True:
    #if MyFace.exterior() == True:
        for Vert in vertices(MyFace):
            lIndex = Vert.index()
            if lIndex == VertexOfDofZero:
                Touched = True

if Touched == False:
    print "I have not touched the vertex for dof 0"

    I have executed the program a couple of times with different processors and it always reports that the vertex for dof 0 is not being processed when looping though the mesh. It is processed as soon as I turn off ghost.

  25. Chris Richardson

    What you say makes sense. If you loop over "regular" vertices, there will be some ghost dofs which you do not touch. Change "regular" to "all" and you will reach all dofs.

  26. Stephan Schmidt reporter

    Thank you for the great support and feedback! I took heart of your suggestions and now I loop over "all" facets instead of the "regular" facets. But I am still seeing paradox behavior. Please have a look at the program below:

    1) I am now looping over "all" facets and mark their vertices if the facet is exterior by the vertex global index.

    2) A CG 1 function is created and the result is used as the vector() of the function

    3) I loop over the function dofs and check if all marked vertices are really on the exterior and all unmarked vertices are really in the interior. This is done independently by checking their coordinates.

    The remarkable thing: The check in 3) reports all is okay: All interior vertices are marked, all exterior vertices are unmarked and vice versa.

    But if I open the output file in Paraview, atleast on my system, the file has exactly one pair of wrongly marked vertices: One interior is 1, one exterior is 0. It is as if these two are mixed-up...

    from dolfin import *
import numpy
parameters["ghost_mode"] = "shared_facet"

def VertOnSurface(vert):
    if abs(vert.x(0)) < DOLFIN_EPS: return True
    if abs(vert.x(0)-1.0) < DOLFIN_EPS: return True
    if abs(vert.x(1)) < DOLFIN_EPS: return True
    if abs(vert.x(1)-1.0) < DOLFIN_EPS: return True
    return False

mesh = UnitSquareMesh(100, 100)
mesh.init()

dim = mesh.topology().dim()
Touched = numpy.zeros(mesh.topology().size_global(0))
for MyFace in facets(mesh, "all"):
    if MyFace.exterior() == True:
        for Vert in vertices(MyFace):
            VertIndex = Vert.global_index()
            Touched[VertIndex] = 1.0

NSpaceS = FunctionSpace(mesh, "CG", 1)
d2v = dof_to_vertex_map(NSpaceS)
N = Function(NSpaceS)
dmin, dmax = N.vector().local_range()

LocValues = numpy.zeros(dmax-dmin)

for dof in range(dmax-dmin):
    Vert = Vertex(mesh, d2v[dof])
    gIndex = Vert.global_index()
    value = Touched[gIndex]
    LocValues[dof] = value
N.vector().set_local(LocValues)
N.vector().apply("")

#test if all marked vertices are really on the surface:
for dof in range(dmax-dmin):
    Vert = Vertex(mesh, d2v[dof])
    if abs(N.vector()[dof] - 1.0) < DOLFIN_EPS and VertOnSurface(Vert) == False:
        print "An interior vertex was marked!"
    if abs(N.vector()[dof]) < DOLFIN_EPS and VertOnSurface(Vert) == True:
        print "An exterior vertex was not marked!"

File("./output/N.pvd", "compressed") << N
  27. Stephan Schmidt reporter

    Here is the result of the above program when I run it on two processors in parallelN.png

  30. Stephan Schmidt reporter

    I now had the opportunity to test this further with the following minimal example that works entirely without global indices or connectivity:

    from dolfin import *
import numpy
parameters["ghost_mode"] = "shared_facet"

mesh = UnitSquareMesh(100,100)

NSpaceS = FunctionSpace(mesh, "CG", 1)
d2v = dof_to_vertex_map(NSpaceS)
N = Function(NSpaceS)
dmin, dmax = N.vector().local_range()

LocValues = numpy.zeros(dmax-dmin)

for dof in range(dmax-dmin):
    Vert = Vertex(mesh, d2v[dof])
    LocValues[dof] = Vert.x(0)
N.vector().set_local(LocValues)
N.vector().apply("")

File("./output/X1.pvd", "compressed") << N
XDMFFile(mpi_comm_world(), "./output/X1.xdmf") << N

    This little demo is just supposed to turn the x-coordinate of the (x,y) points into a CG 1 function. And it always fails to me if more than one processor is used, no matter if I use xdmf or pvd, OS X or debian. If you look closely, you will notice there is always one node that is not correctly colored. On my OS X machine, the odd vertex is on the lower edge, on my debian machine, it is exactly the lower right corner...

    Debian_Jessie.pngOSX10_10.png

  31. Chris Richardson

    Looks like dof_to_vertex_map is not right with ghosts, then.

    mesh = UnitSquareMesh(3, 3)
Q = FunctionSpace(mesh, "CG", 1)

d2vi = []
for dof in dof_to_vertex_map(Q):
    d2vi.append(Vertex(mesh, dof).global_index())

vi = []
for v in vertices(mesh, "all"):
    vi.append(v.global_index())

print sorted(d2vi), sorted(vi)

    One dof is repeated, and one missed out in the dof_to_vertex_map.

  32. Jan Blechta

    I reminf, dof_to_vertex_map does not map unowned because dofs where originally globally indexed. With "ghosts" it should be similar story. Let me look into this later this week.

  34. Stephan Schmidt reporter

    I am glad we narrowed this down! Do you have any idea of a work-around? This bug is pretty much a show stopper for my projects with FEniCS...

  38. Log in to comment
Assignee
Type
bug
Priority
critical
Status
invalid
Component
–
Milestone
–
Version
dev
Votes
3
Watchers
7
Jira: the preferred issue tracker for Bitbucket. Join the team!