Add global dofs

[KnutAM]

Solves #1009 following @lijas's suggestion.

Currently the new functionality is not exported, and I'm considering if it makes sense to add this to devdocs until we battle-test it a bit before committing to this interface.

[termi-official]

Shouldn't celldofs[!] return these additional dofs, too?

[KnutAM]

No, the global dofs in this PR have no coupling with the cells, and this can be added for example via AffineConstraints.

For "the other type of global dofs", I think we should go with something like master...kam/GlobalDofs, which allows you to do

dh = DofHandler(grid)
sdh1 = SubDofHandler(dh, domain1)
add!(sdh1, :u, Lagrange{RefTriangle,1}())
add!(sdh1, :g, GlobalInterpolation(n)) # Add n extra dofs that are connected to all cells in sdh1
sdh2 = SubDofHandler(dh, domain2)
add!(sdh2, :u, Lagrange{RefTriangle,1}())
close!(dh)

If you only want the global dof on some subdofhandlers.

[termi-official]

[...] the global dofs in this PR have no coupling with the cells [...]

Okay, then I do not understand what the purpose of feature. If you have an independent dof, then you have effectively two independent systems. Hence we can just have a separate dof handler with the additional dofs and assemble two separate systems. Can you point me to the detail which I am missing out?

[...] this can be added for example via AffineConstraints.

Not sure how this fits in the statement above. From what I originally understood in the related discussion in the linked issue is that you just want to have some mechanism to add "simple constraints". Why should the dof handler manage this and why can't we add have some LagrangeMultiplierConstraint which is managed by the constraint handler instead? Or what exactly is the advantage of this design, having the dof handler managing these constraints (partially)?

[KnutAM]

Okay, then I do not understand what the purpose of feature. If you have an independent dof, then you have effectively two independent systems.

By having the DofHandler aware of the total number of dofs, the constraint handler can rely on ndofs(dh) to return the number of rows and columns in the system matrix to which it should add constraints. In #1009 I first proposed passing this info directly to the ConstraintHandler, but this PR introduces the alternative solution that passes this information via the DofHandler. With this design, one can do

dh = DofHandler(grid)
add!(dh, ...)
add_global_dofs!(dh, ...)
close!(dh)
ch = ConstraintHandler(dh) ...
add!(ch, ...)
close!(ch)
K = allocate_matrix(dh, ch)

whereas with my original suggestion you need

dh = DofHandler(grid)
add!(dh, ...)
close!(dh)
n, global_dof_info = setup_global_dofs(...) # User defined data structures
ndofs_total = ndofs(dh) + n
ch = ConstraintHandler(dh; ndofs_total) # New feature that we would need to add. 
add!(ch, ...)
close!(ch)
sp = SparsityPattern(ndofs_total, ndofs_total)
add_sparsity_entries!(sp, dh, ch)
K = allocate_matrix(sp)

Or what exactly is the advantage of this design, having the dof handler managing these constraints (partially)?

The DofHandler doesn't really manage these constraints, it is just aware of the extra dofs. Currently, my workaround is doing dh.ndofs += n after close!(dh), which tricks the setup to work, but provides no info about the additional dofs, this I have to manage myself. The main advantage with adding fields to the dofhandler is IMO that you can directly request these from the DofHandler, and that ndofs(dh) gives the correct number of the degrees of freedom in the system.

[KnutAM]

Why should the dof handler manage this and why can't we add have some LagrangeMultiplierConstraint which is managed by the constraint handler instead?

It might also be possible to do

dh = DofHandler(grid)
add!(dh, ...)
add_global_dofs!(dh, ...)
close!(dh)
ch = ConstraintHandler(dh) ...
add!(ch, ...)
add!(ch, LagrangeMultiplierConstraint(:lambda, n, AffineConstraintInfos)) # Extends with n number of dofs
close!(ch)
K = allocate_matrix(dh, ch)

One advantage that I see with this setup is that it potentially could make it easier to "deactivate" some multipliers, keeping the sparsity pattern the same (nice for e.g. contact simulations).

A disadvantage I see is that it is less flexible then the suggestions above.

[termi-official]

So we are really talking about just adding some constraint $\sum_i u_i = 0$? Or why can you skip the setup_global_dofs to build the connectivity in the the examples 1. and 3.?

A disadvantage I see is that it is less flexible then the suggestions above.

I see this as an advantage, because is is clear what to do in all scenarios. Whereas the in proposed design it seems heavily underdesigned, as it is absolutely unclear to my how the interaction between sparsity pattern, constraints, assembly and distributed parallelization should be (and it is also not clear from going over the code). So maybe it is easier to have a special dof handler for such specific corner cases?

[KnutAM]

So we are really talking about just adding some constraint $\sum_i u_i = 0$?

Yes, the case I hit was adding an extra dof. Then I add many AffineConstraints which constrain other dofs to this extra dof.
(The reason for having it in the system, as opposed to fixing it by hard-coding it in the affine constraints is that I need this for computing correct residuals and sensitivities later, and there are cases when I don't constrain it, but still want the AffineConstraints to the extra dof).

Whereas the in proposed design it seems heavily underdesigned, as it is absolutely unclear to my how the interaction between sparsity pattern, constraints, assembly and distributed parallelization should be (and it is also not clear from going over the code).

The only difference apart from the possibility to get information from the DofHandler, is that it increases the value of ndofs(dh), implying that the total equation system is larger. This allows you to add constraints to these extra dofs, which then will affect e.g. the sparsity system as always. It does not change anything for the assembly. For distributed parallelization it should have no effect either, the effect comes if you tie two dofs together via Affine constraints, which is not different with this PR from before.

So to summarize, an absolute minimal change with this pr would be to add a field n_global_dofs::Int to the dofhandler (allowing you to add a given number of extra dofs, but without any naming etc. of these), and redefine ndofs(dh) = dh.ndofs + dh.n_global_dofs (or update dh.ndofs to this value during close!.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.22%. Comparing base (2c8b751) to head (6b8029f).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1012      +/-   ##
==========================================
+ Coverage   93.20%   93.22%   +0.01%     
==========================================
  Files          38       38              
  Lines        5816     5832      +16     
==========================================
+ Hits         5421     5437      +16     
  Misses        395      395              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.