dealii / dealii

The development repository for the deal.II finite element library
https://www.dealii.org
Other
1.32k stars 736 forks source link

FiniteElement::(restriction|prolongation)_is_implemented() is broken #16201

Open tamiko opened 10 months ago

tamiko commented 10 months ago

When doing something like this:

#include <deal.II/base/config.h>
#include <deal.II/fe/fe_hermite.h>
int main()
{
  constexpr int dim = 2;
  dealii::FE_Hermite<dim> finite_element(1); // does not implement restriction/prolongation
  std::cout << finite_element.restriction_is_implemented() << std::endl;
  return 0;
}

I expect the program to simply print a 0 for false. Instead it throws an exception:

--------------------------------------------------------
An error occurred in line <345> of file </home/tamiko/workspace/dealii/source/fe/fe.cc> in function
    const dealii::FullMatrix<double>& dealii::FiniteElement<<anonymous>, <anonymous> >::get_restriction_matrix(unsigned int, const dealii::RefinementCase<dim>&) const [with int dim = 2; int spacedim = 2]
The violated condition was: 
    restriction[refinement_case - 1][child].n() == this->n_dofs_per_cell()
Additional information: 
    You are trying to access the matrices that describe how to restrict a
    finite element function from the children of one cell to the finite
    element space defined on their parent (i.e., the 'restriction' or
    'projection' matrices). However, the current finite element can either
    not define this sort of operation, or it has not yet been implemented.

This is due to the fact that we call this function in FiniteElement<dim, spacedim>::restriction_is_implemented():

761         // make sure also the lazily initialized matrices are created            
762         get_restriction_matrix(c, RefinementCase<dim>(ref_case));

because we base the decision whether we have restriction/prolongation matrices on the fact that they are initialized... This is a problem because get_restriction_matrix(...) has to be able to signal an error condition (which it happily does with a call to Assert(...).).

tamiko commented 10 months ago

I am opening this issue to get some consensus on how to fix this - before I touch about every finite element to switch over to the Lazy<> wrapper for restriction/prolongation matrices. :smiley:

My gut feeling is that this method actually has very little utility... (and the documentation clearly states this as well).

masterleinad commented 10 months ago

I'm fine with just deprecating these functions. it's easy enough to duplicate the implementation in the tests if we really want it.

bangerth commented 10 months ago

The functions don't appear to be used widely. I can only find these places:

1/dealii> grep -r prolongation_is_implemented tests
tests/fe/fe_prolongation_common.h:      if (fe.isotropic_prolongation_is_implemented())
1/dealii> grep -r restriction_is_implemented tests
tests/fe/fe_restriction_common.h:      if (fe.isotropic_restriction_is_implemented())
tests/bits/fe_tools_08.cc:  if (!fe2.isotropic_restriction_is_implemented())
tests/mpi/fe_tools_extrapolate_common.h:  if (!fe2.isotropic_restriction_is_implemented())
tests/mpi/fe_tools_extrapolate_common.h:  if (!fe2.isotropic_restriction_is_implemented())

So yes, deprecation sounds reasonable if it's not easy to fix them.