Autodiff on sub-domains

[ConnorMallon]

There are some cases when using Gridap(Embedded) and Gridap's Autodiff machinery together does not work. Here we have a problem where a FEFunction only exists in certain parts of the domain. It means that if we do something like get_cell_dof_values we get back an array for which some of the entries are empty. The issue arises then when we try to find the element type of this object. Since the current type inference occurs by querying the type of the first element (testitem) then we might get back the type representing an empty array in the case which the FEFunction does not exist in the first cell. I have tried to set some of the array types manually but I have gotten stuck (e.g. gridap/Gridap.jl#653). I dont think what I have done so far is the best way to go anyway. I am guessing that what is needed to be done here is to treat the indicies which are non-empty separately from those that are empty by exploiting the tools in Arrays/PosNegReindex.jl and then having a special map for dealing with these dual number arrays which have some empty entries. The following MWE would then hopefully work. I think this is probably more of an issue for Gridap not GridapEmbedded but I put it here because the test is only a slight change from GridapEmbeddedTests/BimaterialPoissonCutFEMTests.jl

module GridapEmbeddedAutodiffTests

using Gridap
import Gridap: ∇
using GridapEmbedded
using Test
using Gridap.FESpaces
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.MultiField
using GridapEmbedded
using GridapEmbedded.LevelSetCutters

# Select geometry
const R = 0.7
geo1 = disk(R)
geo2 = !geo1
n = 10
domain = (-1,1,-1,1)
partition = (n,n)

# Setup background model
bgmodel = simplexify(CartesianDiscreteModel(domain,partition))

# Cut the background model
cutgeo = cut(bgmodel,union(geo1,geo2))

# Setup models
model1 = DiscreteModel(cutgeo,geo1)
model2 = DiscreteModel(cutgeo,geo2)

# Setup integration meshes
Ω1 = Triangulation(cutgeo,geo1)
Ω2 = Triangulation(cutgeo,geo2)

# Setup Lebesgue measures
order = 1
degree = 2*order
dΩ1 = Measure(Ω1,degree)
dΩ2 = Measure(Ω2,degree)

# Setup FESpace
V2 = TestFESpace(model1,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
V1 = TestFESpace(model2,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
U1 = TrialFESpace(V1)
U2 = TrialFESpace(V2)
V = MultiFieldFESpace([V1,V2])
U = MultiFieldFESpace([U1,U2])

uh = FEFunction(V1,ones(num_free_dofs(V1)))
j(u)=∫(u)dΩ1
∇(j)(uh)

uh = FEFunction(V2,ones(num_free_dofs(V2)))
j(u)=∫(u)dΩ2
∇(j)(uh)

j((u1,u2)) = ∫(u1)dΩ1 + ∫(u2)dΩ2
js=∇(j)(uh)
assemble_vector(js,V)

end # module