Implement computation of dual basis for enriched element

[Jan Blechta]

By pull request #19 incorrect implementation of dual basis of enriched element is disabled (such that evaluate_dof[s] and tabulate_coordinates throw error). Proper dual basis \psi_j has to satisfy

\psi_j (\phi^i) = \delta_j^i

where \phi^i are basis functions from the enriched space. Such a dual basis can be used to calculate DOFs f_i of function f = f_i \phi^i (summation implied) by formula f_i = \psi_i (f) which is useful in DOLFIN for restrict_as_ufc_function.

Computing \psi_j requires a solve which could be done in ufc::element constructor or during code generation (or earlier in FIAT). The tricky thing is that user may request enriched element which is not unisolvent and solve could fail then. For example (rather trivial)

B = FiniteElement("Bubble", mesh, 3)
element = B + B # !!!

(Are there non-trivial examples and can we catch them?)

Proposed changes

1. Implement solver for searching for \psi_j^{new} = T_{jk} \psi_k s.t. \psi_j^{new} (\phi^i) = \delta_j^i into ir generation or ufc::element code (with "solve on first use" idiom). Handling unsolvable cases (related to element unisolvency) needs more thought.

2. Abstraction for representing \psi_j^{new} being linear combination of Dirac deltas. FIAT.functional.Functional seems appropriate and ready for this.

3. Rewrite ir and code generation with support for linear combination of Dirac deltas. Note: tabulate_coordinates will probably stay undefined for such a functionals. Or does it have a sense?

[Jan Blechta]

Fixed typo in LaTeX.

[Jan Blechta]

Add summary of proposed changes.

[Jan Blechta]

Reminder: don't forget to revert DOLFIN workaround 231662e when the proper implementation is in place.

[Jan Blechta]

I got swayed to recompute rather primal basis while keeping dual basis intact. This is in fact done in FIAT for every element. Thus elements constructed elsewhere (in FFC) may suffer from a lack of the duality relation. This applies to enriched element and maybe also symmetric tensor elements (see DOLFIN issue 516).

I have not yet decided how to fix it. Options are following:

1. Put all elements construction to FIAT and do it consistently, with a proper computation of a basis.

2. Don't extend FIAT. Put more of complicated logic to FFC.

[Jan Blechta]

Enriched element should be moved to FIAT, possibly by Firedrake merge.

[Prof Garth Wells]

Removing milestone: 1.8 (automated comment)

[Corrado Maurini]

I found this regression very annoying. EnrichedElements as Bubble elements are an important feature of FEniCS in my view. Is there any short term plan (the issue dates back 1.5 years) to fix it?

If a proper fix is complex, as it seems form the comment above, can any workaround be implemented meanwhile? Do you have any suggestions/instructions to continue to use EnrichedElements?

[Jan Blechta]

this regression

It is not a regression but rather a fix. Nodal interpolation has never been working for enriched and emitted code was rubbish.

Is there any short term plan

Yes. The planned fix includes moving the implementation from FFC to FIAT. It is quite easy but it is taking such a long time because the plan was in conflict with the Firedrake additions. FIAT part is so far in https://bitbucket.org/fenics-project/fiat/branch/jan/nodal-enriched-element-2 but I still need to find a way how to overcome how tensor product element handles dual basis. (Which is that it does not bother with dual basis at all as Firedrake apparently does not use dual basis.)

suggestions/instructions to continue to use EnrichedElements

Yes. Just avoid doing operations which need nodal interpolation. In DOLFIN it is:

• interpolation to enriched
• projection from enriched to another mesh (generally using enriched coefficient on non-matching mesh in a form)
• using Expression or Constant values in Dirichlet BCs on enriched space

On the other hand following operations should be working

• projection to enriched
• interpolation from enriched to non-enriched on matching mesh
• using enriched Function value in Dirichlet BC (with matching mesh)

[Corrado Maurini]

Thanks a lot @blechta. The projection on BC seems to do the trick and I can have the code with bubbles working again.

Corrado

[Jan Blechta]

NodalEnrichedElement implemented in 529b807 and ef4732b. Check demo/NodalMini.ufl.