Problem with FENothing and DoFTools::map_dofs_to_support_points()

[Davit Gyulamiryan]

Hi all,

I have encountered a problem when trying to obtain the support points from a dof_handler.

The Finite Element used is an hp::FECollection<dim> with two elements: 1) FE_System<dim> fe(FE_Q<dim>(1), dim, FE_Q<dim>(1), 1), and 2) fe_nothing(FE_Nothing<dim>(), dim, FE_Nothing<dim>(), 1). The domain is split into two regions - active and inactive - hence the use of FE_Nothing.

Then, calling DoFTools::map_dofs_to_support_points(MappingQ1<dim>(), dof_handler) gives an error in debug mode, but no errors in release mode!

The error message is as follows:

--------------------------------------------------------

Calling MPI_Abort now.
To break execution in a GDB session, execute 'break MPI_Abort' before running. You can also put the following into your ~/.gdbinit:
  set breakpoint pending on
  break MPI_Abort
  set breakpoint pending auto

--------------------------------------------------------
An error occurred in line <2791> of file </proj/petsc/dealii-9.5.2/source/fe/fe_values.cc> in function
    dealii::FEValuesBase<dim, spacedim>::FEValuesBase(unsigned int, unsigned int, dealii::UpdateFlags, const dealii::Mapping<dim, spacedim>&, const dealii::FiniteElement<dim, spacedim>&) [with int dim = 3; int spacedim = 3]
The violated condition was:
    n_q_points > 0
Additional information:
    There is nothing useful you can do with an FEValues object when using
    a quadrature formula with zero quadrature points!

Stacktrace:
-----------
#0  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: dealii::FEValuesBase<3, 3>::FEValuesBase(unsigned int, unsigned int, dealii::UpdateFlags, dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3> const&)
#1  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: dealii::FEValues<3, 3>::FEValues(dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3> const&, dealii::Quadrature<3> const&, dealii::UpdateFlags)
#2  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: dealii::FEValues<3, 3>::FEValues(dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3> const&, dealii::hp::QCollection<3> const&, dealii::UpdateFlags)
#3  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: std::_MakeUniq<dealii::FEValues<3, 3> >::__single_object std::make_unique<dealii::FEValues<3, 3>, dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3> const&, dealii::hp::QCollection<3> const&, dealii::UpdateFlags const&>(dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3> const&, dealii::hp::QCollection<3> const&, dealii::UpdateFlags const&)
#4  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: dealii::hp::FEValuesBase<3, 3, dealii::FEValues<3, 3> >::select_fe_values(unsigned int, unsigned int, unsigned int)
#5  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: void dealii::hp::FEValues<3, 3>::reinit<false>(dealii::TriaIterator<dealii::DoFCellAccessor<3, 3, false> > const&, unsigned int, unsigned int, unsigned int)
#6  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2:
#7  /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: std::map<unsigned int, dealii::Point<3, double>, std::less<unsigned int>, std::allocator<std::pair<unsigned int const, dealii::Point<3, double> > > > dealii::DoFTools::map_dofs_to_support_points<3, 3>(dealii::Mapping<3, 3> const&, dealii::DoFHandler<3, 3> const&, dealii::ComponentMask const&)
#8  test: main

--------------------------------------------------------

Which I assume has to do with the presence of FE_Nothing.

Minimum code to reproduce the issue:

int main(int argc, char *argv[])

