Partial Static Condensation


I am trying to test the preconditioner in the paper ([1801.04707] Preconditioning of a hybridized discontinuous Galerkin finite element method for the Stokes equations , Preconditioning of a Hybridized Discontinuous Galerkin Finite Element Method for the Stokes Equations - Journal of Scientific Computing). There, only the momentum block is condensed. I am trying to replicate the same by changing CouplingType of the local_dof corresponding to pressure unknowns. However, when I solve the resulting system, I get a solution which does not make sense. Is this the intended way? Or should I be doing it in a more manual way using multiple bilinear forms?

Best regards,
Abdullah Ali Sivas

Hi Abdullah,
having hybrid variables for the pressure you should be able to condense also the element-wise pressure. Yes, you do so by setting the coupling-type.

for debugging, you can dump and check element-matrices as follows :

SetTestoutFile (“test.out”)
bf = BilinearForm(…, printelmat=True)


Hi Joachim,

Thanks for your reply. I now know that I am on the correct route! Also printelmat will definitely be useful for debugging purposes.

However, what I mean is if I try to solve the system described in Eqs. (54)-(55), i.e. eliminating both element momentum and elemnent pressure, -for test purposes using a direct solver- I get the correct solution. However, if I try to solve the system as described in Eqs. (65)-(66) (by setting the coupling type of element pressure unknowns to interface dofs for example), I get a solution which does not make sense.

Should I share my code? Would it better that way?

Best wishes,
Abdullah Ali Sivas

Hi Abdullah,
if you attach a small test example we can have a look at it.

Hi Christopher,

Attached I am sending an implementation of the method. I am using the variable eliminate_up to switch between condensing both u and p and condensing only u. I am setting the coupling type in lines 74-76. I appreciate your time.

Thanks a lot,
Abdullah Ali Sivas


Hi Abdullah,

after changing the coupling-types, you have to update internal tables by calling X.FinalizeUpdate()

Then both versions produce the same result,

best, Joachim

Awesome! Just tried, thanks a lot for your time! I am new to NGSolve and I didn’t notice that function even existed.

Abdullah Ali Sivas