Simplex grid refinement: How to avoid bad mesh quality

[kronbichler]

Together with @dominiktassilostill and @nfehn we found out that a mesh created by the following code:

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>

int main()
{
  using namespace dealii;
  constexpr int dim = 3;
  Triangulation<dim> tria;
  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1, 0, 1);
  for (unsigned int i = 0; i < 4; ++i)
    {
      FE_SimplexP<dim> fe(1);
      MappingFE<dim> mapping(fe);
      QGaussSimplex<dim> quad(3);
      const Vector<double> aspect_ratios = GridTools::compute_aspect_ratio_of_cells(mapping, tria, quad);
      DataOut<dim> data_out;
      data_out.attach_triangulation(tria);
      data_out.add_data_vector(aspect_ratios, "aspect_ratio");
      data_out.build_patches();
      std::ofstream file("grid-" + std::to_string(i) + ".vtk");
      data_out.write_vtk(file);

      tria.refine_global();
    }
}

creates increasingly worse mesh qualities as the mesh gets refined. The reason is that

dealii/source/grid/tria.cc

Lines 6745 to 6751 in acb749f

	static constexpr dealii::ndarray<unsigned int, 6, 2>
	new_line_vertices_tet = {{{{6, 8}},
	{{X, X}},
	{{X, X}},
	{{X, X}},
	{{X, X}},
	{{X, X}}}};

hardcodes the vertex connection 6 <-> 8, which then leads to the (well-documented) problem of sliver elements. The solution appears to be to also check the other two possibilities 5 <-> 7 and 4 <-> 9, and then select the best combination (which we believe is picking the line with the smallest length, because this avoids very long edges in one direction). We are in the process of figuring out the right combinations in the tables for setting up the new quads, the quad-line, quad-vertex, hex-quad, hex-line, hex-vertex connectivities for these cases (there is a large number of hardcoded tables, which in itself is already a burden). If anyone has insights or better ideas, we would be happy to discuss. Otherwise, I hope that we soon find a solution.

[nfehn]

I agree that we should try to make the implementation more "self-documenting". In this context, I suggest to avoid numbers like "4,5,6,7,8,9" for new vertices in all those arrays and tables, but to introduce variable names like v_01 (v=vertex, 01 = connecting vertices 0 and 1 of the original/unrefined tetrahedron) for the new vertices created during refinement. A new line should be called l_01_12 (l=line, connecting v01 with v12) or l_01_23 (connection v01 with v23) in my opinion.

Comments should be added to the arrays/tables explaining in which cases the ordering of certain numbers is essential, and in which cases the order is irrelevant for the code to be functional.

Principally, we should try to use common code paths for both hex and tet implementations, but only where this does not impact code readability in a negative way. Currently, I have difficulties understanding the tet-related implementation. (Ideally, one should not be able to extract from the implementation that hex has been first, and tet was implemented afterwards in my opinion.)

It could make sense to do this refactoring of the code before trying to generalize the functionality. If we now introduce new arrays/tables for the alternative variants of cutting the tetrahedron, we will further increase the burden of realizing such refactoring later.

[dominiktassilostill]

Pull request #16755 implements the 2 further options of refining a tetrahedron, which leads to improved aspect ratios of the refined elements. Sadly, I found no better way to represent the relations between the vertices, lines, faces and tetrahedrons than by hard-coding the aforementioned tables. At least the length of the last table, describing the corresponding vertices of each face of the refined tetrahedrons is shortened by the use of the reference element.

I tried adding more comments in the code to explain the meaning and the origin of the numbers in the tables to document the code, yet there is room for improvement. For clarity it may be better to refactor more code in a further pull request, in which the integration with the hex elements may be also improved by choosing common names, like 'faces' instead of 'quads'.