Generalize DataOut::write_vtk for simplex mesh

[peterrum]

In the case of triangles and tetrahedrons, the input argument n_subdivisions in DataOut::build_patches() gets a new meaning: the cells are not subdivided. Instead, we switch between linear (5/22) and quadratic (10/24) VTK cell types. See also here (p. 9).

The proposed change enables:

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

#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_isoparametric.h>

#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>

#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>

#include <deal.II/simplex/fe_lib.h>
#include <deal.II/simplex/grid_generator.h>

using namespace dealii;

template <int dim>
class RightHandSideFunction : public Function<dim>
{
public:
  RightHandSideFunction(const unsigned int n_components)
    : Function<dim>(n_components)
  {}

  virtual double
  value(const Point<dim> &p, const unsigned int component = 0) const
  {
    return p[component % dim] * p[component % dim];
  }
};

template <int dim, int spacedim = dim>
void
test(const FiniteElement<dim, spacedim> &fe, const unsigned int n_components)
{
  Triangulation<dim, spacedim> tria;
  Simplex::GridGenerator::subdivided_hyper_cube(tria, 2);

  DoFHandler<dim> dof_handler(tria);

  dof_handler.distribute_dofs(fe);

  Vector<double> solution(dof_handler.n_dofs());

  MappingIsoparametric<dim> mapping(Simplex::FE_P<dim>(1));

  VectorTools::interpolate(mapping,
                           dof_handler,
                           RightHandSideFunction<dim>(n_components),
                           solution);

  DataOut<dim> data_out;

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


  data_out.build_patches(mapping, /*1 or*/ 2);

  std::ofstream output("test.vtk");
  data_out.write_vtk(output);
}

int
main()
{
  {
    const unsigned int dim = 2;
    test<dim>(Simplex::FE_P<dim>(2) /*=degree*/, 1);
    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/), dim), dim);
    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/),
                            dim,
                            Simplex::FE_P<dim>(1 /*=degree*/),
                            1),
              dim + 1);
  }
  {
    const unsigned int dim = 3;
    test<dim>(Simplex::FE_P<dim>(2) /*=degree*/, 1);
    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/), dim), dim);
    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/),
                            dim,
                            Simplex::FE_P<dim>(1 /*=degree*/),
                            1),
              dim + 1);
  }
}

The code will not compile at the moment because it depends on features from the other PRs.

FYI @DanielPaukner

[agrayver]

the input argument n_subdivisions in DataOut::build_patches() gets a new meaning: the cells are not subdivided

Depending on VtkFlags, the quads/hexes may not be subdivided, but are written as VTK high-order Lagrange cells (#6994 ).

I wonder if this should also be done for simplices as well? Then you have support for arbitrary order cells/elements and the logic of the class will remain in line with hex/quad.

I can help with the VTK part of the implementation if needed.

[peterrum]

I can help with the VTK part of the implementation if needed.

That would be great and very much appreciated! My suggestion is that we get a first working simplex implementation in deal.II (with my current implementation and all the other PRs) and then based on that we can make additional improvementsl

[agrayver]

No problem, I just wanted to bring this up and let you know.

[peterrum]

@agrayver Thanks! I will!

[peterrum]

I have rebased this branch. I will add a test once #10659 is merged. I need the mapping for creating the patches.

[kronbichler]

This looks good overall, but I have a few comments (and wait for the promised tests).

[kronbichler]

Can you explain what 'further' means?

[kronbichler]

Suggested change

	// special treatment of simplices since they are not subdivided, s.t,
	// special treatment of simplices since they are not subdivided, such that

[kronbichler]

Will that change once we support higher order simplex elements, i.e., is the if(...) continue code a temporary solution that waits for further code?

[peterrum]

Personally, I would like to keep on having the specialization since it also works (for mixed meshes) and older versions of Paraview.

[kronbichler]

You mean the specialization for the quadratic version of those or the plain linear one? It is clear that we cannot use the tensor product variant below for simplices. But then my comment simply boils down to the question whether we can restructure the code and properly indent by if / else rather than if ... continue.

[peterrum]

Done!

[peterrum]

@kronbichler I have added the test and applied the suggestions!

[peterrum]

/rebuild