Merging DofHandler and MixedDofHandler

[fredrikekre]

Background

Ferrite.jl has, for quite some time, maintained two different DoF handlers: the original DofHandler and the MixedDofHandler. The DofHandler is limited to i) single element type (e.g. only triangles, or only quads) and ii) to uniform dof distribution (e.g. all fields live on all elements). MixedDofHandler enables both these use cases, so functionality-wise it it is a superset of DofHandler. Naturally it would be good to only have a single DoF handler instead of maintaining two.

Performance wise the DofHandler is generally better, since it doesn't have to deal with the Union element type (when combining different elements). A lot of work, mainly by @kimauth has been done to reduce the performance gap, and I think that we are now at a point where we can merge them into a new DoF handler that can support all usecases with good performance.

DofHandler interface

Setup

dh = DofHandler(grid)
add!(dh, :u, 2, u_interp)
add!(dh, :p, 1, p_interp)
close!(dh)

Constraints

ch = ConstraintHandler(dh)
add!(ch, Dirichlet(...))
close!(ch)

Assembly loop

for c in CellIterator(dh)
    # ...
end

MixedDofHandler interface

(I have never used this, so please correct me if the stuff below is incorrect or atypical.)

Setup

dh = MixedDofHandler(grid)
# Domain 1: u
fh1 = FieldHandler([Field(:u, u_interp, 2)], set1)
add!(dh, fh1)
# Domain 2: u and p
fh2 = FieldHandler([Field(:u, u_interp, 2), Field(:p, p_interp, 1)], set2)
add!(dh, fh2)
close!(dh)

Constraints

ch = ConstraintHandler(dh)
add!(ch, fh1, Dirichlet(...)) # add!(ch, Dirichlet(...)) also exist now though
close!(ch)

The add!(ch, fh, dbc) method is a bit strange, since you can already specify the domain selection in the Dirichlet constructor. There is no need to be able to do it twice, IMO.

Assembly loop

# Domain 1
for c in CellIterator(dh, fh1.cellset)
    # ...
end
# Domain 2
for c in CellIterator(dh, fh2.cellset)
    # ...
end

NewDofHandler interface

These are my ideas for a new common interface that I have also briefly discussed with @kimauth and @lijas

Setup

dh = NewDofHandler(grid)
# Domain 1: u
sdh1 = SubDofHandler(grid, set1)
add!(sdh1, :u)
add!(dh, sdh1)
# Domain 2: u and p
sdh2 = SubDofHandler(grid, set2)
add!(sdh2, :u)
add!(sdh2, :p)
add!(dh, sdh2)
close!(dh)

Constraints

ch = ConstraintHandler(dh)
add!(ch, Dirichlet(...))
close!(ch)

Assembly loop

# Domain 1
for c in CellIterator(sdh1)
    # ...
end
# Domain 2
for c in CellIterator(sdh2)
    # ...
end

This is mostly the same as the current MixedDofHandler, with FieldHandler now called SubDofHandler. I don't really like the name FieldHandler because what it really is is a "region-handler". I think it makes sense to rename, and then also include the parent dofhandler as a field. Then the subdofhandlers can be drop-in replacement in most cases. In addition, we can put the concrete celltype as a type-paramter on the subdofhandler, which we can use to help the compiler when requestion cells from the union-vector.

In addition (not shown above) the exact same interface as for the current DofHandler will be supported. In that case, when adding the first field, an internal SubDofHandler will be set up. In this mode, e.g. dof_range will work directly on the dofhandler, but it will error in the previous case (there you will be forced to call it on your subdofhandler).

Implementation plan

Initially the plan was to simply remove DofHandler and then rename MixedDofHandler to DofHandler. However, I think the following might be easier to execute:

• "Fork" NewDofHandler from MixedDofHandler
• Update NewDofHandler to new common interface
• Compare and test performance and functionality between NewDofHandler and DofHandler/MixedDofHandler
• Delete DofHandler and MixedDofHandler, rename NewDofHandler to DofHandler.

[fredrikekre]

Okay, now you may read.

[lijas]

Regarding the name FieldHandler. Maybe DofRegion, CellGroup, Region or similar could work? If we rename it SubDofHandler, it sounds like it has the same complexity and functionality as the DofHandler, like distributing dofs, keeping track of celldofs etc (which it cant).

I like the idea of pushing seperate fields to the "SubDofHandler" though.

