interpolation onto Discontinuous Lagrange Trace silently gives wrong answer

[Marie Elisabeth Rognes]

Consider the below code snippet interpolating the constant 1.0 onto a Discontinuous Lagrange Trace. I would expect resulting printed array to be [1, 1, 1, 1, 1], and its norm to be sqrt(5), but it gives a bunch of zeros. Two problems with this:

1. Wrong answer computed silently
2. Lots of "*** FFC warning: evaluate_dof(s) for enriched element not implemented" warnings emitted, unclear what enriched elements have to do with this.

(Probably a problem further down the chain, but convenient to test with interpolate.)

from fenics import *
mesh = UnitSquareMesh(1, 1)
V = FiniteElement("Discontinuous Lagrange Trace", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, V)
v0 = Expression("1.0", degree=0)
v0 = interpolate(v0, V)
print v0.vector().array()
print "|v0| = ", v0.vector().norm("l2")

[Marie Elisabeth Rognes]

Note that this used to work fine at at least one changeset between 2016.1 and and 2016.2 :-) @miklos1 @blechta any ideas?

[Jan Blechta]

When FIAT does not provide nodal basis, evaluate_dof needed for interpolate is not generated. This was so far only the case for enriched, that's why misleading message.

FIAT.HDivTrace (the replacement of DiscontinuousLagrangeTrace) is not nodal according to its specifications. Dual basis orthogonal to primal basis is not available. Somebody who knows the mathematics here, can it be added?

I can fix the misleading warning message. But I don't regarding the element itself.

I also think that the warning message should be an error.

[Jan Blechta]

There might be more serious problem as v0 = project(v0, V) (which would normally be a working workaround for non-nodal element) fails when FFC generates some nans.

[Lawrence Mitchell]

If you project on the trace, it should be fine. But project will use a quadrature rule on the cell, and the basis is not defined in the cell, only on its trace.

[Marie Elisabeth Rognes]

The dual basis did exist for this element at some point in the past (as far as I remember and as far as results with old version shows).

I'm not surprised that default

project

gives nan's, that's fine.

[Miklós Homolya]

Does this patch help?

diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index 4e40019..2d78fb4 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -275,6 +275,10 @@ class HDivTrace(FiniteElement):
         """Return number of members of the expansion set."""
         raise NotImplementedError("get_num_members not implemented for the trace element.")

+    @staticmethod
+    def is_nodal():
+        return True
+

 def construct_dg_element(ref_el, degree):
     """Constructs a discontinuous galerkin element of a given degree

[Lawrence Mitchell]

WFM. [Edit: With Miklos' patch, I should say]

[Marie Elisabeth Rognes]

@miklos1 Yes, much better. Could you push?

[Jan Blechta]

I insist here being careful and push that fix only if there is a test for nodality. See <FIAT>/tets/unit/test_fiat.py.

[Jan Blechta]

Unless nobody is able to write the test, you can always do monkey-patching in your application code, @meg.

[Marie Elisabeth Rognes]

Was this ever fixed?

[Miklós Homolya]

@meg: I made a patch above, then Jan insisted on also adding a test for nodality, then nothing happened. So the fix is known, but it's not in master.

[Marie Elisabeth Rognes]

Ok! The DGT seems to be overall broken (again), so I'll revisit this shortly.

[Marie Elisabeth Rognes]

@blechta I'd like to get this fixed now cf. https://bitbucket.org/fenics-project/fiat/branch/meg/fix-dolfin-issue-871. I can add the DOLFIN test in the report or alternatively could you contribute with the test you request?

[Marie Elisabeth Rognes]

Fixed via FIAT.