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

Putting MeshFunctions in XDMF causes segfaults on reading

Issue #12 resolved
Patrick Farrell created an issue

Consider the following code:

from dolfin import *

mesh = UnitIntervalMesh(2)
facets = MeshFunction("size_t", mesh, 0)
xdmf = XDMFFile("test.xdmf")
xdmf << mesh
xdmf << facets # comment this line to make it work
mesh = Mesh("test.xdmf")

Reading in the mesh from the XDMF file fails if the facets are dumped:

$ python readxdmf.py 
Traceback (most recent call last):
  File "readxdmf.py", line 15, in <module>
    mesh = Mesh("test.xdmf")
  File "/data/pfarrell/fenics-trunk/lib/python2.7/site-packages/dolfin/mesh/meshes.py", line 68, in __init__
    cpp.Mesh.__cppinit__(self, *args, **kwargs)
  File "/data/pfarrell/fenics-trunk/lib/python2.7/site-packages/dolfin/cpp/mesh.py", line 1872, in __init__
    _mesh.Mesh_swiginit(self,_mesh.new_Mesh(*args))
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to complete call to function operator>>().
*** Reason:  Assertion xdmf_topology failed.
*** Where:   This error was encountered inside /data/pfarrell/dolfin/git-dolfin/dolfin/io/XDMFFile.cpp (line 311).
*** Process: 0
*** -------------------------------------------------------------------------

but works fine if the marked line is commented out.

Comments (21)

  2. Patrick Farrell reporter

    This also doesn't work, for a different reason:

    from dolfin import *

mesh = UnitIntervalMesh(2)
facets = MeshFunction("size_t", mesh, 0)

xdmf = XDMFFile("mesh.xdmf")
xdmf << mesh
mesh = Mesh("mesh.xdmf")

xdmff = XDMFFile("facets.xdmf")
xdmff << mesh
xdmff << facets
facets = MeshFunction("size_t", mesh, "facets.xdmf")

    Traceback (most recent call last):
  File "readxdmf.py", line 16, in <module>
    xdmff << facets
Exception: std::bad_alloc

    How should you use facets with XDMF/HDF5 files?

  3. Chris Richardson

    The XDMF interface was originally intended for visualisation, and it uses the GenericFile interface with >> and << which is not well suited to fine-grained control. For example, when reading a mesh from XDMF, which mesh should it read when there are multiple meshes in the HDF5 file?

    Nevertheless, it should bail out a bit more gracefully, when failing... I'll have a look at fixing it.

    You may find it more reliable to use HDF5File instead:

    hdf5 = HDF5File("example.h5","w")
hdf5.write(mesh, "my_mesh")
hdf5.write(meshfunction, "my_meshfunction")



mesh2=Mesh()
hdf5 = HDF5File("example.h5","r")
hdf5.read(mesh2, "my_mesh")

    etc.

  4. Patrick Farrell reporter

    I'm still having trouble with facets, even with the HDF5File interface:

    from dolfin import *

m = UnitIntervalMesh(2)
f = MeshFunction("size_t", m, 0)
h = HDF5File("test.h5", "w")
h.write(m, "mesh")
h.write(f, "facets")

m = Mesh()
h = HDF5File("test.h5", "r")
h.read(m, "mesh")
f = MeshFunction("size_t", m, 1) # can't do f = MeshFunction()
h.read(f, "facets") # fails

    [pef@aislinn:/tmp]$ python bugs.py 
Traceback (most recent call last):
  File "bugs.py", line 15, in <module>
    h.read(f, "facets") # fails
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to read meshfunction topology.
*** Reason:  Cell dimension mismatch.
*** Where:   This error was encountered inside HDF5File.cpp.
*** Process: 0
*** -------------------------------------------------------------------------

    I tried looking for demo code in the demo/undocumented directory, but there doesn't seem to be anything to do with HDF5 there:

    [pef@aislinn:~/src/dolfin/dolfin/demo/undocumented]$ grep -r HDF5 .
[pef@aislinn:~/src/dolfin/dolfin/demo/undocumented]$ grep -r XDMF .
[pef@aislinn:~/src/dolfin/dolfin/demo/undocumented]$

    Are there any demos anywhere?

  5. Chris Richardson

    This should work:

    from dolfin import *

mesh = UnitSquareMesh(2,2)
facets = MeshFunction("size_t", mesh, 1)

xdmf = XDMFFile("mesh.xdmf")
xdmf << mesh
xdmff = XDMFFile("facets.xdmf")
xdmff << facets