{

  try

  {

    Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

    MPI_Comm mpi_communicator(MPI_COMM_WORLD);

    const unsigned int dim = 3;

   

    parallel::distributed::Triangulation<dim> triangulation(mpi_communicator);

    GridGenerator::hyper_cube(triangulation, 0, 1);

    triangulation.refine_global(4);

    FESystem<dim>            fe(FE_Q<dim>(1), dim, FE_Q<dim>(1), 1);

    FESystem<dim>            fe_nothing(FE_Nothing<dim>(), dim, FE_Nothing<dim>(), 1);

    hp::FECollection<dim>    fe_collection;

    fe_collection.push_back(fe);

    fe_collection.push_back(fe_nothing);

    DoFHandler<dim> dof_handler(triangulation);

    for (const auto &cell : dof_handler.active_cell_iterators())

    {

      if (cell->is_locally_owned())

      {

        if (cell->center()[0] < 0.5)

          cell->set_active_fe_index(0);

        else

          cell->set_active_fe_index(1);

      }

    }

    dof_handler.distribute_dofs(fe_collection);

    MappingQ1<dim> mapping;

    std::map<types::global_dof_index, Point<dim> > support_points;

    support_points = DoFTools::map_dofs_to_support_points(mapping, dof_handler);

  }

  catch (std::exception &exc)

  {

    std::cerr << std::endl

              << std::endl

              << "----------------------------------------------------"

              << std::endl;

    std::cerr << "Exception on processing: " << std::endl

              << exc.what() << std::endl

              << "Aborting!" << std::endl

              << "----------------------------------------------------"

              << std::endl;

    return 1;

  }

  catch (...)

  {

    std::cerr << std::endl

              << std::endl

              << "----------------------------------------------------"

              << std::endl;

    std::cerr << "Unknown exception!" << std::endl

              << "Aborting!" << std::endl

              << "----------------------------------------------------"

              << std::endl;

    return 1;

  }

  return 0;

}

I would appreciate any help! Let me know if there is any other information required about the case.

Best wishes,

Davit Gyulamiryan

[Bruno Turcksin]

Hello Davit,

This does not work because FE_Nothing does not have a support points. I needed something similar so I implemented it in my code, you can take a look at https://github.com/adamantine-sim/adamantine/blob/efde0a1818ae4f2c476aee56c55a187c7b3ece08/source/experimental_data_utils.cc#L21-L65 You may need to slightly adapt it for your use case.

Best,

Bruno

[Wolfgang Bangerth]

> Then, calling DoFTools::map_dofs_to_support_points(MappingQ1<dim>(),

> dof_handler) gives an error in *debug mode*, but *no *errors in *release mode*!

There are ~14,000 checks in the library that only run in debug mode. They
typically verify that input arguments of functions are reasonable. You ran
into one of them. But just because you don't get an error in release mode
(because the check is not run in release mode) does not mean that the result
you get is reasonable.

>
> The error message is as follows:
> --------------------------------------------------------
> Calling MPI_Abort now.
> To break execution in a GDB session, execute 'break MPI_Abort' before running.
> You can also put the following into your ~/.gdbinit:
>   set breakpoint pending on
>   break MPI_Abort
>   set breakpoint pending auto
>
> --------------------------------------------------------
> An error occurred in line <2791> of file
> </proj/petsc/dealii-9.5.2/source/fe/fe_values.cc> in function
>     dealii::FEValuesBase<dim, spacedim>::FEValuesBase(unsigned int, unsigned
> int, dealii::UpdateFlags, const dealii::Mapping<dim, spacedim>&, const
> dealii::FiniteElement<dim, spacedim>&) [with int dim = 3; int spacedim = 3]
> The violated condition was:
>     n_q_points > 0
> Additional information:
>     There is nothing useful you can do with an FEValues object when using
>     a quadrature formula with zero quadrature points!

> [...]

I'd consider this a bug. I opened a github issue here:
https://github.com/dealii/dealii/issues/18279
I don't know when someone will get around to fixing it, but at least the issue
is recorded.

> Minimum code to reproduce the issue:

Great little test, thanks for coming up with a minimum piece of code that
shows the issue!

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Wolfgang Bangerth]

On 3/23/25 10:48, Wolfgang Bangerth wrote:
>
> I'd consider this a bug. I opened a github issue here:
> https://github.com/dealii/dealii/issues/18279
> I don't know when someone will get around to fixing it, but at least the issue
> is recorded.

It turns out that this was easy enough to fix:
https://github.com/dealii/dealii/pull/18280

[Davit Gyulamiryan]

Great! Thank you for your reply.

Best,

Davit