Merging simplex triangulations - GridTools::consistently_order_cells

[Gordan Segon]

Hi,

Merging triangulations using GridGenerator::merge_triangulations has a call to  GridTools::consistently_order_cells(cells) before constructing the triangulation. The implemented case is only for all cells being hypercubes. Having a quick look at the implementation of `is_consistent()`, it is using GeometryInfo which is specific to hyper_cubes like line_to_cell_vertices so the same call would not work with triangles.

Now, I use merge_triangulations for both quads and triangles, and both work in my test cases. However, since in the case of a call to the function with all-triangles Triangulation<>, the call to GridTools::consistently_order_cells(cells) is simply skipped, is this something that's missing that could cause issues down the line or perhaps it is not needed in the first place?

Best,

~Gordan

[Gordan Segon]

To clarify, the `test cases` I mention are not unit tests (or otherwise) testing the merging functionality but rather just application examples where a valid triangulation is broken down in 2 parts, 1 part is transformed and a new triangulation is made by merging the 2 parts. 

~G

[Gordan Segon]

Addendum: When running GridGenerator::merge_triangulations, with quads, with copy_boundary_ids=true and if the cells are inconsistently ordered (i.e. GridTools::unnamed::reorient(cells) was called), an exception is thrown in my case with:

--------------------------------------------------------
An error occurred in line <6208> of file </home/gordan/foss-dev/dealii/cmake-build-debug/install/include/deal.II/grid/tria_accessor.h> in function
    void dealii::TriaAccessor<<anonymous>, <anonymous>, <anonymous> >::set_boundary_id(dealii::types::boundary_id) const [with int structdim = 1; int dim = 2; int spacedim = 2; dealii::types::boundary_id = unsigned int]
The violated condition was:
    boundary_ind != numbers::internal_face_boundary_id
Additional information:
    You are trying to set the boundary_id to an invalid value
    (numbers::internal_face_boundary_id is reserved).

I suspect this is due to face ordering being changed during reorientation as the code to copy the boundary ids relies on the assumption that the face ordering (indexing) is identical between resultant merged triangulation and the input triangulations. However, when `GridTools::consistently_order_cells(cells);` is commented out altogether, the merging passes without errors and the solution on the merged triangulation looks correct, or rather it doesn't look incorrect. I'll make tests to confirm that the solution is in fact correct.

To rephrase my question: Is consistent cell ordering necessary to guarantee correct downstream computations? Would the solution differ between an inconsistently and consistently ordered triangulation for the same problem in general? My guess is that the consistent ordering is utilized for performance optimizations and does not necessarily have an impact on the solution.

Best,
~G