DataOut::write_vtu_in_parallel+ DataOutBase::VtkFlags::write_higher_order_cells + dim==3: distorted cells in Paraview

[mschreter]

While using DataOut::write_vtu_in_parallel (with #processes>1) together with DataOutBase::VtkFlags::write_higher_order_cells=true, n_subdivisons>1 and dim==3, I obtain distorted cells in Paraview (v5.10.0), see right picture in the following figure. To me it looks like a numbering issue, but I am not sure. Does someone have a suspicion what might go wrong here or have similar problems? It would be really great to make use of this flag, since otherwise the file size of the output is almost doubled.

I've used the following code to produce the results:

#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>

#include <deal.II/base/function_signed_distance.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>

#include <iostream>

using namespace dealii;

using VectorType = LinearAlgebra::distributed::Vector<double>;

template <int dim>
void test(const unsigned int fe_degree,
          const bool write_higher_order_cells = true) {
  parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);

  GridGenerator::hyper_cube(triangulation, 0, 1);
  triangulation.refine_global(2);

  MappingQ1<dim> mapping;
  FE_Q<dim> fe(fe_degree);

  DoFHandler<dim> dof_handler(triangulation);
  dof_handler.distribute_dofs(fe);

  VectorType solution(dof_handler.n_dofs());

  IndexSet locally_relevant_dofs;
  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);

  solution.reinit(dof_handler.locally_owned_dofs(), locally_relevant_dofs,
                  MPI_COMM_WORLD);

  dealii::VectorTools::interpolate(
      mapping, dof_handler, Functions::SignedDistance::Sphere<dim>(), solution);

    DataOutBase::VtkFlags flags;
    flags.write_higher_order_cells = write_higher_order_cells;

    DataOut<dim> data_out;
    data_out.set_flags(flags);

    data_out.add_data_vector(dof_handler, solution, "solution");

    data_out.build_patches(mapping, fe_degree);

    // write vtu in parallel
    std::string output("test_vtu_parallel_dim" + std::to_string(dim) + "_" +
                       std::to_string(fe_degree) + "_" +
                       std::to_string(write_higher_order_cells) + ".vtu");
    data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
}

int main(int argc, char *argv[]) {
  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

  test<3>(2);
  test<3>(2, false /*write_higher_order_cells*/);

  return 0;
}

FYI @peterrum @nmuch