WeakDirichlet 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

Attention: Patch coverage is 95.12195% with 6 lines in your changes missing coverage. Please review.

Project coverage is 93.77%. Comparing base (5f20e7b) to head (fbc2a53).
Report is 4 commits behind head on master.

Files with missing lines	Patch %	Lines
src/Dofs/ConstraintHandler.jl	95.49%	5 Missing ⚠️
src/interpolations.jl	90.90%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1151      +/-   ##
==========================================
+ Coverage   93.69%   93.77%   +0.07%     
==========================================
  Files          39       39              
  Lines        6237     6361     +124     
==========================================
+ Hits         5844     5965     +121     
- Misses        393      396       +3     

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

[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.

[KnutAM]

Currently, we don't have any H(div) or H(curl) interpolations with different number of shape functions associated with different facets. But theoretically, we could, hence this addition which is not tested.

Removing it will case a singular matrix if such are implemented in the future, alternatively we could check and throw if size(K,1) != length(shape_nrs).

[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)