step-20 Mixed Poisson with different quadrature rule

[Eldar Khattatov]

Hello,

I am relatively new to Deal.II and was trying to implement the mixed finite element method for the Darcy equation in 2d using BDM1 x Q0 spaces. The only non-standard thing that I am trying to do is to use the trapezoid quadrature rule for the LHS terms in variational formulation, actually I need this quadrature on the (u,v) term only, but for simplicity I used it on the (p, div(v)) term as well, so I didn't do also any changes to the code in step-20 example except this:

    QGauss<dim> quadrature_formula(degree+2);

changed to 

    QTrapez<dim> quadrature_formula;

in assemble_system() implementation.

Due to a paper by Wheeler and Yotov this method shows the 1st order of convergence in the velocity, however running the step-20 code in Deal.II with the only change shown above leads to 0.5 order of convergence in the velocity.

Could you please point to my mistake in the implementation of the desired change in quadrature rule?

Thank you in advance

[Eldar Khattatov]

I might have provided too little details, so that no one replied, but I am still looking for help with this question. 

I am using BDM1 - Q0 pair of finite elements for velocity-pressure variables. By now I have implemented the use of trapezoid quadrature rule in a(*,*) term in the variational formulation only, while the rest of the bilinear forms are integrated using Gaussian quadrature rule of order 3 (I tried >= 3 with no improvement). The terms in RHS involving the integral coming from integration by parts (imposing pressure BC) are also integrated using Gaussian quad. rule of order 3, so that the implementation agrees with the method discussed here: http://www.math.pitt.edu/~yotov/research/publications/mfmfe.pdf

Except for the changes mentioned above, the code I use is the one from step-20.cc taken from deal.II website. 

The visualized solutions show that the pressure is being computed correctly (same can be seen from errors and rates in pressure variable) while both components of the velocity show instabilities near the boundary as seen on the pictures attached to this post. I cannot see where these instabilities come from, as my change of quadrature rule does not (?) affect the way boundary conditions are imposed. 

Any ideas would be much appreciated, thank you.

[Eldar Khattatov]

I tried running the code for the test case with homogeneous pressure boundary conditions p = 0 on the boundary. This lead me to the correct solutions with theoretically predicted rates and without any instabilities in velocity. So it seems that the way I am changing the quadrature rule is somehow affecting the way natural boundary conditions are imposed, unfortunately I can not see how. I attached the code to this post if someone could take a look at the changes in the assemble_system() that I made and point me to my mistake. Thank you.

[Eldar Khattatov]

Turns out there is indeed a necessity to treat the natural bc in a special way in this method, this was discussed in the other paper by same authors, so it is not an implementation or Deal.II issue.