Writing DG vectors (dim >= 3) to XDMF gives node-centered fields

[Tianyi Li]

Correct behavior when writing a scalar DG field to a XDMF file

...
<Attribute Name="Test" AttributeType="Scalar" Center="Cell">
    <DataItem Dimensions="8 1" Format="HDF">u.h5:/VisualisationVector/0</DataItem>
</Attribute>
...

However writing a vector DG field returns an interpolated nodal field

...
<Attribute Name="Test" AttributeType="Vector" Center="Node">
    <DataItem Dimensions="8 3" Format="HDF">u.h5:/VisualisationVector/0</DataItem>
</Attribute>
...

MWE

from fenics import *
import numpy as np

mesh = UnitSquareMesh(2, 2)
num_cells = 8

# Write a scalar DG field to a XDMF file
# Correct behavior, Center="Cell" in XDMF
V = FunctionSpace(mesh, "DG", 0)
u = Function(V, name="Test")
values = np.arange(num_cells)
u.vector()[:] = values.flatten()

with XDMFFile("u.xdmf") as f:
    f.write(u, 0)

# Write a vector (dim >= 3) DG field to a XDMF file
# With dim = 2, correct behavior, Center="Cell" in XDMF
# With dim >= 3, wrong behavior, Center="Node" in XDMF
dim = 3
V = VectorFunctionSpace(mesh, "DG", 0, dim=dim)
v = Function(V, name="Test")
values = np.zeros((num_cells, dim))
values[:, 0] = np.arange(num_cells)
values[:, 1] = -np.arange(num_cells)

v.vector()[:] = values.flatten()

with XDMFFile("v.xdmf") as f:
    f.write(v, 0)

Suspect

bool XDMFFile::has_cell_centred_data(const Function& u)
{
  // Test for cell-centred data
  std::size_t cell_based_dim = 1;
  for (std::size_t i = 0; i < u.value_rank(); i++)
    cell_based_dim *= u.function_space()->mesh()->topology().dim();
  return (u.function_space()->dofmap()->max_element_dofs() == cell_based_dim);
}