Saving component of a mixed function to XML causes exception

[Christian Waluga]

Hello,

I just ran into another problem, which you can reproduce using the following code:

from dolfin import *

mesh = UnitSquareMesh(10,10)

V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'DG', 0)
W = V*Q

w = Function(W)
v, q = w.split()

File('v-1.xml') << interpolate(v, V)
File('v-2.xml') << v

When I use the workaround with the interpolation, everything works fine. However, the second line (v-2.xml) causes an exception.

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Python(43896) malloc: *** error for object 0x7f8f01432de0: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug
[feynman:43896] *** Process received signal ***
[feynman:43896] Signal: Abort trap: 6 (6)
[feynman:43896] Signal code:  (0)
[feynman:43896] [ 0] 2   libsystem_c.dylib                   0x00007fff8c9fd94a _sigtramp + 26
[feynman:43896] [ 1] 3   libsystem_c.dylib                   0x00007fff8ca20002 asl_file_fetch_pos + 408
[feynman:43896] [ 2] 4   libsystem_c.dylib                   0x00007fff8ca289b9 free + 392
[feynman:43896] [ 3] 5   libdolfin.1.2.dylib                 0x0000000110be43a3 _ZNSt6vectorISt4pairIiiESaIS1_EE13_M_insert_auxEN9__gnu_cxx17__normal_iteratorIPS1_S3_EERKS1_ + 405
[feynman:43896] [ 4] 6   libdolfin.1.2.dylib                 0x0000000110be23e8 _ZN6dolfin15XMLFunctionData24build_global_to_cell_dofERSt6vectorIS1_ISt4pairIiiESaIS3_EESaIS5_EERKNS_13FunctionSpaceE + 1196
[feynman:43896] [ 5] 7   libdolfin.1.2.dylib                 0x0000000110be27f6 _ZN6dolfin15XMLFunctionData5writeERKNS_8FunctionEN4pugi8xml_nodeE + 418
[feynman:43896] [ 6] 8   libdolfin.1.2.dylib                 0x0000000110bcb1e5 _ZN6dolfin7XMLFilelsERKNS_8FunctionE + 113
[feynman:43896] [ 7] 9   _io.so                              0x0000000117f33497 _wrap_File___lshift__ + 24471
[feynman:43896] [ 8] 10  Python                              0x000000010f4128b1 PyObject_Call + 97
[feynman:43896] [ 9] 11  Python                              0x000000010f41e526 instancemethod_call + 502
[feynman:43896] [10] 12  Python                              0x000000010f4128b1 PyObject_Call + 97
[feynman:43896] [11] 13  Python                              0x000000010f46cbf0 call_maybe + 288
[feynman:43896] [12] 14  Python                              0x000000010f46a234 slot_nb_lshift + 308
[feynman:43896] [13] 15  Python                              0x000000010f40f574 binary_op1 + 212
[feynman:43896] [14] 16  Python                              0x000000010f40f1d7 PyNumber_Lshift + 23
[feynman:43896] [15] 17  Python                              0x000000010f4a498a PyEval_EvalFrameEx + 13226
[feynman:43896] [16] 18  Python                              0x000000010f4a1596 PyEval_EvalCodeEx + 1990
[feynman:43896] [17] 19  Python                              0x000000010f4a0dc6 PyEval_EvalCode + 54
[feynman:43896] [18] 20  Python                              0x000000010f4c7e3e PyRun_FileExFlags + 174
[feynman:43896] [19] 21  Python                              0x000000010f4c79a0 PyRun_SimpleFileExFlags + 768
[feynman:43896] [20] 22  Python                              0x000000010f4db7a8 Py_Main + 2952
[feynman:43896] [21] 23  libdyld.dylib                       0x00007fff92ff37e1 start + 0
[feynman:43896] *** End of error message ***
Abort trap: 6

MacOS X 10.8 using the current nightly binary of FEniCS. I have observed similar issues with previous versions though.

Best, Christian

[Jan Blechta]

What is the purpose of storing subfunction to XML? I can't see any.

[Christian Waluga]

For example if you have two components, say velocity and pressure, and you only want to store one. Of course I could always project then or save the other component as well... Then again, if the code design suggests that you can do this, why does it crash? In any case I think the code should be either throw an error or show the expected behavior.

[Jan Blechta]

So you expect that you store function in space W.sub(0) to XML and later you initialize some subfunction by code

w = Function(W)
v, q = w.split()
File('v-2.xml') >> v

I doubt that this initialization is supported (I don't claim it should not be) so your output is useless. You could not initialize anything else than W.sub(0) as this is different space than V.

[Christian Waluga]

Later I want to do the following:

v = Function(V, 'v-2.xml')

As I said, I can always work around this, but in my opinion this should work.

I mean, you can also do the following without problems...

from dolfin import *

mesh = UnitSquareMesh(10,10)

V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'DG', 0)
W = V*Q

w = Function(W)
v, q = w.split()

File('v-1.pvd') << v

If this works for the VTK output, why doesn't it work for the XML output?

[Jan Blechta]

You want to store W.sub(0) function and load V function but this cannot work because W.sub(0) and V are different spaces. What you call workaround is logical practice - you create V function from W.sub(0) function, store V function and load V function.

VTK output is different. It is used for visualization and not for initializing FEM functions.

I don't say that XML I/O for subfunctions can't work. But it would need to be used consistently, not mixed with non-mixed functions like you suggest.

[Christian Waluga]

I don't see where W.sub(0) and V are -mathematically- that different. Sure, the internal DOF-mapping etc. is different, but those are IMHO the details that dolfin tries to keep away from the user. Anyway, I was just wondering why this doesn't work. If it is not a bug but a feature, I am cool with that :-). In that case you can consider this issue resolved.

[Jan Blechta]

Watch issue 84 as this could affect related business.

I can't close this issue. You could do it to reduce developers' load.

[Prof Garth Wells]

It doesn't make sense to save a shallow sub-function since it's just a view into a larger data structure. However, a sensible error message should be printed, which does not presently seem to be the case.

[Prof Garth Wells]

Fixes issue #93.

→ <<cset 1f390285a591>>