stadic condensation for LDG

[Guosheng Fu]

Hello,

I am using LDG to solve the Laplace equation.
I want to do stadic condensation to eliminate the flux variable resulting a global linear system for the scalar only.

To do so, I manually set scalar dofs to be interface_dof,

for el in V.components[1].Elements():
    for i in range(len(el.dofs)):
          V.SetCouplingType(ndof+el.dofs[i],COUPLING_TYPE.INTERFACE_DOF)
          V.FreeDofs()[ndof+el.dofs[i]] = True
          V.FreeDofs(True)[ndof+el.dofs[i]] = True

but it seems this makes stadic-condensed matrix to be singular...

The sample code is put in the following post.

Best,
Guosheng

[Anonymous]

[Anonymous]

[Anonymous]

[Guosheng Fu]

the code

[schruste]

Dear Guosheng,

What you are trying to do is not supported in NGSolve right now. The static condensation only works for (completely) element-based assemblies. Only after the element-volume and element_boundary (cf. HDG formulations) contributions of one element the static condensation is executed. In your case that means that the skeleton-integrals are not treated correctly which leads to singular matrices.

If you have a formulation suited to the usual structure -- formulated by this type of element-volume - element_boundary-type integrals -- you are safe. Note: element_boundary-type integrals are "volume"-integrals internally in NGSolve and thus static condensation behaves as for the volume-terms.

When using DG (dgjumps) or -- more importantly -- skeleton-type integrals the structure is different. You would need to compute the facet integrals for all facets of one element first in order to do that and exploit the a priori knowledge that there is no coupling between the LOCAL_DOFs of neighboring elements. This is not provided in NGSolve right now.

You would need to reformulate your discretization in terms of element_boundary-type integrals in order to apply static condensation.

This is much easier with Hybrid methods as the element_boundary-integrals are more natural. Further, the dependencies to neighboring elements are gone and the CouplingTypes are clear a priori...

If you don't want to change the formulation, what you can also do is the following: Assemble the full system (without static condensation) and do the static condensation by hand afterwards (e.g. with scipy, try to exploit the block structure.... ).

Hope this explains NGSolve's behavior and that you can work with that.

Best,
Christoph

[Guosheng Fu]

Hi Christoph,

I see, it make sense now.

Thanks,
Guosheng