del xdmf
del xdmff

mesh = Mesh("mesh.xdmf")
facets = MeshFunction("size_t", mesh, 1)

xdmff = File("facets.xdmf")
xdmff >> facets
  6. Chris Richardson

    Good point - there should be demos. You can find some examples of usage in test/unit/io/python for now...

  8. Chris Richardson

    Also, in your hdf5 example above, there really is a cell dimension mismatch - should be "0" for both MeshFunctions... try changing it, and I think it will work.

  9. Patrick Farrell reporter

    So it does! I had it in my head that it was the default value, but it's actually not .. Thanks for pointing that out.

  10. Chris Richardson

    I think I have now fixed:

    f = MeshFunction('size_t', mesh, 'filename.xdmf')

    in my fork "chris_richardson/phdf5" - I guess I'll try and get it merged soon.

  12. Patrick Farrell reporter

    Running the original code in the bug report with the latest dolfin trunk:

    [pef@caoimhe:/tmp/issue12]$ cat report1.py 
from dolfin import *

mesh = UnitIntervalMesh(2)
facets = MeshFunction("size_t", mesh, 0)
xdmf = XDMFFile("test.xdmf")
xdmf << mesh
xdmf << facets # comment this line to make it work
mesh = Mesh("test.xdmf")
[pef@caoimhe:/tmp/issue12]$ python report1.py 
Traceback (most recent call last):
  File "report1.py", line 8, in <module>
    mesh = Mesh("test.xdmf")
  File "/data/pfarrell/fenics-trunk/lib/python2.7/site-packages/dolfin/mesh/meshes.py", line 68, in __init__
    cpp.Mesh.__cppinit__(self, *args, **kwargs)
  File "/data/pfarrell/fenics-trunk/lib/python2.7/site-packages/dolfin/cpp/mesh.py", line 1872, in __init__
    _mesh.Mesh_swiginit(self,_mesh.new_Mesh(*args))
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to complete call to function operator>>().
*** Reason:  Assertion xdmf_topology failed.
*** Where:   This error was encountered inside /data/pfarrell/dolfin/git-dolfin/dolfin/io/XDMFFile.cpp (line 311).
*** Process: 0
*** -------------------------------------------------------------------------

[pef@caoimhe:/tmp/issue12]$

    Running the code in my first comment:

    [pef@caoimhe:/tmp/issue12]$ cat report2.py 
from dolfin import *

mesh = UnitIntervalMesh(2)
facets = MeshFunction("size_t", mesh, 0)

xdmf = XDMFFile("mesh.xdmf")
xdmf << mesh
mesh = Mesh("mesh.xdmf")