[fredrikekre]

keeping track of celldofs etc (which it cant).

I was thinking it should do this though, and would enable cell iteration over the subdofhandler directly, with all the same queries as you can do on the standard dofhandler.

[termi-official]

Thanks for writing this up Fredrik (and also for the hard work from Kim and Elias!).

Several questions here.

1. For the user-facing portion of the API, independent of concrete naming, wouldn't it make sense to have a common function add!(::AbstractDofHandler, ::Symbol, ::Interpolation, ::SomeSubdomainType) to define a field on some subdomain of a dof handler?
2. Maybe DomainHandler and/or SubdomainHandler instead of SubDofHandler? This way we can also internally control the dispatches.
3. What is the advantage the SubDofHandler carrying multiple fields+the dof handler?
4. Was there something wrong with the old cell iterator design where we had the subdomain/fieldhandler/... as a second argument? I.e.

for c in CellIterator(dh, subdomainhandler1)
    # ...
end

To concretize a bit on the data structures, I was thinking about something along these lines

struct DomainHandler{IP <: Interpolation}
   name::Symbol
   ip::IP
end

struct SubDomainHandler{IP <: Interpolation}
   name::Symbol
   ip::IP
   subdomain # view on cell array or array of cell indices?
end

Note that DomainHandler is similar to the current data structure Field.

[lijas]

I was thinking it should do this though, and would enable cell iteration over the subdofhandler directly, with all the same queries as you can do on the standard dofhandler.

Hmm, I dont see very much benefit with this I think 🤔 . The way we do it now feels more clear/direct:

for c in CellIterator(dh, cellset)
    # use dofhandler directly
end

[fredrikekre]

1. Maybe, but then I think we still need to partitione the domain based on the input, instead of asking the user to tell it explicitly. E.g., imagine if you do

   add!(dh, :u, 3) # all cells
   add!(dh, :p, 1, somesubdomain)

   then I think we would need to split the :u part into subsubdomain and fulldomain - subdomain to simplify things. In addition, at least for the usecase of multiple element types, you would need different interpolations anyway. Migth not be worth having this shortcut for the case of uniform element type and different fields on different regions.

2. I wasn't thinking to make SubDofHandler <: AbstractDofHandler, but maybe it would simplify things. Otherwise I don't understand why the name matters for the dispatch.

3. I don't understand the question, advantage over what?

4. (+ Elias comment). That probably works just fine, but there is no guarantee what celltype is included in that cellset, and you would still need to send around that cellset instead of a subdofhandler. My thinking was that we have

   struct SubDofHandler{CT}
       dh::DofHandler
       cellset::Set{Int}
   end
   

   where we can compute CT at construction (since each subdomain is required to have a concrete type, already now). Then we now that when iterating over a SubDofHandler is must return elements of type CT.

[fredrikekre]

So I guess either you pass around a cellset, or you pass around a subdofhandler. For the latter you can also query ndofs_per_cell and stuff like that, which is pretty nice IMO.

[koehlerson]

4. (+ Elias comment). That probably works just fine, but there is no guarantee what celltype is included in that cellset, and you would still need to send around that cellset instead of a subdofhandler. My thinking was that we have
   ```
   struct SubDofHandler{CT}
       dh::DofHandler
       cellset::Set{Int}
   end
   ```
   where we can compute `CT` at construction (since each subdomain is required to have a concrete type, already now). Then we now that when iterating over a `SubDofHandler` is must return elements of type `CT`.

I'd vote for this as well

[termi-official]

Answering in order again.

1. For this specific use case you need to partition the domain in both approaches into subdomain1 and subdomain2. But maybe I am missing the point why one approach should be better (or worse) than the other here.
2. Sorry, I was not specific here. The name suggestion was more about clarifying the role of the type. SubDofHandler suggests that it might be an AbstractDofhandler subtype. I mean, I think we can explore different designs and see what the advantages/disadvantages are. Regarding finer control over the dispatches, we can have a path for the full domain DomainHandler which does not require a full copy of all cells and a path with finer SubDomainHandler (over having the SubDofHandler carrying e.g. a set of all indices - that is at least my current understanding of the proposed design). Maybe my question here is rather if you could clarify how a field defined on the full domain is intended to be handled?
   3&4. I will come back to this later. Will try to throw together some prototype of my idea to concretize.

[fredrikekre]

