Gradient of Integration domain contribution involving SkeletonTriangulation

[kishore-nori]

Currently, I am unable to take the gradient of integration DomainContribution with respect to the FEFunction's degrees of freedom, when the integral involves SkeletonTriangulation. The following is the MWE:

using Gridap

model = CartesianDiscreteModel((0.,1.,0.,1.),(3,3))

Λ = SkeletonTriangulation(model)
dΛ = Measure(Λ,5)

reffe = ReferenceFE(lagrangian,Float64,2)
V = TestFESpace(model,reffe,conformity=:L2)

u(x) = 0.0
U = TrialFESpace(V,u)

uh = FEFunction(U,rand(num_free_dofs(U)))

f(uh) =  ∫(mean(uh))*dΛ  # ∫(jump(uh)*jump(uh))*dΛ

fuh = f(uh)
∑fuh = ∑(fuh)

Gridap.gradient(f,uh)

The error Gridap.gradient(f,uh) call produces is as follows:

ERROR: This function is not yet implemented
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] macro expansion
   @ ~/.julia/packages/Gridap/kce3i/src/Helpers/Macros.jl:21 [inlined]
 [3] _compute_cell_ids(uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, ttrian::SkeletonTriangulation{1, 2, BoundaryTriangulation{1, 2, Gridap.Geometry.BodyFittedTriangulation{1, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, Gridap.Geometry.GridView{1, 2, Gridap.Geometry.UnstructuredGrid{1, 2, Float64, Gridap.Geometry.NonOriented, Nothing}, Vector{Int64}}, Vector{Int64}}, Gridap.Geometry.FaceToCellGlue{Vector{Int64}, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Reindex{Vector{Int8}}, 1, Tuple{Base.OneTo{Int64}}}, Int8, 1, Tuple{Vector{Int64}}}, FillArrays.Fill{Int8, 1, Tuple{Base.OneTo{Int64}}}}}})
   @ Gridap.FESpaces ~/.julia/packages/Gridap/kce3i/src/FESpaces/FEAutodiff.jl:40
 [4] _gradient(f::Function, uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}}, fuh::Gridap.CellData.DomainContribution)
   @ Gridap.FESpaces ~/.julia/packages/Gridap/kce3i/src/FESpaces/FEAutodiff.jl:23
 [5] gradient(f::typeof(f), uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
   @ Gridap.FESpaces ~/.julia/packages/Gridap/kce3i/src/FESpaces/FEAutodiff.jl:5

[amartinhuertas]

Hi @kishore-nori ... perhaps we can remove the terms ∫(∇(uh)⋅∇(uh))*dΩ + ∫((uh-u)*(uh-u))*dΓ from the functional to have an even more reduced MWE. Am I correct?

[kishore-nori]

Sure Alberto, I ll make the changes and update now!

[amartinhuertas]

Hi @kishore-nori,

I have came up with a piece of code that I think partially addresses the issue. Can you check that it produces the correct result? I say partially above because there are still some issues (marked with comments below) which we still have to think how to address. The most concerning one is the one tagged with "MAJOR ISSUE".

Take a look and let me know if you have doubts, etc. We can talk tomorrow about this during the meeting, etc.

using Gridap
using Gridap.Arrays
using ForwardDiff
using FillArrays

# ISSUE: we have to overload jump for type T below as the default
#        behaviour for SkeletonPair assumes that jump has been called with
#        the normal component of a Cell Field, i.e., jump(uh⋅nΛ), and thus
#        is (wrongly for our scenario) defined as a.plus+a.minus
T=Gridap.CellData.SkeletonPair{<:Gridap.CellData.CellFieldAt}
Gridap.CellData.jump(a::T) = a.plus-a.minus

function Gridap.FESpaces._change_argument(op, f,
                                          trian::SkeletonTriangulation,
                                          uh::Gridap.FESpaces.SingleFieldFEFunction)
  U      = Gridap.FESpaces.get_fe_space(uh)
  strian = get_triangulation(U)
  @assert Gridap.CellData.is_change_possible(strian,trian) "ERROR"
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(trian,Val(D))
  function g(cell_u_plus,cell_u_minus)
    cf_plus  = CellField(U,cell_u_plus).plus
    cf_minus = CellField(U,cell_u_minus).minus
    # ISSUE: Hard-coded calls to change_domain_ref_ref,
    # To-think. when should we call change_domain_phys_phys?
    cf_plus=Gridap.CellData.change_domain_ref_ref(cf_plus,trian,sglue,tglue)
    cf_minus=Gridap.CellData.change_domain_ref_ref(cf_minus,trian,sglue,tglue)
    # MAJOR ISSUE: At present only jump() can be called with an SkeletonPair
    #        Any other term in f(cf) (e.g., integration on boundary or interior)
    #        will trigger an error. To-Think.
    cf=Gridap.Geometry.SkeletonPair(cf_plus,cf_minus)
    cell_grad_plus = f(cf)
    Gridap.CellData.get_contribution(cell_grad_plus,trian)
  end
  g
end

function Gridap.FESpaces._compute_cell_ids(uh,ttrian::SkeletonTriangulation)
  tcells_plus  = Gridap.FESpaces._compute_cell_ids(uh,ttrian.plus)
  tcells_minus = Gridap.FESpaces._compute_cell_ids(uh,ttrian.minus)
  Gridap.Geometry.SkeletonPair(tcells_plus,tcells_minus)
end

function Gridap.Arrays.autodiff_array_gradient(a,i_to_x,j_to_i::Gridap.Geometry.SkeletonPair)
  # Work for plus side
  i_to_xdual_plus = lazy_map(Gridap.Arrays.DualizeMap(ForwardDiff.gradient),i_to_x)
  j_to_ydual_plus = a(i_to_xdual_plus,i_to_x)
  j_to_x_plus = lazy_map(Reindex(i_to_x),j_to_i.plus)
  j_to_cfg_plus = lazy_map(Gridap.Arrays.ConfigMap(ForwardDiff.gradient),j_to_x_plus)
  j_to_result_plus = lazy_map(Gridap.Arrays.AutoDiffMap(ForwardDiff.gradient),
                              j_to_ydual_plus,j_to_x_plus,j_to_cfg_plus)

  # Work for minus side
  i_to_xdual_minus = lazy_map(Gridap.Arrays.DualizeMap(ForwardDiff.gradient),i_to_x)
  j_to_ydual_minus = a(i_to_x,i_to_xdual_minus)
  j_to_x_minus = lazy_map(Reindex(i_to_x),j_to_i.minus)
  j_to_cfg_minus = lazy_map(Gridap.Arrays.ConfigMap(ForwardDiff.gradient),j_to_x_minus)
  j_to_result_minus = lazy_map(Gridap.Arrays.AutoDiffMap(ForwardDiff.gradient),
                              j_to_ydual_minus,j_to_x_minus,j_to_cfg_minus)

  # Assemble on SkeletonTriangulation expects an array of interior of facets
  # where each entry is a 2-block BlockVector with the first block being the
  # contribution of the plus side and the second, the one of the minus side
  lazy_map(Gridap.Fields.BlockMap(2,[1,2]),j_to_result_plus,j_to_result_minus)
end

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

Λ = SkeletonTriangulation(model)
Γ = BoundaryTriangulation(model)


dΛ = Measure(Λ,0)
dΓ = Measure(Γ,0)

reffe = ReferenceFE(lagrangian,Float64,1)
V = TestFESpace(model,reffe,conformity=:L2)

u(x) = 0.0
U = TrialFESpace(V)

# uh = FEFunction(U,[1.0,1.0,1.0,1.0,5.0,5.0,5.0,5.0])
uh = FEFunction(U,rand(8))

f(uh) = ∫(jump(uh)*jump(uh))*dΛ

fuh = f(uh)
∑fuh = ∑(fuh)

grad_fuh=Gridap.gradient(f,uh)
assemble_vector(grad_fuh,U)

[amartinhuertas]

Can you check that it produces the correct result?

For the records, we confirmed that it produces the correct result.

[amartinhuertas]

To explore this alternative interface for the jumps, means etc. (TO-THINK)

f(uh)      = f(uh,uh)
f(uh+,uh-) = ∫(jump(uh+,uh-)*jump(uh+,uh-))*dΛ + ∫(uh+*vh)*dOmega

[amartinhuertas]

Another idea. TO-THINK. Rather than forcing the user to write his functionals differently (i.e., two input arguments instead of one), try to think about a new kind of CellField, say, SkeletonCellField, that holds inside two CellFieldAt{+/-} instances, and then to program the logic that let us support this kind of CellFields in the high level API.

[amartinhuertas]

Another idea. TO-THINK. Rather than forcing the user to write his functionals differently (i.e., two input arguments instead of one), try to think about a new kind of CellField, say, SkeletonCellField, that holds inside two CellFieldAt{+/-} instances, and then to program the logic that let us support this kind of CellFields in the high level API.

Something we have to think in regards to this. By design, CellFields are associated with a Triangulation and either defined on ReferenceDomain or PhysicalDomain. How all this fits into the concept of SkeletonCellField ?

[amartinhuertas]

We came up with the following draft design (to be completed) for SkeletonCellField. Let us see if it sufficient or not, there are no showstoppers, etc.

struct SkeletonCellField{P<:CellField,M<:CellField,T<:SkeletonTriangulation} <: CellField
  cf_plus::P
  cf_minus::M 
  trian::T
  function SkeletonCellField(cf_plus  :: CellField,
                             cf_minus :: CellField, 
                             trian::T) where T<:SkeletonTriangulation

    Gridap.Helpers.@check DomainStyle(cf_plus)==DomainStyle(cf_minus)
    Gridap.Helpers.@notimplementedif !(get_triangulation(cf_plus)===
                                        get_triangulation(cf_minus))

    ptrian=get_triangulation(cf_plus)

    Gridap.Helpers.@check num_dims(ptrian) == num_dims(trian)-1
    Gridap.Helpers.@check get_background_model(ptrian)===get_background_model(trian)

    # We ensure that we will store instances of CellFieldAt within SkeletonCellField
    cf_plus_plus  = cf_plus.plus
    cf_plus_minus = cf_minus.minus

    P = typeof(cf_plus_plus)
    M = typeof(cf_plus_minus)
    T = typeof(trian)
    new{P,M,T}(cf_plus_plus,cf_plus_minus,trian) 
  end 
end 

function Gridap.CellData.DomainStyle(a::SkeletonCellField)
  Gridap.Helpers.@check DomainStyle(a.cf_plus)
end 

function Gridap.CellData.get_triangulation(a::SkeletonCellField)
  a.trian
end 

function Gridap.CellData.change_domain(a::SkeletonCellField,
                       a_ds::DomainStyle,
                       ttrian::Triangulation,
                       o_ds::DomainStyle)
  change_domain(a.cf_plus,a_ds,ttrian,o_ds)
end 

function Gridap.CellData.change_domain(a::SkeletonCellField,
                       a_ds::DomainStyle,
                       ttrian::SkeletonTriangulation,
                       o_ds::DomainStyle)
  # Gridap.Helpers.@check false "Fatal Error: you cannot have an SkeletonCellField as an integrand of an integral over the skeleton"
  Gridap.Helpers.@notimplementedif !(get_triangulation(a)===get_triangulation(ttrian))
  a  
end 

uh->SkeletonTriangulation 
nΛ->SkeletonTriangulation 
uh⋅nΛ -> SkeletonPair !!!!



function Base.getproperty(a::SkeletonCellField,sym::Symbol)
  if sym==:plus 
    a.cf_plus
  else if sym==:minus
    a.cf_minus
  else 
    getfield(a,sym)
  end
end 

[amartinhuertas]

Solved in PRs #797 #799 and #800