- edited description
Prevent 4-vector being interpreted as 2D tensor in XDMF
Issue #1033
resolved
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)
-
reporter -
reporter - changed title to Prevent 4-vector being interpreted as 2D tensor in XDMF
- edited description
-
- changed status to resolved
Fix issue
#1033→ <<cset 8117508fb77f>>
- Log in to comment