Linear constraints

[lijas]

It would be nice to have the possibility to add linear constraints. Linear constraints are described here, and deal ii has a good tutorial with it.

To apply the linear constraints, one can condense the stiffness matrix via: C^T * K * C. However this is pretty inefficient since it creates a new (condensed) stiffness matrix. The way deal ii does it, is that they call a function which adds the correct stiffness terms to the masterdofs in the K-matrix, and puts 1s on the diagonal for the dofs that are constrained (similar to how our apply!(ch, K, f) function works). After solving the system, they then distribute the correct values to the constraint dofs. If the number of linear constraint is low (<<ndofs), this should be much faster.

ch = ConstraintHandler(dh)
add(ch, LinearConstraint(3, [1.0, 2.0], [4,5]); # u3 = 1.0*u4 + 2.0*u5
close!(ch)
...
K = create_sparsity_pattern(dh, ch) #also give info about constraind dofs
...

apply!(ch, K, f)
a = K\f
distribute!(a, ch) # or apply!(...)

I have tried creating the function that does the C^T * K * C operation (but without creating a new matrix), but it is a bit trickier than I thought. It would be nice the get some ideas/help here.

[termi-official]

Maxi and I had a prototype of assembling C explicitly. Maybe @koehlerson can drop a reference to the repository since I cannot find it. We used [1] to implement it, but care, this is very different from how deal ii handles it.

In numerical analysis the operation you describe is sometimes called RAP-operator (for restriction-action-prolongation). Basically what you can do here is to define a matrix type which carries C and K along with its action (i.e. implementing the multiplication function [2]). This way you can apply iterative solvers without explicitly constructing C^T K C. However it is ongoing research to construct efficient preconditioners for these kind of problems.

I have a question here: Do you by any chance know who came up with what is called "The Master-Slave Method" in your reference? Turns out I worked on basically the same method for some time. :)

References [1] Cerveny, Jakub, Veselin Dobrev, and Tzanio Kolev. "Nonconforming mesh refinement for high-order finite elements." SIAM Journal on Scientific Computing 41.4 (2019): C367-C392. [2] URL.https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.mul!

[dkarrasch]

A walk-by comment: For lazy matrix composition, you may want to take a look at LinearMaps.jl. Following the implemented types there, you might consider subtyping LinearMap there to get all the functionality for free, or just use the existing LinearMap(C)'K*C::CompositeMap approach.

[lijas]

Thanks for the feedback. I already opened a PR for linear constraints in #401 . In this, I used the same algorithm as in deal ii. Basically it does the C'KC operation directly on K, so should be quite efficient if the number of constraints is "much" less than the number of dofs.