Face integration with MultiFESpace

[bmxam]

A major bug happens when a face integral is computed without an explicit dependance on all the unknowns.

The following script produces an error, no matter the mesh dimension or the function space degree. The error also happens for only two FESpaces, but it's easier to debug with three.

using Bcube

degree = 0

mesh = line_mesh(3)

fs = FunctionSpace(:Lagrange, degree)
_U = TrialFESpace(fs, mesh)
_V = TestFESpace(_U)

U = MultiFESpace(_U, _U, _U)
V = MultiFESpace(_V, _V, _V)

dΓ = Measure(InteriorFaceDomain(mesh), 2 * degree + 1)

# NOTE THAT V3 IS NOT PRESENT
l((v1, v2, v3)) = ∫(side_n(v1) + side_p(v2))dΓ

@show assemble_linear(l, V)

The problems comes from this line :

Bcube.jl/src/assembler/assembler.jl

Line 392 in 9f09919

values_i = map(identity, map(identity, map(first, unwrapValues)))

We try to apply first to the every members of the Tuple unwrapValues. Because v3 is not explicitely present, this Tuple looks like this:

unwrapValues = (((1.0,), Bcube.LazyOperators.NullOperator()), (Bcube.LazyOperators.NullOperator(), (1.0,)), Bcube.LazyOperators.NullOperator())

so the quoted line tries to apply first on Bcube.LazyOperators.NullOperator(), leading to an error. If v3 is explicitely present in the integral, the Tuple looks something like this (the difference is at the end):

unwrapValues = (((1.0,), Bcube.LazyOperators.NullOperator()), (Bcube.LazyOperators.NullOperator(), (1.0,)), ((1.0,), Bcube.LazyOperators.NullOperator()))

and everything is fine.

Of course a functionnal & brutal & non satisfying solution is to define:

Base.first(a::Bcube.LazyOperators.NullOperator) = a
Base.last(a::Bcube.LazyOperators.NullOperator) = a

I don't really know the better place / way to fix this problem. @ghislainb , you will know better than me.

[ghislainb]

This is not really a bug. The starting assumption of the first implemented version is that the function l(v) to be assembled with assemble_linear must depend linearly and explicitly on all variables in v. If it does not work with ...+0*side_n(v3), then there's a bug.
What you propose is an improvement, but it raises the question of what to do with constant terms that don't depend on v1, v2, v3 at the time of assembly.

[bmxam]

For me it's definitively a bug (maybe we knew about it, but it's still a bug!). If you do the same with a volume integral there is no problem. As for constant terms we need (but this is not urgent) to display an understandable error message.