Meaning of cell_index and cell_dof_index in xml Function file

Question

Hi,

I am trying to write a Matlab script that creates a dolfin xml file containing a function. I have successfully written a script to write an xml mesh file, but I am struggling with the function file (can't find any documentation on this, sorry for the question if I have missed it). To try to understand the contents of an xml function file, I created a simple file as follows, which just writes the constant 8.5 on a 2x2 square mesh (i.e. 3x3=9 nodes).

from dolfin import *
mesh = UnitSquareMesh(2,2)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
c = Constant(8.5)
f.interpolate(c)
File('func.xml') << f

The output file func.xml is as follows:

<?xml version="1.0"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <function_data size="9">
    <dof index="0" value="8.5" cell_index="5" cell_dof_index="1" />
    <dof index="1" value="8.5" cell_index="1" cell_dof_index="1" />
    <dof index="2" value="8.5" cell_index="4" cell_dof_index="2" />
    <dof index="3" value="8.5" cell_index="0" cell_dof_index="0" />
    <dof index="4" value="8.5" cell_index="0" cell_dof_index="2" />
    <dof index="5" value="8.5" cell_index="6" cell_dof_index="2" />
    <dof index="6" value="8.5" cell_index="0" cell_dof_index="1" />
    <dof index="7" value="8.5" cell_index="2" cell_dof_index="2" />
    <dof index="8" value="8.5" cell_index="2" cell_dof_index="1" />
  </function_data>
</dolfin>

I assume 'dof index' refers to the 9 nodes, 'value' is obvious, but I don't understand the meaning of 'cell_index' or 'cell_dof_index'. If I wanted to create func.xml from outside of fenics, how would I determine these values? I tried setting these all to zero, but that seems to cause problems when reading/plotting the function.

Thanks.

Answer 1

The cell_index refers to the cell in the mesh and the cell_dof_index to the degree of freedom with respect to that particular cell. The DofMapper/Mesh then e.g. know that vertex 2 of cell 8 is vertex 42 of the global mesh which corresponds to DOF 23. The format now requires that you don't only write out absolute but relative DOF indices here. I guess this is to avoid problems when the internal ordering is changed for some reason (as it did recently).

Say you have a CG1 function living on a triangulation which is given in MATLAB as a connectivity matrix in the following format (I'm completely ignoring the Fortran-style indexing of MATLAB in the following... this should of course not be ignored later)

t = [0, 1, 2; 1, 3, 2 ...]

Then each line corresponds to a cell_index and each entry of that line to a DOF-index (assuming you order your DOFs according to your vertices). In cell 0, you would have the cell_dof_index equal to the global DOF indices, but in cell 1, local index 0 corresponds to global DOF 1, local index 1 to global DOF 3 etc.

What you should be aware of are two things:

• MATLAB indexing starts with 1, while it starts with 0 in Python.
• The internal DOF-ordering in DOLFIN may differ from the node-ordering. Here a vertex_to_dof_map might need to be used.

By the way: maybe you are reinventing the wheel. You might want to have a look at: http://github.com/FEniCS/dolfin/tree/master/utils/matlab

Comments

Thanks Christian. That makes sense. One thing that confused me was that although the mesh had 8 cells, cell_index in the xml function file only referred to 6 of them (cells 3 and 7 are not there). However, each vertex belongs to multiple cells, and since each dof index is only linked with one cell in the xml file, I can see how this is possible.

Thanks for pointing me to the Matlab utils. There are m-files for writing meshes, vectors and matrices, but as far as I can tell, there isn't one for writing a function xml similar to the one above.

Answer 2

I've come up with a much easier solution than writing an xml function file (this seems to be quite difficult as it's not obvious which cell_index should be chosen for a given vertex... there are multiple possibilities and choosing the wrong one seems to create 'holes' in the function).

The solution is just to write out the values into a simple text file in the same vertex ordering as in the mesh file. Then import into a numpy array via numpy.loadtxt and write this to the function vector. Very important is to set parameters["reorder_dofs_serial"] = False. An example is given below.

from dolfin import *
parameters["reorder_dofs_serial"] = False
mesh = Mesh('mesh.xml')
V = FunctionSpace(mesh, 'CG', 1)
f = Function(V)
a = loadtxt('ao_u.txt')
f.vector()[:] = a

Comments

This is nice solution but what if you work with some space where DOFs do not correspond to vertices e.g.

V = FunctionSpace(mesh, 'CG', 2)

?