New DL Trace fails with "Could not find unique facet to tabulate on"

[Marie Elisabeth Rognes]

The below code fails with

File "/home/meg/local/fenics/master/lib/python2.7/site-packages/FIAT/hdiv_trace.py", line 142, in tabulate raise TraceError("Could not find a unique facet to tabulate on.") TraceError: Could not find a unique facet to tabulate on. .

from dolfin import *

n = 16
mesh = UnitSquareMesh(n, n)

M = VectorElement("DG", mesh.ufl_cell(), 1)
V = FiniteElement("DG", mesh.ufl_cell(), 0)
Q = FiniteElement("Discontinuous Lagrange Trace", mesh.ufl_cell(), 1)

W = FunctionSpace(mesh, MixedElement([M, V, Q]))
(tau, v, q) = TestFunctions(W)
(sigma, u, p) = TrialFunctions(W)

n = FacetNormal(mesh)
a0 = inner(sigma, tau)*dx + u*div(tau)*dx + div(sigma)*v*dx

A = assemble(a0)

I do believe that it at least assembled with the old trace elements, could you take a look @TGibson1120?

[Thomas Gibson]

Sure thing, I'll have a look at this

[Thomas Gibson]

Hi Marie, Sorry for the delay.

So it seems the TraceError is working as it should, the problem is that the logic is a bit strict in the backwards compatibility component of hdiv_trace.py. I presume when you assemble a cell integral on a mixed element, tabulate is called for all component spaces. The problem here is that your quadrature rule is of the cell-type. None of the quadrature points will be on facet entities, hence the TraceError being raised.

If I recall correctly, before the changes the tabulation dict is initialized to have all zeros in its entries. In the case that a facet is not found, nothing happens and the zeros are silently propagated through everything else. At the end, the "correct" numerical result is obtained. However, now I'm being quite strict in asserting that the points are in fact lying on cell-facets.

How does this sound: instead of raising the TraceError, I can instead return a dict of NaNs. This would be mathematically sound as trace functions are not defined on the cell interiors. This should work I think, as the trace contributions should not end up in the generated code to begin with. Does this sound like a better approach? I will begin writing a patch for this.

[Thomas Gibson]

@meg I reproduced your error. And I went ahead and wrote a patch in FIAT doing exactly as I described above. It seems to work just fine now (no trace errors halting anything and FFC compiles the form a0). Are you happy with the approach? If so, I will file a FIAT PR asap.

[Miklós Homolya]

@meg: Is it now resolved?

[Marie Elisabeth Rognes]

Yeth.