Bug: Autodiff when variable is not defined on the whole integration domain.

[JordiManyer]

The following script fails:

using Gridap
using Gridap.CellData

model = CartesianDiscreteModel((0,1,0,1),(2,2))

Ω_space = Triangulation(model,[1,2])

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(Ω_space,reffe)

Ω = Triangulation(model)
dΩ = Measure(Ω,2)
j(u) = ∫(u)dΩ

uh = zero(V)
c = get_contribution(gradient(j,uh),Ω)

collect(c) # Fails

The issue happens when the variable being differentiated is not defined on the whole integration domain. Although a niche thing, this kind of scenario does happen quite often in embedded scenarios.

It is caused by how autodiff works, and how cell-wise arrays are reindexed to match eachother's length.
I will be working on a fix.

@zjwegert

[JordiManyer]

Possible fix:

function FESpaces._compute_cell_ids(uh,ttrian)
  strian = get_triangulation(uh)
  if strian === ttrian
    return collect(IdentityVector(Int32(num_cells(strian))))
  end
  @check is_change_possible(strian,ttrian)
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  @notimplementedif !isa(sglue,FaceToFaceGlue)
  @notimplementedif !isa(tglue,FaceToFaceGlue)
  scells = IdentityVector(Int32(num_cells(strian)))
  mcells = extend(scells,sglue.mface_to_tface)
  tcells = lazy_map(Reindex(mcells),tglue.tface_to_mface)
  # <-- Remove collect to keep PosNegReindex
  return tcells
end

# New dispatching
function Arrays.lazy_map(k::Reindex,ids::Arrays.LazyArray{<:Fill{<:PosNegReindex}})
  k_posneg = ids.maps.value
  posneg_partition = ids.args[1]
  pos_values = lazy_map(Reindex(k.values),k_posneg.values_pos)
  pos_values, neg_values = Geometry.pos_neg_data(pos_values,posneg_partition)
  lazy_map(PosNegReindex(pos_values,neg_values),posneg_partition)
end

function Arrays.evaluate!(result,k::AutoDiffMap,ydual,x,cfg::ForwardDiff.GradientConfig{T}) where T
  @notimplementedif ForwardDiff.chunksize(cfg) != length(x)
  @notimplementedif length(result) != length(x)
  !isempty(x) && ForwardDiff.extract_gradient!(T, result, ydual) # <-- Watch for empty cell contributions
  return result
end

But I am worried about performance since:

• We are not collecting the cell ids, which means they are going to be computed on demand multiple times
• We are adding some branching into a very low-level evaluate function, which is not ideal.

[zjwegert]

The collect is still broken if computed using the LazyArrays machinery:

function lazy_collect(a)
    c = array_cache(a)
    for i in eachindex(a)
        getindex!(c,a,i)
    end
end

Throws a type-instability error...