Interface triangulations and CellFields

# some more code ...
# The heart of the problem: mean and jump operators
@inline _jump(w1, w2) = w1.⁺ - w2.⁻
@inline _mean(w1, w2) = w1.⁺ + w2.⁻

# this fails - namely, the calculated concentration field is discontinous
# even though this is used in terms used to enforce continuity
@inline _mean(u1, u2, n) = (D_field.⁺ ⋅ ∇(u1.⁺) ⋅ n.⁺ + D_field.⁻ ⋅ ∇(u2.⁻) ⋅ n.⁻)
# strangely enough, this works as expected:
@inline _mean(u1, u2, n) = (D1 ⋅ ∇(u1.⁺) ⋅ n.⁺ + D2 ⋅ ∇(u2.⁻) ⋅ n.⁻)
```

As implied above, if I define the mean operator on the interface as in the first case, I obtain discontinuities in the concentration field, which is not the case using the second approach. For clarity, I now also append the preparation of FESpaces, weak forms and relevant operators.

```julia

# preparation of FESpaces
order = 2
degree = 2*order

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model, tags="Casing")
dΓ = Measure(Γ, degree)
n_Γ = get_normal_vector(Γ)
dΩ = Measure(Ω, degree)

Ω1 = Interior(model, tags="face1")
Ω2 = Interior(model, tags="face2")
dΩ1 = Measure(Ω1, degree)
dΩ2 = Measure(Ω2, degree)

Γ_if = InterfaceTriangulation(Ω1, Ω2)
d_Γ_ifs = Measure(Γ_if, degree)
n_Γ_ifs = get_normal_vector(Γ_if)

reffe = ReferenceFE(lagrangian, Float64, order)
V1 = TestFESpace(Ω1, reffe, conformity=:H1)
V2 = TestFESpace(Ω2, reffe, conformity=:H1)
U1 = TransientTrialFESpace(V1)
U2 = TransientTrialFESpace(V2)

V = MultiFieldFESpace([V1, V2])
U = TransientMultiFieldFESpace([U1, U2])

D1 = 1.
D2 = 1.

C1 = 1.
C2 = 1.

flux(t) = t <= 1 ? sin(π * t) : 0.

# Diffusion constants

tags_field = CellField(domain_tags, Ω)
D_from_tag(tag) = tag == 1 ? D1 : D2 # very simple mapping function in this case
D_field = D_from_tag ∘ tags_field

# weak form

function _prefactors(triangulation, measure, param)

dim = num_cell_dims(triangulation)
_hmat = get_array(∫(1)measure)
_α = CellField(lazy_map(h -> param/h.^dim, _hmat), triangulation)

return _α
end
_α = _prefactors(Γ_if, d_Γ_ifs, 0.1)

# Note that, since we no longer have a single global FE space, we need to provide
# left and right functions to define jumps and averages at the interface.
@inline _jump(w1, w2) = w1.⁺ - w2.⁻
@inline _mean(w1, w2) = w1.⁺ + w2.⁻

# this works
@inline _mean(u1, u2, n) = (D1 ⋅ ∇(u1.⁺) ⋅ n.⁺ + D2 ⋅ ∇(u2.⁻) ⋅ n.⁻)
# this doesn't - a discontinuity in concentration is produced;
# @inline _mean(u1, u2, n) = (D_field.⁺ ⋅ ∇(u1.⁺) ⋅ n.⁺ + D_field.⁻ ⋅ ∇(u2.⁻) ⋅ n.⁻)
# note that in regular volume and surface integrations, I can work with CellFields just fine.

# for some reason, this works better with no symmetrization term in the Nitsche scheme
res(t, (u1, u2), (v1, v2)) = begin

∫(∂t(u1) ⋅ v1)dΩ1 + ∫(∂t(u2) ⋅ v2)dΩ2 +
∫(-flux(t) ⋅ v1)dΓ + # only on one side this time
∫((D_field ⋅ ∇(u1)) ⋅ ∇(v1))dΩ1 + ∫((D_field ⋅ ∇(u2)) ⋅ ∇(v2))dΩ2 -
∫(_jump(v1, v2) ⋅ _mean(u1, u2, n_Γ_ifs))d_Γ_ifs +
∫(_α ⋅ _jump(u1, u2) ⋅ _jump(v1, v2))d_Γ_ifs

end

jac(t, (u1, u2), (du1, du2), (v1, v2)) = begin

∫((D_field ⋅ ∇(du1)) ⋅ ∇(v1))dΩ1 + ∫((D_field ⋅ ∇(du2)) ⋅ ∇(v2))dΩ2 -
∫(_jump(v1, v2) ⋅ _mean(du1, du2, n_Γ_ifs))d_Γ_ifs +
∫(_α ⋅ _jump(du1, du2) ⋅ _jump(v1, v2))d_Γ_ifs

end

jac_t(t, (u1, u2), (dtu1, dtu2), (v1, v2)) = begin

∫(dtu1 ⋅ v1)dΩ1 + ∫(dtu2 ⋅ v2)dΩ2

end

# Operators

init_guess = interpolate_everywhere([c1, c2], U(0.))
op = TransientFEOperator(res, jac, jac_t, U, V)

# As said above, I then use some convenience routines to solve this via DifferentialEquations.jl routines

```

Can someone reproduce this behavior or would be willing to look into it? I would be happy to cooperate on this issue. Interestingly, this same issue does not appear if I use global FESpaces U and V (hence if I do not break the problem into subdomains), there the CellFields work nicely on the interfaces. The reason I am bringing this up is that I would eventually want to treat more complex diffusion cases with possible discontinuities in certain quantities (think of generalized diffusion with continous chemical potentials across interfaces and consequently discontinous concentration fields).

Also, the issue with interfaces and AD still persists (it probably boils down to some missing signature) - namely, one needs to manually add the Jacobians if working with interfaces in transient problems, as otherwise AD will complain about some function belonging to an interface definition.