Step-44: Use FE_DGP instead of FE_DGPMonomial?

[jppelteret]

After a discussion with @peterrum, I was compelled to try substituting the FE_DGPMonomial used for step-44's pressure and dilatation fields by FE_DGP. After re-reading the documentation for these two classes, it appeared to me that FE_DGP should work just fine. In order to satisfy the LBB condition, the requirement of the fields is that they are discontinuous, of one degree lower than the displacement field, and contain the complete polynomials of the FE space (see Hughes "The Finite Element Method" p295/296). For example, in 2d a bi-quadratic displacement must be paired with a piece-wise linear (not bilinear) pressure and dilatation. Both the FE_DGP and FE_DGPMonomial classes claim to do this, with the main difference between them being the following:

The basis functions for this element [FE_DGPMonomial] are chosen to be the monomials listed above. Note that this is the main difference to the FE_DGP class that uses a set of polynomials of complete degree p that form a Legendre basis on the unit square. Thus, there, the mass matrix is diagonal, if the grid cells are parallelograms. The basis here does not have this property; however, it is simpler to compute. On the other hand, this element has the additional disadvantage that the local cell matrices usually have a worse condition number than the ones originating from the FE_DGP element.

So it would seem that we can expect both to work, and should observe a difference in condition number of the linear system. I've tested this with order 4 displacement (and corresponding order 3 pressure and displacement fields) and this does seem to be the case. The FE_DGP seem to give superior performance.
step-44-Q(4)_DGP(3)_DGP(3).txt
step-44-Q(4)_DGPMonomial(3)_DGPMonomial(3).txt

Unfortunately, I cannot recall why we chose the FE_DGPMonomial element as opposed to the FE_DGP. @mac-a, do you have any recollection as to why we did this? I propose that we investigate changing the element type for these discontinuous fields, on the basis of an increase in performance with no apparent decrease in the accuracy of the element.

step-44.cc.txt

[peterrum]

FYI @mschreter

[mac-a]

@jppelteret I have no recollection why we chose FE_DGPMonomial. We did explore both options at the time though. I agree though, it seems it could improve conditioning. Worth investigating.

[jppelteret]

We did explore both options at the time though.

Thanks for your reply! I recall the same, but I thought I might be imagining it! Maybe there was some confounding error in the code when we did so. I'll do a small parametric study to see if such a change would be robust.

[mac-a]

I seem to recall that we had an issue setting J = 1 as the initial condition for the dilatation.

[bangerth]

I would expect FE_DGP to lead to a better condition number. The fact that you see that experimentally makes great sense to me.

[jppelteret]

I've attached some logs below for tests using poly orders 1-4 with the grid refinement level decreased as p increases. There's no change in the number of linear solver iterations for p=1, but a significant improvement at p=4.

DGP_Monomial
step-44-dgp_monomial-P1_Q2_GR_4.txt
step-44-dgp_monomial-P2_Q3_GR_3.txt
step-44-dgp_monomial-P3_Q4_GR_2.txt
step-44-dgp_monomial-P4_Q5_GR_2.txt

DGP
step-44-dgp-P1_Q2_GR_4.txt
step-44-dgp-P2_Q3_GR_3.txt
step-44-dgp-P3_Q4_GR_2.txt
step-44-dgp-P4_Q5_GR_2.txt