FE_DGP/FE_FaceP and hp::DoFHandler

[sml.i...@gmail.com]

Dear all,

As indicated in the documentation, the finite elements FE_DGP and FE_FaceP are not based on nodal interpolation (as eg. FE_DGQ and FE_FaceQ), but on projection. Apparently, this prevents these elements from being used in combination with a hp::DoFHandler, since there seems to be no easy way (or even no way?) to implement the functions hp_*_identities for FE_FaceP (which is needed to support hp) when no nodal interpolation (ie. no unit_support_points) is available.

I would really like to use the FE_DGP/FE_FaceP combination in a hp FEM computation, as this requires fewer DoFs than the FE_DGQ/FE_FaceQ couple.

Does anyone have an idea how to fix this problem? I imagine that either:
- one could implement hp_*_identities, somehow, based on the existing code with projection, or
- one could implement a different version of FE_FaceP with nodal interpolation
but i don't really know whether these approaches could work and where to start.

Best,
Samuel

[Wolfgang Bangerth]

> As indicated in the documentation
> <http://dealii.org/developer/doxygen/deal.II/classFE__DGP.html>, the finite

The interpolation property is not actually necessary. What is necessary for
using things in the hp framework is that the space on one side is *larger*
than the one on the other side, and consequently that you can constrain one to
the other by a set of linear constraints. That is what
make_hanging_node_constraints() does.

The hp_*_identities actually do a subset of that, namely identify degrees of
freedom on one side with those on the other side. For interpolation, that is
simple. But it doesn't have to be. Think of shape functions
1, x, x^2
along an edge (-1,1), for example. These are moments, not interpolatory. Now
on a neighboring cell, you may have
1, x, x^2, x^3
and for continuity, you'll have to have that on that cell, the coefficient for
the last shape function is zero and the coefficients for the first three equal
the same coefficients for the three shape functions on the first of the two
cells. This identity is what the hp_*_identities() functions record. In other
words, for that face, you should get the identities
(0,0) // DoF 0 on one side equals DoF 0 on the other side
(1,1)
(2,2)
and later on, make_hanging_node_constraints() would set U_3 on the second
cell's face to zero.

Does that make sense?

Best
W.


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

[sml.i...@gmail.com]

Yes, that makes sense to me. But now if we have neighboring faces with differrent refinement levels (e.g. the face element on one side is only half as long as the one on the other side) the situation gets a little more complicated, right? Since then you cannot just copy the coefficients 1 to 1, or do I get that wrong?

So what you suggest is that one could actually go on and implement hp_*_identities for FE_FaceP, in the way you describe? Do you know of some other element that already implements this, not using interpolation, which one could use as reference?

Best,
Samuel

[Wolfgang Bangerth]

On 08/04/2017 07:07 AM, sml.i...@gmail.com wrote:
> Yes, that makes sense to me. But now if we have neighboring faces with
> differrent refinement levels (e.g. the face element on one side is only half
> as long as the one on the other side) the situation gets a little more
> complicated, right? Since then you cannot just copy the coefficients 1 to 1,
> or do I get that wrong?

Yes. But then make_hanging_node_constraints() knows what to do with this. The
hp_*_identities() functions aren't used for this.

> So what you suggest is that one could actually go on and implement
> hp_*_identities for FE_FaceP, in the way you describe? Do you know of some
> other element that already implements this, not using interpolation, which one
> could use as reference?

I don't know for sure if that element implements the hp constraints, but the
modal RT element has a similar structure in that the degrees of freedom are
located on the faces and are defined via modes or integrals, rather than
interpolation. (There is also an interpolatory RT element; both of them are
available in deal.II.)

[sml.i...@gmail.com]

I have tried the following implementation (for 2D computations only hp_line_dof_identities has nonempty constraints):

template <int dim, int spacedim>
std::vector<std::pair<unsigned int, unsigned int> >
FE_FaceP<dim, spacedim>::
hp_line_dof_identities (const FiniteElement<dim, spacedim> &fe_other) const
{
  // this element is continuous only for the highest dimensional bounding object
  if (dim !=2)
    return
      std::vector<std::pair<unsigned int, unsigned int> > ();
  else
    {
      if (const FE_FaceP<dim,spacedim> *fe_p_other
          = dynamic_cast<const FE_FaceP<dim,spacedim>*>(&fe_other))
        {
          std::vector<std::pair<unsigned int, unsigned int> > identities;

          for (unsigned int i=0; i<std::min(this->dofs_per_face, fe_p_other->dofs_per_face); ++i)
            identities.emplace_back (i, i);

          return identities;
        }
      else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != nullptr)
        {
          // the FE_Nothing has no degrees of freedom, so there are no
          // equivalencies to be recorded
          return std::vector<std::pair<unsigned int, unsigned int> > ();
        }
      else if (fe_other.dofs_per_face == 0)
        {
          // if the other element has no elements on faces at all,
          // then it would be impossible to enforce any kind of
          // continuity even if we knew exactly what kind of element
          // we have -- simply because the other element declares
          // that it is discontinuous because it has no DoFs on
          // its faces. in that case, just state that we have no
          // constraints to declare
          return std::vector<std::pair<unsigned int, unsigned int> > ();
        }
      else
        {
          Assert (false, ExcNotImplemented());
          return std::vector<std::pair<unsigned int, unsigned int> > ();
        }
    }
}

so it just adds identities up to the minimum number of shape functions involved. But somehow it doesn't work as expected (in my model computation I still see the skeleton solution not agreeing on both sides of faces with hanging nodes). Do you see some apparent mistake? Am I wrong in assuming that the shape functions are ordered by ascending polynomial degrees?

I have appended the output of a testcase (I just adapted the one which D. Arndt wrote for FE_FaceQ here: https://github.com/dealii/dealii/pull/4577)

[sml.i...@gmail.com]

Can anyone help me with that matter?

Best, Samuel