Deformed cell for 3D mesh deformation in release mode

[anne-glerum]

@Kaili270 and I ran into a curious error that only occurs in Release mode, only in 3D, and only for certain combinations of X mesh extent and X mesh repetitions.

dealii version: 9.5.2 ASPECT commit (May 29 2024): 249618d20bb5c8601226c4bf48632883d4029558

Error message:

-----------------------------------------------------------------------------
--     . version 2.6.0-pre (lithosphere_with_rift_IC_cookbook, bd49ff69a)
--     . using deal.II 9.5.2
--     .       with 32 bit indices and vectorization level 3 (512 bits)
--     . using Trilinos 13.2.0
--     . using p4est 2.3.2
--     . using Geodynamic World Builder 0.6.0
--     . running in OPTIMIZED mode
--     . running with 1 MPI process
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-- For information on how to cite ASPECT, see:
--   https://aspect.geodynamics.org/citing.html?ver=2.6.0-pre&sha=bd49ff69a&src=code
-----------------------------------------------------------------------------
Number of active cells: 20 (on 1 levels)
Number of degrees of freedom: 1,248 (891+60+297)

Number of mesh deformation degrees of freedom: 180
---------------------------------------------------------
TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped
to avoid a possible deadlock.
---------------------------------------------------------


----------------------------------------------------
Exception 'ExcMessage( "The Newton iteration to find the reference point " "did not converge in 10 iterations. Do you have a " "deformed cell? (See the glossary for a definition " "of what a deformed cell is. You may want to output " "the vertices of your cell.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <977> of file </home/projects/bbp00039/dealii/dealii_v9.5.2/candi_built/tmp/unpack/deal.II-v9.5.2/source/grid/manifold.cc> in function
    dealii::Tensor<1, spacedim> dealii::FlatManifold<dim, spacedim>::normal_vector(const typename dealii::Triangulation<dim, spacedim>::face_iterator&, const dealii::Point<spacedim>&) const [with int dim = 3; int spacedim = 3; typename dealii::Triangulation<dim, spacedim>::face_iterator = dealii::TriaIterator<dealii::TriaAccessor<2, 3, 3> >]
The violated condition was: 
    iteration < 10
Additional information: 
    The Newton iteration to find the reference point did not converge in
    10 iterations. Do you have a deformed cell? (See the glossary for a
    definition of what a deformed cell is. You may want to output the
    vertices of your cell.

Stacktrace:
-----------
#0  /home/projects/bbp00039/dealii/dealii_v9.5.2/candi_built/deal.II-v9.5.2/lib/libdeal_II.so.9.5.2: dealii::FlatManifold<3, 3>::normal_vector(dealii::TriaIterator<dealii::TriaAccessor<2, 3, 3> > const&, dealii::Point<3, double> const&) const
#1  /home/projects/bbp00039/dealii/dealii_v9.5.2/candi_built/deal.II-v9.5.2/lib/libdeal_II.so.9.5.2: void dealii::VectorTools::internal::map_dofs_to_normal_vectors_and_normal_fluxes<3, 3>(dealii::DoFHandler<3, 3>::cell_iterator const&, unsigned int, std::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&, std::map<unsigned int, dealii::Function<3, double> const*, std::less<unsigned int>, std::allocator<std::pair<unsigned int const, dealii::Function<3, double> const*> > > const&, dealii::hp::FEFaceValues<3, 3>&, unsigned int, dealii::IndexSet const&, unsigned int, std::multimap<dealii::VectorTools::internal::VectorDoFTuple<3>, std::pair<dealii::Tensor<1, 3, double>, dealii::DoFHandler<3, 3>::cell_iterator>, std::less<dealii::VectorTools::internal::VectorDoFTuple<3> >, std::allocator<std::pair<dealii::VectorTools::internal::VectorDoFTuple<3> const, std::pair<dealii::Tensor<1, 3, double>, dealii::DoFHandler<3, 3>::cell_iterator> > > >&, std::map<dealii::VectorTools::internal::VectorDoFTuple<3>, dealii::Vector<double>, std::less<dealii::VectorTools::internal::VectorDoFTuple<3> >, std::allocator<std::pair<dealii::VectorTools::internal::VectorDoFTuple<3> const, dealii::Vector<double> > > >&)
#2  /home/projects/bbp00039/dealii/dealii_v9.5.2/candi_built/deal.II-v9.5.2/lib/libdeal_II.so.9.5.2: void dealii::VectorTools::internal::compute_nonzero_normal_flux_constraints_active_or_level<3, 3>(dealii::DoFHandler<3, 3> const&, unsigned int, std::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&, std::map<unsigned int, dealii::Function<3, double> const*, std::less<unsigned int>, std::allocator<std::pair<unsigned int const, dealii::Function<3, double> const*> > > const&, dealii::AffineConstraints<double>&, dealii::Mapping<3, 3> const&, dealii::IndexSet const&, unsigned int)
#3  /home/projects/bbp00039/dealii/dealii_v9.5.2/candi_built/deal.II-v9.5.2/lib/libdeal_II.so.9.5.2: void dealii::VectorTools::compute_no_normal_flux_constraints<3, 3>(dealii::DoFHandler<3, 3> const&, unsigned int, std::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&, dealii::AffineConstraints<double>&, dealii::Mapping<3, 3> const&)
#4  /home/bbpanneg/software/aspect/build__initial_conditions_continental_rift/aspect-release: aspect::MeshDeformation::MeshDeformationHandler<3>::make_initial_constraints()
#5  /home/bbpanneg/software/aspect/build__initial_conditions_continental_rift/aspect-release: aspect::MeshDeformation::MeshDeformationHandler<3>::setup_dofs()
#6  /home/bbpanneg/software/aspect/build__initial_conditions_continental_rift/aspect-release: aspect::Simulator<3>::setup_dofs()
#7  /home/bbpanneg/software/aspect/build__initial_conditions_continental_rift/aspect-release: aspect::Simulator<3>::run()
#8  /home/bbpanneg/software/aspect/build__initial_conditions_continental_rift/aspect-release: void run_simulator<3>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool)
#9  /home/bbpanneg/software/aspect/build__initial_conditions_continental_rift/aspect-release: main
--------------------------------------------------------

Aborting!

The error can be reproduced by taking the cookbook crustal_deformation/crustal_model_3D.prm and changing the Geometry model settings from

    set X extent  = 128e3
    set Y extent  =  96e3
    set Z extent  =  16e3
    set X repetitions  = 8 
    set Y repetitions  = 3

to

    set X extent = 125e3
    set Y extent = 100e3
    set Z extent = 25e3
    set X repetitions = 5 
    set Y repetitions = 4

To get rid of the error, one only has to change the X extent to 128e3. It seems like a very specific edge case where the normal vector to the surface cannot be computed. For example,

    set X extent  = 128e3
    set Y extent  =  96e3
    set Z extent  =  16e3
    set X repetitions  = 8 
    set Y repetitions  = 6

does work, even though it also has repetitions that lead to square surface elements.

The error occurs for any mesh deformation plugin (as it occurs in make_initial_constraints in the mesh_deformation interface) and we also tried several different setups. For a more complex one, the same error occurs but during make_constraints instead of make_initial_constraints in the mesh deformation interface.