Periodicity issues in deal.II version 9.3

[maurice rohracker]

Dear deal.II community,

we recently decided to update to deal.II 9.3 (9.3.1). After some tests of our periodic homogenization code, we observed some strange mesh-dependent behavior in our serial and distributed code.

In order to get to the bottom of this issue, we created a small mwe that reads the mesh, assigns boundary IDs, collects the periodic face pairs, prints those, and adds the periodicity to the triangulation.

Attached you will find the code as well as two exemplary meshes (one working with version 9.0.1 and 9.3.1 for triangulation and pdt [a100_r20_t5_fine_20_32_Mesh.inp], and another which results in an assertion while calling tria.add_periodicity() [a100_r20_t5_coarse_16_24_Mesh.inp] for the pdt using dealii 9.3). The tria type can be changed in line 280 of the code and we compiled the code with dealii 9.0.1 and 9.3.1 installed via spack.

We provide additionally the output of the periodic face pairs for all cases as well (periodicFaceCenter.txt); 2 meshes * 2 tria types * 2 dealii versions = 8 outputs. The output is the same for each mesh regardless of the version or the tria type, with which we want to show that the collection of the periodic face pairs is working properly.

However, we get the following assertion with dealii 9.3.1 on the pdt using the mesh (a100_r20_t5_coarse_16_24_Mesh.inp):

```

--------------------------------------------------------
An error occurred in line <2067> of file </tmp/ltmadmin/spack-stage/spack-stage-dealii-9.3.1-aowzkeyukhrxyzsosrunew52i4dpwd72/spack-src/source/distributed/tria.cc> in function
    bool dealii::parallel::distributed::{anonymous}::enforce_mesh_balance_over_periodic_boundaries(dealii::parallel::distributed::Triangulation<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
The violated condition was:
    topological_vertex_numbering[j] == j
Additional information:
    This exception -- which is used in many places in the library --
    usually indicates that some condition which the author of the code
    thought must be satisfied at a certain point in an algorithm, is not
    fulfilled. An example would be that the first part of an algorithm
    sorts elements of an array in ascending order, and a second part of
    the algorithm later encounters an element that is not larger than the
    previous one.
   
    There is usually not very much you can do if you encounter such an
    exception since it indicates an error in deal.II, not in your own
    program. Try to come up with the smallest possible program that still
    demonstrates the error and contact the deal.II mailing lists with it
    to obtain help.

Stacktrace:
-----------
#0  /opt/dealii-9.3.1/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.4.0/dealii-9.3.1-aowzkeyukhrxyzsosrunew52i4dpwd72/lib/libdeal_II.g.so.9.3.1:
#1  /opt/dealii-9.3.1/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.4.0/dealii-9.3.1-aowzkeyukhrxyzsosrunew52i4dpwd72/lib/libdeal_II.g.so.9.3.1: dealii::parallel::distributed::Triangulation<2, 2>::prepare_coarsening_and_refinement()
#2  /opt/dealii-9.3.1/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.4.0/dealii-9.3.1-aowzkeyukhrxyzsosrunew52i4dpwd72/lib/libdeal_II.g.so.9.3.1: dealii::parallel::distributed::Triangulation<2, 2>::copy_local_forest_to_triangulation()
#3  /opt/dealii-9.3.1/opt/spack/linux-ubuntu20.04-x86_64/gcc-9.4.0/dealii-9.3.1-aowzkeyukhrxyzsosrunew52i4dpwd72/lib/libdeal_II.g.so.9.3.1: dealii::parallel::distributed::Triangulation<2, 2>::add_periodicity(std::vector<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > >, std::allocator<dealii::GridTools::PeriodicFacePair<dealii::TriaIterator<dealii::CellAccessor<2, 2> > > > > const&)
#4  ./verify-mesh: void Check_meshPeriodicity<2u, dealii::parallel::distributed::Triangulation<2, 2>, dealii::ConditionalOStream>(dealii::parallel::distributed::Triangulation<2, 2>&, double const&, dealii::ConditionalOStream&)
#5  ./verify-mesh: main
--------------------------------------------------------

```

For the attached meshes the program is called with: ./verify-mesh ../path/to/mesh/file 100 with the specimen length 100 as also described in the README.md.

Do you have any idea where this issue comes from? The meshes have been generated with Gmsh (4.4.1).

Thanks for your help in advance.

Best regards,
Maurice Rohracker

FAU Erlangen-Nürnberg

[Paras Kumar]

Dear All,

I am trying to revive interest in this query which was sent in by my colleague some time ago.

The issue persists still, even with dealii 9.4.0 and even with meshes generated with a newer Gmsh version (4.10.0); see the two files in attachment. The mesh RVEMultiInc2D_a500_r25_Vf13_H7.5_GMSH4.10_Mesh.inp does not lead to any assertions for both serial triangulation as well as pdt but the one in RVEMultiInc2D_a500_r25_t5_Vf13_H7.5_GMSH4.10_Mesh.inp leads to the same assertion as above, with pdt. The specimen length is 500 for each of the meshes.

The fact the same mesh leads to different behaviors with serial tria and pdt, seems to point towards a possible bug in dealii.

Any leads on resolving this issue would be helpful.

Thanks and best regards,

Paras Kumar