Dealii: subdivided_hyper_rectangle create grids with wrong material_id

[Sebastian-47]

The function subdivided_hyper_rectangle creates grids with wrong material_id.

In 2D subdivided_hyper_rectangle creates a grid with material_id = 3 and in 3D subdivided_hyper_rectangle creates a grid with material_id = 7.

Here a minimal failing example:

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>

#include <fstream>
#include <iostream>

using namespace dealii;

template <int dim>
void
generate_grid(Triangulation<dim> &triangulation)
{
  std::vector<unsigned int> n_subdivisions;
  if (dim == 2)
    {
      n_subdivisions = {2, 2};
    }
  else if (dim == 3)
    {
      n_subdivisions = {2, 2, 2};
    }
  else
    {
      Assert(false, ExcNotImplemented());
    }

  const Point<dim> left_edge =
    dim == 2 ? Point<dim>(0.0, 0.0) : Point<dim>(0.0, 0.0, 0.0);
  const Point<dim> right_edge =
    dim == 2 ? Point<dim>(1.0, 1.0) : Point<dim>(1.0, 1.0, 1.0);

  GridGenerator::subdivided_hyper_rectangle(
    triangulation, n_subdivisions, left_edge, right_edge, true);
}

int
main()
{
  { // dim = 2;
    Triangulation<2> triangulation;
    generate_grid<2>(triangulation);

    // print the grid:
    std::ofstream output_file("Grid-2d.vtu");
    GridOut().write_vtu(triangulation, output_file);
  }

  { // dim = 3;
    Triangulation<3> triangulation;
    generate_grid<3>(triangulation);

    // print the grid:
    std::ofstream output_file("Grid-3d.vtu");
    GridOut().write_vtu(triangulation, output_file);
  }

  return 0;
}

[jpthiele]

I also ran this code and can confirm this problem. However in grids.cc it seems the id is set to 0, so no idea where this error comes from.

[peterrum]

Let me check later the implementation!

[peterrum]

The reason are these lines (during coloring):

https://github.com/dealii/dealii/blob/b42193f72751984d28e11b249e47e3aebf6685f2/source/grid/grid_generator.cc#L1310-L1320

Don't ask me what the intention here is ;)

[masterleinad]

This behavior is documented at https://www.dealii.org/current/doxygen/deal.II/namespaceGridGenerator.html#ac76417d7404b75cf53c732f456e6e971

[...]
Additionally, material ids are assigned to the cells according to the octant their center is in: being in the right half plane for any coordinate direction xi adds 2i (see the glossary entry on colorization). For instance, the center point (1,-1,1) yields a material id 5 (this means that in 2d only material ids 0,1,2,3 are assigned independent from the number of repetitions).

@Sebastian-47 Why do you think this is a bug?

[Sebastian-47]

I just missed that part in the documentation.

I used the function hyper_rectangle previously and the function hyper_rectangle creates a grid with material_id = 0 (also with colorize = true). Therefore I just did not expected that the material_id will change, which lead me to thinking that this behavior was caused by a bug.

Still the behavior of subdivided_hyper_rectangle seems to be for a very special use case only.
At the moment, as a sort of workaround, I have added an extra loop to set the material_id back to zero.

[masterleinad]

Let's close this then. The behavior is documented and we can't change it without potentially breaking other code.