diff --git a/src/ReferenceFEs/Pullbacks.jl b/src/ReferenceFEs/Pullbacks.jl index d339d446b..f9f8ad81f 100644 --- a/src/ReferenceFEs/Pullbacks.jl +++ b/src/ReferenceFEs/Pullbacks.jl @@ -31,7 +31,7 @@ end function Arrays.lazy_map( ::Broadcasting{typeof(gradient)}, a::LazyArray{<:Fill{Broadcasting{Operation{<:Pushforward}}}} ) - cell_ref_fields, args = a.args + cell_ref_fields, args... = a.args cell_ref_gradient = lazy_map(Broadcasting(∇),cell_ref_fields) return lazy_map(a.maps.value,cell_ref_gradient,args...) end @@ -123,15 +123,16 @@ function Arrays.lazy_map( ::typeof(evaluate),k::LazyArray{<:Fill{<:Pullback}},ref_cell_fields::AbstractArray ) pb = k.maps.value - phys_cell_dofs, pf_args = k.args + phys_cell_dofs, pf_args... = k.args phys_cell_fields = lazy_map(pb.pushforward,ref_cell_fields,pf_args...) return lazy_map(evaluate,phys_cell_dofs,phys_cell_fields) end function evaluate!( - cache, k::PullBack, σ_phys::MomentBasedDofBasis, args... + cache, k::Pullback, σ_phys::AbstractVector{<:Dof}, args... ) - Broadcasting(Operation(k))(f_phys,args...) + pf(f_ref) = evaluate(k.pushforward,f_ref,args...) # TODO: Can we avoid this? + return Arrays.OperationMap(σ_phys,pf) end # InversePullback @@ -145,7 +146,7 @@ where - V̂* is a dof space on the reference cell K̂ and - V* is a dof space on the physical cell K. Its action on reference dofs ̂σ : V̂ -> R is defined in terms of the pushforward map F* as - σ = F**(̂σ) := ̂σ∘(F*)^-1 : V -> R + σ = (F**)^-1(̂σ) := ̂σ∘(F*)^-1 : V -> R """ struct InversePullback{PF} <: Map pushforward::PF @@ -162,11 +163,18 @@ function Arrays.lazy_map( ::typeof(evaluate),k::LazyArray{<:Fill{<:InversePullback}},phys_cell_fields::AbstractArray ) pb = inverse_map(k.maps.value) - ref_cell_dofs, pf_args = k.args + ref_cell_dofs, pf_args... = k.args ref_cell_fields = lazy_map(inverse_map(pb.pushforward),phys_cell_fields,pf_args...) return lazy_map(evaluate,ref_cell_dofs,ref_cell_fields) end +function evaluate!( + cache, k::InversePullback, σ_ref::AbstractVector{<:Dof}, args... +) + ipf(f_phys) = evaluate(inverse_map(k.pushforward),f_phys,args...) # TODO: Can we avoid this? + return Arrays.OperationMap(σ_ref,ipf) +end + # ContraVariantPiolaMap struct ContraVariantPiolaMap <: Pushforward end diff --git a/test/pullbacks.jl b/test/pullbacks.jl index a59f851b3..ab9573288 100644 --- a/test/pullbacks.jl +++ b/test/pullbacks.jl @@ -1,7 +1,8 @@ using FillArrays using Gridap -using Gridap.ReferenceFEs, Gridap.FESpaces +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.CellData +using Gridap.Fields using Gridap.ReferenceFEs: Pullback @@ -29,15 +30,25 @@ Arr = lazy_map(evaluate,σ_ref,φ_ref)[1] pf = ContraVariantPiolaMap() +# Inverse Pushforward ipf = inverse_map(pf) f_ref = lazy_map(ipf,φ_phys,Jt,sign) Brr = lazy_map(evaluate,σ_ref,f_ref)[1] Brr == Arr φ_ref_x = lazy_map(evaluate,φ_ref,pts)[1] f_ref_x = lazy_map(evaluate,f_ref,pts)[1] -f_ref_x == φ_ref_x +f_ref_x ≈ φ_ref_x +# Pullback pb = Pullback(pf) θ_ref = lazy_map(pb,σ_phys,Jt,sign) +θ_ref_x = lazy_map(evaluate,θ_ref,φ_ref) +σ_ref_x = lazy_map(evaluate,σ_ref,φ_ref) +θ_ref_x ≈ σ_ref_x +# Inverse Pullback ipb = inverse_map(pb) +θ_phys = lazy_map(ipb,σ_ref,Jt,sign) +θ_phys_x = lazy_map(evaluate,θ_phys,φ_phys) +σ_phys_x = lazy_map(evaluate,σ_phys,φ_phys) +θ_phys_x ≈ σ_phys_x