Refactoring dof distrubtion to allow using generic local celldof orderings

[termi-official]

This PR essentially introduces the following:

1. I defined the new dof interface in Improve dof distribution interface #627 , but it was not complete, as the local ordering was fixed. This PR now allows arbitrary local permutations for the dofs. This simplifies implementing high order interpolations in two ways. A) The reasoning about facet permutations is more sane in 3D and 4D. B) We can easily setup tensor-product elements (e.g. Generic arbitrary order interpolations #626 ).
2. An algorithm parameter, so we can in the future allow users to use different dof distribution algorithms, e.g. to allow dof numbering follow the node numbering or to unlock p-refinement. This is also useful for FerriteDistributed if we modularize the algorithms further, so I can significantly reduce code duplication.

TODO

• Split off dof distribution algorithm interface
• test of permuting the indices for an existing interpolation, and validating that all Ferrite-tests pass for those
• remove permdof2DLagrange2Tri345 This is a task for another day.

[codecov]

Codecov Report

❌ Patch coverage is 92.50000% with 6 lines in your changes missing coverage. Please review.
✅ Project coverage is 93.08%. Comparing base (f99e206) to head (ed1999a).
⚠️ Report is 27 commits behind head on master.

Files with missing lines	Patch %	Lines
src/Dofs/DofHandler.jl	91.04%	6 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1188      +/-   ##
==========================================
- Coverage   94.19%   93.08%   -1.11%     
==========================================
  Files          40       40              
  Lines        6662     6764     +102     
==========================================
+ Hits         6275     6296      +21     
- Misses        387      468      +81     

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

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[termi-official]

May supersede #829 .

[KnutAM]

Cool!

introduces two things.

Is there are reason for discussing these in the same PR? Seems orthogonal to me?

1

This PR now allows arbitrary local permutations for the dofs

How does this affect the periodicity logic and boundary condition application (I believe the fixed ordering is assumed there). And what are the exact requirements to ensure that cells become compatible then (i.e. can I have different permutations for different faces on a cell, does the distribution correct for this?)

The reasoning about facet permutations is more sane in 3D and 4D

Could you elaborate on how this generalization makes facet permutation reasoning more sane?

We can easily setup tensor-product elements

Is that really not possible with the current interface, just providing a wrapper to give the correct numbering? If one could keep the current fixed numbering interface, I think that makes it easier to implement new features.
Then the tensor-product interpolations can remain an advanced feature that we don't need to account for in the dof distribution.

2

An algorithm parameter, so we can in the future allow users to use different dof distribution algorithms, e.g. to allow dof numbering follow the node numbering or to unlock p-refinement

Doesn't it make more sense to add such an interface once we actually implement a new distribution algorithm?

Specifically regarding nodal distriubtion: I think (but cannot find) this was discussed earlier, and the conclusion was that adding a nodal distribution to improve the dof-distribution performance is only negligible (compared to solving the system), so it is not worth having more code for this case, and that we should rather provide the tools necessary to work with non-nodal-based dof-distributions that work in general.

[termi-official]

introduces two things.

Is there are reason for discussing these in the same PR? Seems orthogonal to me?

1

I put them together to avoid the discussion on "what do we need this for".

This PR now allows arbitrary local permutations for the dofs

How does this affect the periodicity logic and boundary condition application (I believe the fixed ordering is assumed there). And what are the exact requirements to ensure that cells become compatible then (i.e. can I have different permutations for different faces on a cell, does the distribution correct for this?)

I think periodicity is broken with this

Ferrite.jl/src/Dofs/ConstraintHandler.jl

Line 340 in 0b130ef

push!(local_facet_dofs, (fdof - 1) * field_dim + d + offset)

(actually looks fine upon closer inspection). Dirichlet should not be affected, as it already uses the dof tables internally.

We can easily setup tensor-product elements

Is that really not possible with the current interface, just providing a wrapper to give the correct numbering? If one could keep the current fixed numbering interface, I think that makes it easier to implement new features. Then the tensor-product interpolations can remain an advanced feature that we don't need to account for in the dof distribution.

2

The stuff which I implement here has to be done in one way or another anyway. I would argue that it is easiest as in this PR, where we just use the tables to access dof information and distribute the dofs directly w.r.t. these tables.

Specifically regarding nodal distriubtion: I think (but cannot find) this was discussed earlier, and the conclusion was that adding a nodal distribution to improve the dof-distribution performance is only negligible (compared to solving the system), so it is not worth having more code for this case, and that we should rather provide the tools necessary to work with non-nodal-based dof-distributions that work in general.

Why should we not allow for flexibility here? I think making the dof distribution nodal by default will save us some ime in slack and it also makes things like VTK import more easy.

[termi-official]

I can also open a PR in FerriteInterfaceElements.jl to show how this relaxed assumption can lead to simpler interpolations.

[KnutAM]

(actually looks fine upon closer inspection). Dirichlet should not be affected, as it already uses the dof tables internally.

Nice! Perhaps related to

what are the exact requirements to ensure that cells become compatible then (i.e. can I have different permutations for different faces on a cell, does the distribution correct for this?)

The stuff which I implement here has to be done in one way or another anyway. I would argue that it is easiest as in this PR, where we just use the tables to access dof information and distribute the dofs directly w.r.t. these tables