xdmff = XDMFFile("facets.xdmf")
xdmff << mesh
xdmff << facets
facets = MeshFunction("size_t", mesh, "facets.xdmf")
[pef@caoimhe:/tmp/issue12]$ python report2.py 
[caoimhe:14939] *** Process received signal ***
[caoimhe:14939] Signal: Segmentation fault (11)
[caoimhe:14939] Signal code: Address not mapped (1)
[caoimhe:14939] Failing at address: 0xf0
[caoimhe:14939] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0xfcb0) [0x7f14839d4cb0]
[caoimhe:14939] [ 1] /data/pfarrell/dolfin/git-dolfin/build/dolfin/libdolfin.so.1.2(_ZNK6dolfin12MeshTopology4sizeEm+0x32) [0x7f147cc97cb2]
[caoimhe:14939] [ 2] /data/pfarrell/dolfin/git-dolfin/build/dolfin/libdolfin.so.1.2(_ZNK6dolfin8HDF5File34reorder_vertices_by_global_indicesERKNS_4MeshE+0x6f) [0x7f147cbbf81f]
[caoimhe:14939] [ 3] /data/pfarrell/dolfin/git-dolfin/build/dolfin/libdolfin.so.1.2(_ZN6dolfin8HDF5File5writeERKNS_4MeshEmSs+0x218) [0x7f147cbbfc38]
[caoimhe:14939] [ 4] /data/pfarrell/dolfin/git-dolfin/build/dolfin/libdolfin.so.1.2(_ZN6dolfin8HDF5File19write_mesh_functionImEEvRKNS_12MeshFunctionIT_EESs+0xfe) [0x7f147cbd672e]
[caoimhe:14939] [ 5] /data/pfarrell/dolfin/git-dolfin/build/dolfin/libdolfin.so.1.2(_ZN6dolfin8XDMFFile19write_mesh_functionImEEvRKNS_12MeshFunctionIT_EE+0x268) [0x7f147cc60c78]
[caoimhe:14939] [ 6] /data/pfarrell/fenics-trunk/lib/python2.7/site-packages/dolfin/cpp/_io.so(+0x4215d) [0x7f145d4a815d]
[caoimhe:14939] [ 7] python(PyObject_Call+0x36) [0x4e9f36]
[caoimhe:14939] [ 8] python() [0x4ec11a]
[caoimhe:14939] [ 9] python(PyObject_Call+0x36) [0x4e9f36]
[caoimhe:14939] [10] python() [0x478de8]
[caoimhe:14939] [11] python() [0x539cde]
[caoimhe:14939] [12] python() [0x42587d]
[caoimhe:14939] [13] python(PyEval_EvalFrameEx+0x20ac) [0x499cac]
[caoimhe:14939] [14] python(PyEval_EvalCodeEx+0x1a0) [0x49f1c0]
[caoimhe:14939] [15] python(PyRun_FileExFlags+0xe1) [0x4a9081]
[caoimhe:14939] [16] python(PyRun_SimpleFileExFlags+0x1d1) [0x4a9311]
[caoimhe:14939] [17] python(Py_Main+0x55d) [0x4aa8bd]
[caoimhe:14939] [18] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed) [0x7f14826e776d]
[caoimhe:14939] [19] python() [0x41b9b1]
[caoimhe:14939] *** End of error message ***
Segmentation fault (core dumped)

    So I don't think it should be closed yet.

  13. Chris Richardson

    Case 1 is failing at an assert statement, which might be expected, as it is trying to find a suitable Mesh to read in from an XDMF file with multiple meshes, and probably finding the facet data instead... This should be more robustly checked for.

    Case 2 is not so clear - possibly caused by trying to read a file which is still open for writing... I'll see if I can make it behave more gracefully.

    Generally, trying to read in XDMF files with multiple datasets in them is asking for trouble. Use HDF5File instead.

  14. Chris Richardson

    I've added some more error messages to pick up "report1.py" type errors.

    For your second case, you need to rewrite your code like this:

    from dolfin import *

mesh = UnitIntervalMesh(20)
facets = MeshFunction("size_t", mesh, 0)

xdmf = XDMFFile("mesh.xdmf")
xdmf << mesh
# delete XDMFFile object, so it will close the associated HDF5 file
del xdmf

xdmff = XDMFFile("facets.xdmf")
xdmff << facets
# delete XDMFFile object, so it will close the associated HDF5 file
del xdmff

# reopen XDMF files in read mode
mesh = Mesh("mesh.xdmf")
facets = MeshFunction("size_t", mesh, "facets.xdmf")

    The reason is that XDMFFile keeps a reference to its associated H5 file, and you were effectively opening the H5 file again for read, without closing it for write beforehand.

  15. Patrick Farrell reporter

    Ah, I see. Is there any way to check for that in the code, to stop people making the same mistake?

  16. Chris Richardson

    I think it would require some kind of file locking - which may be OS dependent. Maybe it should be an enhancement, but perhaps low priority...

    Your code, does however, also expose another kind of bug. Try this:

    from dolfin import *
mesh = UnitSquareMesh(10,10)
meshfunction = MeshFunction('size_t', mesh, 1, 0)

mesh = Mesh()
print meshfunction.mesh()

    [rook:13568] *** Process received signal ***
[rook:13568] Signal: Segmentation fault (11)
[rook:13568] Signal code: Invalid permissions (2)
[rook:13568] Failing at address: 0x7ff5b91c6798
[rook:13568] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0xfcb0) [0x7ff5b9af6cb0]
[rook:13568] [ 1] /lib/x86_64-linux-gnu/libc.so.6(+0x3b8798) [0x7ff5b91c6798]
[rook:13568] *** End of error message ***
Segmentation fault (core dumped)

    Should this happen?

  17. Prof Garth Wells

    MeshFunction should use a shared_ptr to store the Mesh. This would fix the above issue pointed out by Chris.

  18. Johan Hake

    That would be nice. We should probably go over the whole interface to check that we use shared ptr where it makes sense. J

  22. Log in to comment
Assignee
–
Type
bug
Priority
major
Status
resolved
Component
–
Milestone
–
Version
–
Votes
0
Watchers
2
Jira: the preferred issue tracker for Bitbucket. Join the team!