Wrong sparsity pattern for vector-valued dG (IPDG) methods with fluxes over faces with (FE_DGQ) and deal.II v9.0

[Uwe Köcher]

Dear all,

I'm facing a (new?) problem for deal.II v9.0.x by getting an obviously wrong sparsity pattern for IPDG methods

for vector valued problems (i.e. linear elasticity, elastic wave propagation etc.).

Fortunately, I'll get an internal Assert coming from a Trilinos function whenever I want to distribute my face terms into the global matrix with flux terms.

Please note that I could not check the current developer version for the this issue so far.

I'll open an issue on github soonly and this thread should support this by having the opportunity to include figures.

History: for certain reasons I usually use FE_DGQArbitraryNodes along with a Q-Gauss

quadrature, so maybe the reported issue is not really new, but I think this comes from

v8.5.1 to v9.0.x

Problem: the flux sparsity pattern is incorrect, if quadratures having 0 and 1 in the list of q-points.

Note that this thread reads as a top-level problem description without studying deeper the deal.II library where a (or the) error(s) occur.

The correct pattern for a locally polynomial degree p=1 (using QGauss(2)) on a tria with 4 elements for dim=spacedim=2,

without reordering by its components, looks like this:

,

where the figure was produced from and FE_System having 2 FE's of the Type dealii::FE_DGQArbitraryNodes<dim>(*fe_quad)

along with a Q-Gauss(2) quadrature and writing out the sparsity pattern;

but, whenever I use Q-GaussLobatto(2),

or (even harder to understand) a user defined quadrature having at least the q-points 0 and 1,

I'll get a sparsity pattern with missing entries:

.

Please note that I'm using a coupling table with dealii::DoFTools::nonzero and dealii::DoFTools::always.

Also note: It looks like that DGQ in a FESystem only couples the dof of the first vector component and not all vector components

(which is barely necessary for elasticity).

A test code can be found on:

https://bitbucket.org/dtmproject/unittest-003

This code is not minimal, but minimized to find out what may goes wrong.

Notes on my current studies:

If a quadrature (having 0 and 1 in the set of q-points), the code

dealii::FE_DGQArbitraryNodes<dim>(*fe_quad)

forces deal.II v9.0 to switch over to an FE_DGQ Element,

even if the local polynomial degree is higher than 2

and may lead to the problem.

Further notes to the developers: (from my understanding of the former releases)

- when I remember correctly, this issue does not occur with v8.5.1 and/or v8.4.x

- the FE_DGQ( QGaussLobatto(4) ) should produce (slightly) different basis functions compared to

the FE_DGQ(3) (having equidistantly spaced local dofs including "face" dof)

- I'm not sure how long the performance switch of DGQ_Arbitrary -> standard DGQ is inside the library

Next week I will have the opportunity to study the problem deeper in the library

(i.e. testing it with the developer version, and older versions,

maybe finding the problem and implementing a patch and test suite programs)

Kind regards

  Uwe

[Uwe Köcher]

as proposed, I've opened an issue on github

https://github.com/dealii/dealii/issues/7390