Gridap.gradient fails when dealing with two domains of different dimensionalities

[wei3li]

Here is an MWE:

using Gridap

reffe = ReferenceFE(lagrangian, Float64, 1)
model = CartesianDiscreteModel((0, 1, 0, 1), (5, 5))
diri_tags, neum_tags = [1, 2, 3, 4, 5, 6, 7], [8]
Ω = Triangulation(model)
Γ = BoundaryTriangulation(Ω, tags=neum_tags)
dΩ, dΓ = Measure(Ω, 3), Measure(Γ, 3)
V0 = TestFESpace(model, reffe; dirichlet_tags=diri_tags)
Qh = FESpace(Γ, reffe)
l(γ, v) = ∫(-1.0 * v)dΩ + ∫(γ * v)dΓ
γh = FEFunction(Qh, rand(num_free_dofs(Qh)))
λh = FEFunction(V0, rand(num_free_dofs(V0)))
Gridap.gradient(γ -> l(γ, λh), γh)

The error traces are:

ERROR: AssertionError: A check failed
Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/Gridap/971dU/src/Helpers/Macros.jl:60 [inlined]
 [2] _compute_cell_ids(uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, ttrian::Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}})
   @ Gridap.FESpaces ~/.julia/packages/Gridap/971dU/src/FESpaces/FEAutodiff.jl:34
 [3] _gradient(f::Function, uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, fuh::Gridap.CellData.DomainContribution)
   @ Gridap.FESpaces ~/.julia/packages/Gridap/971dU/src/FESpaces/FEAutodiff.jl:22
 [4] gradient(f::var"#7#8", uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
   @ Gridap.FESpaces ~/.julia/packages/Gridap/971dU/src/FESpaces/FEAutodiff.jl:5
 [5] top-level scope
   @ ~/Documents/code/adjoint_method.jl:209

In the MWE, if l(γ, v) = ∫(γ * v)dΓ, no errors will occur.

Explanation of the error from @amartinhuertas:

The cause of the problem is that Gridap.gradient under the hood checks whether it is possible the change the domain of γh (i.e., Γ) to the domain associated to the measure dΩ (i.e., Ω), no matter whether γh is part of the integrand of the term or not. Thus, is_change_possible returns false, because γh is defined on a lower-dimensional manifold (i.e., 1D versus 2D).
This in particular happens here

function _gradient(f,uh,fuh::DomainContribution)
  terms = DomainContribution()
  for trian in get_domains(fuh)
    g = _change_argument(gradient,f,trian,uh)
    cell_u = get_cell_dof_values(uh)
    cell_id = _compute_cell_ids(uh,trian)
    cell_grad = autodiff_array_gradient(g,cell_u,cell_id)
    add_contribution!(terms,trian,cell_grad)
  end
  terms
end

In the call to _compute_cell_ids(uh,trian), uh is γh and trian is Ω in the particular context of your MWE.
In other words, it is assuming somehow that all the terms in the functional being differentiated depend on the input argument
I dont think it is correct, the first term in your functional does not depend on the input argument, therefore Gridap should not assume that
I think it is not possible to determine whether a given element in fuh depends on uh given the information available in the local scope of _gradient