My motivation was that not what was the easiest way, but what leads to the easiest code to work with for standard cases and for later implementation. I.e. if there are assumptions that simplifies our code now, which are no longer valid with this generalization. But probably related to above, that I don't fully understand what are now the allowed interface for different interpolations and how things are being ordered (I guess it is not fully flexible, i.e. still it needs to be numbered by the fields being consecutive, otherwise e.g. we don't have the local dof_range being a range?

Why should we not allow for flexibility here? I think making the dof distribution nodal by default will save us some ime in slack and it also makes things like VTK import more easy.

If some features are easy to implement for nodal distributions, do we want to have different code to support these different cases, or do we favour a single code that can handle it in general? I think the user-friendliness goes down a lot if we add features that only work if you use nodal dof-distribution?

[termi-official]

[...] if there are assumptions that simplifies our code now, which are no longer valid with this generalization [...]

The assumption starting with Ferrite 1.0 was that the dof ordering must follow the provided tables. However, we had at least two failed attempts to use the tables in close! without significantly degrading performance and hence we had to put the restriction on the tables, that the ordering has to be in increasing order.

I will add more test coverage for non-increasing dof tables if we reach agreement on the API in this PR, as I would like to move towards stabilizing the API, so users can implement their interpolations as they like (under the remaining invariants for the API).

I don't fully understand what are now the allowed interface for different interpolations

The only changing piece is that the tables for single interpolations can now be renumbered as you please, as intended in the original PR where the API was introduced. E.g. instead of numbering the dofs on your vertices of a quad with 1,2,3,4 for vertices 1,2,3,4 you can now number them e.g. 1,2,4,3 or 4,3,1,2 or any other permutation of 1,2,3,4, if you want.

how things are being ordered (I guess it is not fully flexible, i.e. still it needs to be numbered by the fields being consecutive, otherwise e.g. we don't have the local dof_range being a range?

Correct. Just how the individual interpolations themselves are distributed for a single section in the celldofs changes. There is no "cross-talk" allowed. The dof handler logic should capture if you go outside of your assigned range. The local dof ranges for the fields are the same as before, so dof_range is still returning the same information as before. Just the bijection between the local index into this resulting range and the location of the dof is allowed to change with this PR.

If some features are easy to implement for nodal distributions, do we want to have different code to support these different cases, or do we favour a single code that can handle it in general? I think the user-friendliness goes down a lot if we add features that only work if you use nodal dof-distribution?

Indeed, I agree here. From the experience I gather over the last years, most users seem to prefer the nodal dof distribution if possible. However, we can't make this default. I also do not want to make the distributed hp-adaptivity dof distribution algorithm the default either, as it is slower than what we have. I always prefer flexibility over being restricted for a slightly simper code.

For example, we could right now distribute the dofs field by field and this would simplify the dof distribution algorithm quite a bit. However, this would also degrade cache performance for multi-field problems, since the data is not localized anymore.

For loading VTK data I have a very simple interface á la

vtkdh = DofHandler(VTKFile("my-solution.vtu"))

in mind, where we might also add some options here. The functionality can be put into an extension.

[KnutAM]

Correct. Just how the individual interpolations themselves are distributed for a single section in the celldofs changes. There is no "cross-talk" allowed. The dof handler logic should capture if you go outside of your assigned range. The local dof ranges for the fields are the same as before, so dof_range is still returning the same information as before. Just the bijection between the local index into this resulting range and the location of the dof is allowed to change with this PR.

OK so I can have for a Triangle

vertexdof_interior_indicies = ((1,), (3,), (2,))
edgedof_interiori_indices = ((5,), (6,), (4,))

But not

vertexdof_interior_indicies = ((1,), (6,), (2,))
edgedof_interiori_indices = ((5,), (3,), (4,))

?

[termi-official]

Both fine. You are not allowed to do shenanigans like vertexdof_interior_indicies = ((1,), (7,), (2,)) or vertexdof_interior_indicies = ((-1,), (3,), (2,)) for your quadratic triangle (assuming you have exactly 6 dofs).

[KnutAM]

OK, great.

I think I'm understanding it more now. So the only requirements to be fullfilled are

• The set of dof indices for one interpolation must be equal to Set(1:getnbasefunctions(ip)).
• The dof-location (ref Make constraint handling tools available for users #1158) (i.e. where the base function is one for nodal interpolations) must be ordered in the positive direction along edges (and in the future when that has been defined, according to the positive face ordering)?

[termi-official]

OK, great.

I think I'm understanding it more now. So the only requirements to be fullfilled are

* The set of dof indices for one interpolation must be equal to `Set(1:getnbasefunctions(ip))`.

Yes.

* The dof-location (ref [Make constraint handling tools available for users #1158](https://github.com/Ferrite-FEM/Ferrite.jl/pull/1158)) (i.e. where the base function is one for nodal interpolations) must be ordered in the positive direction along edges (and in the future when that has been defined, according to the positive face ordering)?

Technically a no. We just have this assumption for edges right now to make our lives a bit easier. However, users need to manually override the permute_and_set function in this case right now for high order elements (i.e. more than 1 dof per edge). More than 1 dof per face does still not work with this PR.