Make Gauss-Lobatto points the default support points in FE_Q/FE_DGQ

[kronbichler]

As expected the change in polynomials at degree three and higher resulted in a notable amount of changes in the output. In particular, the output of shape functions, prolongation & restriction matrices etc. resulted in a sizable change. When reviewing the patch, I recommend to skip the patch 'Adjust test output due to change to Gauss-Lobatto points' which is rather mechanical.

I also have a few changes in the test suite where I switched the old Gauss-Lobatto-based FE_Q versions to the default constructor. A notable change is in fe/restrict_prolongation_01.cc where I fixed a bug: The old version would skip a check done done in DoFAccessor::get_interpolated_dof_values, namely not overwriting when the incoming value is zero. This did not happen for equidistant points of FE_Q where there each mother DoF is also on some of the children (and all others are thus zero) but it does happen for non-equidistant points (in other words, the test would have been broken for arbitrary node distributions).

As promised, I wrote documentation in FE_Q and FE_DGQ about this change.

This fixes #2438.

[davydden]

an option could be to explicitly specify equidistant support points in unit tests so that you would not need to change their output and therefore be backward compatible.

[kronbichler]

@davydden That would be a possibility. I added the equidistant points in some of the tests (see e.g. my changes in tests/fe/fe_data_test.cc, tests/fe/fe_tools_01.cc) but I think that testing the default more thoroughly is not wrong either. I know it changes many things (due to copying the new result files also some roundoff changes crept into my patch). The largest change sets are in hp (create_{laplace,mass}_matrix, hp_constraints) where large matrices are output. There it might be worth to change the element type to reduce the change set - a plain text version of the diff file of my last commit is 86M. :-( Any opinions?

[kronbichler]

Just to clarify my intent with testing GL points in a broader context: After my change I discovered a bug in fe/restrict_prolongation_01.cc where I first thought that there was a problem in the restriction matrix of FE_Q on non-equidistant points and only after a few hours I realized that the bug was in the test.

[kronbichler]

This supersedes #2436.

[davydden]

I think that testing the default more thoroughly is not wrong either

fully agree. I would just keep those for future patches to keep diffs reasonable for a single PR. Then you would add a proper fe/restrict_prolongation_01.cc which does the right thing for GL.

In any case, let's see what other guys think...

[bangerth]

I think I would have gone with @davydden 's suggestion of just changing the tests to explicitly use equidistant points, just to keep the diffs small. In fact, one might argue that testing the uncommon path (which is now to use equidistant points) is the better approach because that path is not going to be tested in other codes and/or future tests.

Would it make sense to do something like this:

  class FE_Q
  {
     enum BasisChoice { Gauss_Lobatto, Gauss_Legendre };
     FE_Q (const unsigned int q,
                const BasisChoice basis_choice = ...);

to make selecting one or the other one of the common bases simpler? One can still choose a completely different basis by calling the other constructor that takes a set of quadrature points.

[kronbichler]

I could definitely change things in the hp tests which reduces the changes in tests by 90% to judge from git diff statements. I think I will leave the smaller things as they are but the larger change sets will remain the old behavior. This way we keep a reasonably tested variant for the equidistant points.

Today I'm in a hurry so I'll do it during the next days (unless someone else comes with better suggestions).

When it comes to be basis choice: From a user's perspective, one can write then FE_Q<dim,spacedim>(q,FE_Q<dim,spacedim>::equidistant); rather than FE_Q<dim,spacedim>(QIterated<1>(QTrapez<1>(),q)). Sure, the former is more explicit and clear in what it does. On the other hand, the advantage is not really large given that I also have to take this into account in the get_name() method. I could probably add the latter to the output whenever I detect equidistant points and teach fe_tools to read this. But I'm open to do that if we see enough value in providing it. (For the test suite it's easy enough to insert Qiterated<1>(QTrapez<1>() through a script.

[bangerth]

On 04/06/2016 08:19 AM, Martin Kronbichler wrote:

When it comes to be basis choice: From a user's perspective, one can write
then |FE_Q<dim,spacedim>(q,FE_Q<dim,spacedim>::equidistant);| rather than
|FE_Q<dim,spacedim>(QIterated<1>(QTrapez<1>(),q))|. Sure, the former is more
explicit and clear in what it does. On the other hand, the advantage is not
really large given that I also have to take this into account in the
|get_name()| method. I could probably add the latter to the output whenever I
detect equidistant points and teach fe_tools to read this. But I'm open to do
that if we see enough value in providing it. (For the test suite it's easy
enough to insert |Qiterated<1>(QTrapez<1>()| through a script.
Either way is fine by me. It was just an idea. Maybe documenting the way to
get equidistant points in the docs of the constructor is enough.

[davydden]

to keep PR smaller, we could also leave the idea of enum to a later PR, if it will be decided that it's a good feature to have.

[guidokanschat]

My favorite is two classes: FE_Q for Gauss-Lobatto and FE_Q_Equidistant for the old one. Both only have a constructor and use a common base class.

[kronbichler]

In a very long afternoon, I changed quite a number of tests back to equidistant points in order to reduce the number of changes. It is now down from over 1 million insertions/deletions to around 28k insertions and deletions. For some tests it was not feasible to convert to equidistant points due to construction (fe/fe_restriction_xxx, fe/fe_prolongation_xxx, fe/fe_support_points_xxx, bits/dof_tools, bits/fe_tools) while I did it for many others, in particular the hp folder. This should cover both cases (equidistant points and Gauss-Lobatto points) very well both for FE_Q and FE_DGQ.

[kronbichler]

@guidokanschat That might be a good idea - but let's leave that to another PR. (There was enough manual work in this one already.)

[kronbichler]

I suggest to wait for a few hours. I'm re-running the test suite right now.

[kronbichler]

Except for the problems with manifold reported in #2491 this patch passes the testsuite.