Changed the method for SkeletonCellFieldPair for AD to work with operations inside mean and jump

NEWS.md

@@ -10,6 +10,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added functionality to take gradient of functional involving integration (`DomainContribution`) over Skeleton faces, with respect to the degrees-of-freedom `FEFunction`.  The interface remains the same - `gradient(f,uh)`. Since PR [#797](https://github.com/gridap/Gridap.jl/pull/797)
 - Extended the `MultiField` functional gradient (with respect to degrees-of-freedom of `MultiFieldFEFunction`) to functionals involving Skeleton integration. The interface remains the same `gradient(f,xh)`. Since PR [#799](https://github.com/gridap/Gridap.jl/pull/799)
 
+### Fixed
+- Fixed the behavior of `gradient` for functionals involving operations of `CellFields` inside `mean` and `jump` terms of Skeleton Integration terms. Since PR [#800](https://github.com/gridap/Gridap.jl/pull/800)
+- Fixed the behavior `SkeletonCellFieldPair` at the Boundary integration terms. Since PR [#800](https://github.com/gridap/Gridap.jl/pull/800)
+
 ## [0.17.13] - 2022-05-31
 
 ### Added

src/CellData/SkeletonCellFieldPair.jl

@@ -68,10 +68,6 @@ function get_triangulation(a::SkeletonCellFieldPair)
 end
 
 #=
-The below change_domain has been overloaded for SkeletonCellFieldPair
-to work without errors for integrations over other triangulation - Ω
-(BodyFittedTriangulation) and Γ (BoundaryTriangulation).
-
 We arbitrarily choose plus side of SkeletonCellFieldPair for evaluations, as
 we don't use the DomainContributions associated with Ω and Γ while performing
 the sensitivities for SkeletonTriangulation (Λ) DomainContribution. This
@@ -80,12 +76,28 @@ is passed into it. Ideally, if we could parse and extract only the Skeleton
 integration terms from the functional's julia function form, this fix is not
 required, but this is not trivial to do. On the positive side, since the
 evaluations are all lazy and not used, this doesn't put any noticable memory
-or computational overhead.
+or computational overhead. Ofcourse, it is made sure that the such plus side
+pick doesn't happen when the integration over the SkeletonTriangulation
 =#
-function change_domain(a::SkeletonCellFieldPair,target_trian::Triangulation,target_domain::DomainStyle)
-  change_domain(a.plus,DomainStyle(a),target_trian,target_domain)
+# If SkeletonCellFieldPair is evaluated we just pick the plus side parent
+function evaluate!(cache,f::SkeletonCellFieldPair,x::CellPoint)
+  _f, _x = _to_common_domain(f.plus.parent,x)
+  cell_field = get_data(_f)
+  cell_point = get_data(_x)
+  lazy_map(evaluate,cell_field,cell_point)
+end
+
+# fix for CellFieldAt{T}(parent::OperationCellField) giving right output
+# or else it chooses the SkeletonCellFieldPair directly as the parent
+# resulting in errors and it not the intended behaviour to have the parent
+# as SkeletonCellFieldPair and is not consistent with our getproperty rules
+function CellFieldAt{T}(parent::SkeletonCellFieldPair) where T
+  getproperty(parent,T)
 end
 
+# to handle the evaluation of SkeletonCellFieldPair at BoundaryTriangulation
+get_data(a::SkeletonCellFieldPair) = get_data(a.plus)
+
 function Base.propertynames(a::SkeletonCellFieldPair, private::Bool=false)
   (fieldnames(typeof(a))...,:⁺,:plus,:⁻,:minus)
 end

test/FESpacesTests/FEAutodiffTests.jl

@@ -1,5 +1,6 @@
 module FEAutodiffTests
 
+using Test
 using LinearAlgebra
 using Gridap.Algebra
 using Gridap.FESpaces
@@ -120,6 +121,13 @@ a_Λ(u) = ∫( - jump(u*n_Λ)⊙mean(∇(u))
             - mean(∇(u))⊙jump(u*n_Λ)
             + jump(u*n_Λ)⊙jump(u*n_Λ) )dΛ
 
+# functionals having mean and jump of products of FEFunctions/CellFields
+g_Λ(uh) = ∫(mean(uh*uh))*dΛ
+h_Λ(uh) = ∫(jump(uh*uh*n_Λ)⋅jump(uh*uh*n_Λ))*dΛ
+j_Λ(uh) = ∫(mean(∇(uh)*uh)⋅jump((∇(uh)⋅∇(uh))*n_Λ))*dΛ
+
+@test sum(∫(mean(uh*uh))*dΛ) == 0.5*sum(∫(uh.plus*uh.plus + uh.minus*uh.minus)*dΛ)
+
 function f_uh_free_dofs(f,uh,θ)
   dir = similar(uh.dirichlet_values,eltype(θ))
   uh = FEFunction(U,θ,dir)
@@ -138,4 +146,20 @@ gridapgrada = assemble_vector(gradient(a_Λ,uh),U)
 fdgrada = ForwardDiff.gradient(a_Λ_,θ)
 test_array(gridapgrada,fdgrada,≈)
 
+g_Λ_(θ) = f_uh_free_dofs(g_Λ,uh,θ)
+h_Λ_(θ) = f_uh_free_dofs(h_Λ,uh,θ)
+j_Λ_(θ) = f_uh_free_dofs(j_Λ,uh,θ)
+
+gridapgradg = assemble_vector(gradient(g_Λ,uh),U)
+fdgradg = ForwardDiff.gradient(g_Λ_,θ)
+test_array(gridapgradg,fdgradg,≈)
+
+gridapgradh = assemble_vector(gradient(h_Λ,uh),U)
+fdgradh = ForwardDiff.gradient(h_Λ_,θ)
+test_array(gridapgradh,fdgradh,≈)
+
+gridapgradj = assemble_vector(gradient(j_Λ,uh),U)
+fdgradj = ForwardDiff.gradient(j_Λ_,θ)
+test_array(gridapgradj,fdgradj,≈)
+
 end # module

test/MultiFieldTests/MultiFieldFEAutodiffTests.jl

@@ -156,7 +156,7 @@ dΛ = Measure(Λ,2)
 n_Λ = get_normal_vector(Λ)
 
 g_Λ((uh,ph)) = ∫( mean(uh) + mean(ph) + mean(uh)*mean(ph) )dΛ
-# f_Λ((uh,ph)) = ∫( mean(uh*uh) + mean(uh*ph) + mean(ph*ph) )dΛ
+f_Λ((uh,ph)) = ∫( mean(uh*uh) + mean(uh*ph) + mean(ph*ph) )dΛ
 a_Λ((uh,ph)) = ∫( - jump(uh*n_Λ)⊙mean(∇(ph))
                   - mean(∇(uh))⊙jump(ph*n_Λ)
                   + jump(uh*n_Λ)⊙jump(ph*n_Λ) )dΛ
@@ -175,13 +175,13 @@ end
 # can also do the above by constructing a MultiFieldCellField
 
 g_Λ_(θ) = f_uh_free_dofs(g_Λ,xh,θ)
-# f_Λ_(θ) = f_uh_free_dofs(f_Λ,xh,θ)
+f_Λ_(θ) = f_uh_free_dofs(f_Λ,xh,θ)
 a_Λ_(θ) = f_uh_free_dofs(a_Λ,xh,θ)
 
 θ = get_free_dof_values(xh)
 
 # check if the evaluations are working
-# @test sum(f_Λ(xh)) == f_Λ_(θ)
+@test sum(f_Λ(xh)) == f_Λ_(θ)
 @test sum(a_Λ(xh)) == a_Λ_(θ)
 @test sum(g_Λ(xh)) == g_Λ_(θ)
 
@@ -191,6 +191,10 @@ gridapgradf_Ω = assemble_vector(gradient(f_Ω,xh),X)
 fdgradf_Ω = ForwardDiff.gradient(f_Ω_,θ)
 test_array(gridapgradf_Ω,fdgradf_Ω,≈)
 
+gridapgradf = assemble_vector(gradient(f_Λ,xh),X)
+fdgradf = ForwardDiff.gradient(f_Λ_,θ)
+test_array(gridapgradf,fdgradf,≈)
+
 gridapgradg = assemble_vector(gradient(g_Λ,xh),X)
 fdgradg = ForwardDiff.gradient(g_Λ_,θ)
 test_array(gridapgradg,fdgradg,≈)