Prevent 4-vector being interpreted as 2D tensor in XDMF

The following MWE explains it all

import numpy as np
from fenics import *

mesh = UnitSquareMesh(10, 10)
dim = 4
V = VectorFunctionSpace(mesh, "CG", 1, dim=dim)

u = Function(V)
val = np.arange(dim)
u.interpolate(Constant(val))

# For dim == 2, will show (0, 1, 0) (padded with 0 as expected)
# For dim == 4, will show (0, 1, 0, 2) in ParaView, and not (0, 1, 2, 3)
# For other dimensions > 1, works fine
with XDMFFile("test1.xdmf") as f:
    f.write(u)

# For dim == 4, still getting (0, 1, 0, 2) in the ascii XDMFFile already
with XDMFFile("test2.xdmf") as f:
    f.write(u, XDMFFile.Encoding.ASCII)

# write/read_checkpoint works fine, but the XDMF/HDF5 files are
# unable to be visualized in ParaView for dim >= 4
with XDMFFile("test3.xdmf") as f:
    f.write_checkpoint(u, "u")
v = Function(V)
with XDMFFile("test3.xdmf") as f:
    f.read_checkpoint(v, "u")
for i in range(len(val)):
    assert np.isclose(assemble(v.sub(i) * dx), val[i])

The following diff to XDMFFile.cpp resolves the issue.

+++ b/dolfin/io/XDMFFile.cpp
@@ -2912,7 +2912,8 @@ std::vector<double> XDMFFile::get_point_data_values(const Function& u)

   std::int64_t width = get_padded_width(u);

-  if (u.value_rank() > 0)
+  const std::size_t value_rank = u.value_rank();
+  if (value_rank > 0)
   {
     // Transpose vector/tensor data arrays
     const std::size_t num_local_vertices = mesh->num_entities(0);
@@ -2922,7 +2923,7 @@ std::vector<double> XDMFFile::get_point_data_values(const Function& u)
     {
       for (std::size_t j = 0; j < value_size; j++)
       {
-        std::size_t tensor_2d_offset = (j > 1 && value_size == 4) ? 1 : 0;
+        std::size_t tensor_2d_offset = (j > 1 && value_rank == 2 && value_size == 4) ? 1 : 0;
         _data_values[i*width + j + tensor_2d_offset]
           = data_values[i + j*num_local_vertices];
       }

Comments (3)