Maybe my question here is rather if you could clarify how a field defined on the full domain is intended to be handled?

Just like in the mixed dof handler now; you need to specify it for every subdomain you want it.

[fredrikekre]

Since we (I) want to store the parent dof handler we need something like:

dh = NewDofHandler(grid)
# Domain 1: u
sdh1 = SubDofHandler(dh, set1)
add!(sdh1, :u)
# Domain 2: u and p
sdh2 = SubDofHandler(dh, set2)
add!(sdh2, :u)
add!(sdh2, :p)
close!(dh)

i.e. pass the dofhandler to the subdofhandler constructor. At least you don't need the add!(dh, shd) calls.

[fredrikekre]

Actually, we could figure out the parent type also in the original suggestion, make it a mutable union with Nothing, and insert the parent after.

Another alternative spelling would be

dh = NewDofHandler(sdh1, sdh2, ..)

which could also perhaps imply close!.

[KnutAM]

Very nice to see this progress!
Regarding the "Implementation plan", I would prefer to make the change in one release at least to avoid having to update dependent libraries twice.

dh = NewDofHandler(sdh1, sdh2, ..)

This would be very nice and we could even allow dh = NewDofHandler(grid, fieldname, dim) as a one-liner setup for a single-field case:)
For setting up, we could also have a FieldDomain which is only used for setting up the NewDofHandler,

struct FieldDomain
 cellset::Set{Int}
 fields::Vector{Field}
end

where in Field the interpolation field is changed to Union{Interpolation,Nothing} (where nothing implies using the default interpolation). Then, upon calling dh = NewDofHandler(grid, fd1, fd2, ...) to create and close the NewDofHandler, which contains SubDofHandlers that are created based on the FieldDomains.
This way, the NewDofHandler and the SubDofHandlers can be better typed, e.g. the number of SubDofHandlers (same as above), but also the specific interpolations for each SubDofHandler.

[termi-official]

Implementing the one liner is already possible with something along these lines:

function DofHandler(grid, fieldname, dim)
    dh = DofHandler(grid)
    ip = default_interpolation(getcelltype(dh.grid))
    add!(dh, fieldname, dim, ip)
    close!(dh)
    return dh
end

where we should be able to do something similar for the MixedDofHandler by iterating over the cell types in the grid and adding the default interpolations individually.

@fredrikekre Assuming your SubDofHandler design. Let us take for example a grid with 3 elements, where the we add the field :u1 to element 1, then :u2 to elements 1,2,3 and finally :u3 to element 3. Any idea how to make a functional dof_range that is not too expensive? This case should show that, assuming that we allow overlapping subdomains as suggested above, then in the worst case dof_range needs to be computed (or cached) per element.

[KnutAM]

Implementing the one liner is already possible with something along these lines:

Yes, even simpler: DofHandler(grid::Grid, fieldname::Symbol, dim::Int) = (dh=DofHandler(grid); add!(dh, fieldname, dim); close!(dh)) (But would be a different approach and I have found it a bit strange to have some constructors closing automatically and some that don't)

Edit 2nd comment: Misread

[termi-official]

Indeed. :)

[kimauth]

The subdofhandler design is in fact what the MixedDofHandler does already, with the FieldHandlers representing the subdomains. An important point is that the cellsets must be disjoint, as documented. This is also checked during construction. This results in one dof_range per Fieldhandler for each field. You can cache them, they will never change after constructing the MixedDofHandler, but they're also not that expensive to compute since #621.
So your example is possible already @termi-official, it would look similar to what the first post shows. (And a point here is to find nicer syntax for it.)

[termi-official]

I see. Just as you wrote I also discovered the related code (but really missed the docs). I totally misunderstood Fredrik's answer, but now it makes sense. Thanks for elaborating on this Kim! This also explains my third question above, why we should have multiple fields per FieldHandler (or whatever we want to call it) in this design.

[fredrikekre]

Looks like everyone is on track now then. I think the implementation plan above will be too confusing. Instead I suggest doing what @kimauth originally suggested: we adjust the MixedDofHandler to be the new thing. This can be done in individual PRs making it easy to review. Finally, we delete DofHandler, and rename MixedDofHandler.

[termi-official]

#625 aims to close the performance gap on several calls when having a FieldHandler on the full domain (i.e. what DofHandler functionality provides).

[termi-official]

#627 adds feature parity between the DofHandlers for high order dof distribution.