Replace FE_DGP with custom class avoiding assert

[fdrmrc]

A new class which consists of a copy of FE_DGP in which this assert is not triggered.

template <int dim, int spacedim>
std::vector<std::pair<unsigned int, unsigned int>>
FE_DGP<dim, spacedim>::hp_vertex_dof_identities(
  const FiniteElement<dim, spacedim> &fe_other) const
{
  // there are no such constraints for DGP elements at all
  if (dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe_other) != nullptr)
    return std::vector<std::pair<unsigned int, unsigned int>>();
  else
    {
      DEAL_II_NOT_IMPLEMENTED();
      return std::vector<std::pair<unsigned int, unsigned int>>();
    }
}

(see https://github.com/dealii/dealii/blob/2170eda0660b33e988d89a5816279387702db0dd/source/fe/fe_dgp.cc#L183-L197)

In the FE_DGQ case we never hit this because deal.II simply decides not to throw and returns an empty vector. The replacement I proposed in this PR does the same.

[luca-heltai]

Maybe I would stick to the name FE_AggloDGP following the notation of the other deal.II FiniteElement classes?

[fdrmrc]

Maybe I would stick to the name FE_AggloDGP following the notation of the other deal.II FiniteElement classes?

You're right, thanks. I did as suggested. However, I think the present assert upstream should not be triggered at all. The documentation for that function is the same as the one in the FE_DGQ class, and in that case an empty vector is returned. Do you agree?

[luca-heltai]

Yes, but then why is not working upstream?

[fdrmrc]

Sorry, I'm not sure I understood your question. DGP elements work upstream in deal.II. However, with the present deal.II master, we cannot use DGP elements in our agglomeration codes (in debug mode) because the assert in the next snippet is triggered.

//dealii/source/fe_dgp.cc
template <int dim, int spacedim>
std::vector<std::pair<unsigned int, unsigned int>>
FE_DGP<dim, spacedim>::hp_vertex_dof_identities(
  const FiniteElement<dim, spacedim> &fe_other) const
{
  // there are no such constraints for DGP elements at all
  if (dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe_other) != nullptr)
    return std::vector<std::pair<unsigned int, unsigned int>>();
  else
    {
      DEAL_II_NOT_IMPLEMENTED();
      return std::vector<std::pair<unsigned int, unsigned int>>();
    }
}

For DGQ elements, the implementation does not throw and returns the empty vector. See next snippet.

//dealii/source/fe_dgq.cc
template <int dim, int spacedim>
std::vector<std::pair<unsigned int, unsigned int>>
FE_DGQ<dim, spacedim>::hp_vertex_dof_identities(
  const FiniteElement<dim, spacedim> &fe_other) const
{
    return std::vector<std::pair<unsigned int, unsigned int>>();
}

I find these different implementations of the same function inconsistent because DGP and DGQ elements should behave in the same way.

The documentation for the DGP version (https://www.dealii.org/developer/doxygen/deal.II/classFE__DGP.html#a6937c42e4dd88f4d4f808559d6b49890) is indeed identical to the one of the DGQ case (https://www.dealii.org/developer/doxygen/deal.II/classFE__DGQ.html#aab33a2ac3dd777ae00f00fd7f4b18cb7)

[luca-heltai]

I think that the assertion is a bug in deal.II.

DG elements do not define dofs on faces, so they should always return the empty vector. Can you submit a patch removing the assert in deal.II?

I think this is a leftover for copying the DGP code from continuous P code.

[fdrmrc]

DG elements do not define dofs on faces, so they should always return the empty vector. Can you submit a patch removing the assert in deal.II?

Sure, I'll do that.