ProjectedDirichlet BC for H(div) and H(curl)

[KnutAM]

For H(div) and H(curl) interpolations, the dof is not associated with a point, but with facets and edges, respectively. However, the boundary conditions appear for the facets in the forms described in the docstring. This implies that we cannot get a coordinate where we want to assign the dof, but instead we have to minimize the difference, not sure how much theory to include, but we could add for e.g. H(div) that our "problem" is,

$$\begin{align*} &\boldsymbol{q}(\boldsymbol{x}) \cdot \boldsymbol{n} = f \quad \text{ on } \Gamma, \quad \text{Not possible for arbitrary } f(\boldsymbol{x},...) \\\ &\text{Instead}\\\ &\boldsymbol{q}(\boldsymbol{x}) = \mathrm{argmin} U(\boldsymbol{q}(\boldsymbol{x})), \quad U(\boldsymbol{q}(\boldsymbol{x})) = \int_\Gamma \left[\boldsymbol{q}(\boldsymbol{x}) \cdot \boldsymbol{n} - f\right]^2\ \mathrm{d}\Gamma \\\ &\text{Found as solution to the directional derivative in all possible directions,}\delta\boldsymbol{q}(\boldsymbol{x})\text{. I.e., for }\epsilon \rightarrow 0, \\\ &0 = \frac{\partial}{\partial \epsilon} U(\boldsymbol{q}(\boldsymbol{x}) + \epsilon \delta\boldsymbol{q}(\boldsymbol{x}) = \int_\Gamma \left[\boldsymbol{q}(\boldsymbol{x}) \cdot \boldsymbol{n} - f\right]\delta\boldsymbol{q}(\boldsymbol{x}) \cdot \boldsymbol{n}\ \mathrm{d}\Gamma \\ &\text{Introducing the FE approximation for test and trial, we then get what we have in the docstring},\\\ &\underbrace{\int_{\Gamma^f} \left[\delta\boldsymbol{N}^f_i \cdot \boldsymbol{n}^f\right]\left[\boldsymbol{N}^f_j \cdot \boldsymbol{n}^f\right]\ \mathrm{d}\Gamma}_{K^f_{ij}}\ a_j^f = \underbrace{\int_{\Gamma^f} \left[\delta\boldsymbol{N}^f_i \cdot \boldsymbol{n}^f\right] q_n\ \mathrm{d}\Gamma}_{f^f_i} \end{align*}$$

where the functional dependence, $f = f(\boldsymbol{x}, t, \boldsymbol{n})$, is skipped fore brevity.

I also chose to introduce the internal function_space(::Interpolation)::Val function, as this was useful for dispatching. @termi-official: RannacherTurek and CrouzeixRaviart should be L2, right?

Side-note

An alternative name could be L2Dirichlet, as this could also be applied to $H^1$ spaces as well, and then it would actually be more accurate in the L2 sense. The problem with implementing that is that we get contributions from neighboring facets, and thus it either becomes a challenging problem to implement, or we need to solve a boundary FE problem on the global scale (probably easier and not so bad considering the sparsity) in order to apply the boundary conditions. But in most (probably all) practical cases, the additional "accuracy" in applying the BC is offset by the fact that an FE solution wouldn't be accurate enough if resolving the actual shape of the boundary prescribing function is important. However, if using $H^1$ spaces to describe e.g. fluxes, then using a standard Dirichlet condition would give the wrong total flux value if one uses e.g. linear shape functions but prescribe a quadratic velocity profile (e.g. in fluid problems...).

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 94.02%. Comparing base (0b130ef) to head (3ece145).
Report is 1 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1151      +/-   ##
==========================================
+ Coverage   93.89%   94.02%   +0.13%     
==========================================
  Files          39       39              
  Lines        6436     6560     +124     
==========================================
+ Hits         6043     6168     +125     
+ Misses        393      392       -1     

☔ 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]

I also chose to introduce the internal function_space(::Interpolation)::Val function, as this was useful for dispatching.

Yea, as already discussed before, I think having a trait system will be useful in the future.

@termi-official: RannacherTurek and CrouzeixRaviart should be L2, right?

Yes.

[fredrikekre]

I haven't gone through the math but if it is sound and tests cover it then 👍

[termi-official]

Since I just saw it, should we rename FunctionSpace trait? Although technically correct, it might a bit misleading/ambiguous here. I think we rather want to describe something like conformity or continuity trait, right?

[KnutAM]

Since I just saw it, should we rename FunctionSpace trait? Although technically correct, it might a bit misleading/ambiguous here. I think we rather want to describe something like conformity or continuity trait, right?

Yes, seems like dealii agrees with you: https://dealii.org/current/doxygen/deal.II/classFiniteElementData.html#a0cd5f34c3ab828fac31004f3b52921a1 :)

[KnutAM]

I haven't gone through the math but if it is sound and tests cover it then 👍

I believe it is sound, but would be good if someone can have a 2nd look on the tests here to make sure that I haven't done anything wrong there.

When going through this again now, I realize that I should also test some moment-integrals, as currently only the total integral quantity is tested (e.g. the total flux over the boundary, but not the spatial distribution). Will update with that during the week.

[fredrikekre]

Is it too magical if Dirichlet(::HcurlInterpolation) adds a WeakDirichlet? Probably, but would be kind of nice to not have to change the BC when changing interpolations.

[KnutAM]

Is it too magical if Dirichlet(::HcurlInterpolation) adds a WeakDirichlet? Probably, but would be kind of nice to not have to change the BC when changing interpolations.

1. We might want to add WeakDirichlet to :H1 interpolations in the future (see "side note" in the first comment).
2. We might want to use Dirichlet for the method they have in dealii (@termi-official has some contacts that could help us figure out their exact implementation)

[termi-official]

Really cool work here, thanks!

I was asking myself if this infrastructure can help us with applying normal and tangential boundary conditions (e.g. static sliding contact)?

[KnutAM]

Thanks to @termi-official on Slack - it is not a local projection problem for 3D and H(curl), so will add an error message that this is not implemented.
There is nothing preventing it being implemented, but that gives the same requirements as for H1 that the equation system to be solved is not local.

[KnutAM]

@lijas identified the corresponding functionality in deal.ii!

For nodal and hdiv: https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a747e71d426643718d52eda904d46bb20

And for hcurl: https://www.dealii.org/current/doxygen/deal.II/group__constraints.html#ga190b302dc1ee7c89727c94c083d132d2

I will look closer if I'm missing something since they have factored out the curl-functionality, but I haven't done that here.

[KnutAM]

After going through the dealii docs, I've added the relevant reference here: https://ferrite-fem.github.io/Ferrite.jl/previews/PR1151/topics/boundary_conditions/#ProjectedDirichlet.

For the H(curl), AFAIU this is related to @termi-official's comment regarding 3d, and they use a stepwise algorithm to reduce the equation system size when it couples (first solving dofs associated with edges, and then those associated with faces). This doesn't affect what is implemented in this PR, nor the user interface.

[KnutAM]

Should be good from my side.
Do you have any further comments/suggestions @termi-official or @fredrikekre?

[termi-official]

LGTM

And for hcurl: https://www.dealii.org/current/doxygen/deal.II/group__constraints.html#ga190b302dc1ee7c89727c94c083d132d2

This is quite a smart strategy here. However, we need to implement EdgeValues and the related topological queries to reproduce this efficiently.

[KnutAM]

Thanks for the comments and input @termi-official, @fredrikekre, and @lijas!