Primal HHO for Poisson problem working!

src/GridapAPIExtensions.jl

@@ -626,7 +626,7 @@ function Gridap.Arrays.return_cache(
   A::Gridap.Fields.MatrixBlock,
   B::Gridap.Fields.MatrixBlock)
 
-  Gridap.Helpers.@check size(A.array) == size(B.array)
+  Gridap.Helpers.@check size(A.array)==size(B.array)
 
   nA=findall(A.touched)
   nB=findall(B.touched)
@@ -708,7 +708,7 @@ function Gridap.Arrays.return_cache(
   B::Gridap.Fields.VectorBlock{<:Matrix})
 
   Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
-  Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
 
   nA=findall(A.touched)
   nB=findall(B.touched)
@@ -735,7 +735,53 @@ function Gridap.Arrays.evaluate!(
   B::Gridap.Fields.VectorBlock{<:Matrix})
   r,c=cache
   Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
+  Gridap.Helpers.@check length(findall(A.touched)) == 1
+  Gridap.Helpers.@check length(findall(B.touched)) == 1
+  Gridap.Helpers.@check r.touched[1]
+  ai=testitem(A)
+  bi=testitem(B)
+  r.array[1]=evaluate!(c,k,ai,bi)
+  r
+end
+
+# A [f,f]      , e.g., [2,2] nonzero block (matrix to be inverted)
+# B [f,1]      , e.g., [2,1] nonzero block (several RHS)
+# Return [1]
+function Gridap.Arrays.return_cache(
+  k::typeof(compute_bulk_to_skeleton_l2_projection_dofs),
+  A::Gridap.Fields.MatrixBlock{<:Matrix},
+  B::Gridap.Fields.MatrixBlock{<:Matrix})
+
+  Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
+
+  nA=findall(A.touched)
+  nB=findall(B.touched)
+
+  Gridap.Helpers.@check length(nA)==length(nB)
+  Gridap.Helpers.@check length(nA)==1
+  Gridap.Helpers.@check nA[1][1]==nB[1][1]
+  Gridap.Helpers.@check nA[1][2]==nB[1][1]
+
+  ai=testitem(A)
+  bi=testitem(B)
+  T=Gridap.Arrays.return_type(k,ai,bi)
+  touched  = Vector{Bool}(undef,1)
+  touched .= true
+  r = Vector{T}(undef,1)
+  c = Gridap.Arrays.return_cache(k,ai,bi)
+  Gridap.Fields.ArrayBlock(r,touched), c
+end
+
+function Gridap.Arrays.evaluate!(
+  cache,
+  k::typeof(compute_bulk_to_skeleton_l2_projection_dofs),
+  A::Gridap.Fields.MatrixBlock{<:Matrix},
+  B::Gridap.Fields.MatrixBlock{<:Matrix})
+  r,c=cache
   Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
   Gridap.Helpers.@check length(findall(A.touched)) == 1
   Gridap.Helpers.@check length(findall(B.touched)) == 1
   Gridap.Helpers.@check r.touched[1]

src/HybridAffineFEOperators.jl

@@ -70,8 +70,14 @@ function Gridap.FESpaces.solve(op::HybridAffineFEOperator)
 end
 
 function Gridap.FESpaces.solve!(uh, solver::LinearFESolver, op::HybridAffineFEOperator, cache)
-  # Solve linear system defined on the skeleton
+   # Solve linear system defined on the skeleton.
+   # Some solvers do not support 0x0 matrices. Avoid
+   # calling solve in such an scenario.
+   if (size(op.skeleton_op.op.matrix) != (0,0))
     lh = solve(op.skeleton_op)
+   else
+    lh = FEFunction(op.skeleton_op.trial,zeros(num_free_dofs(op.skeleton_op.trial)))
+   end
 
   # Invoke weak form of the hybridizable system
     u = get_trial_fe_basis(op.trial)

src/LocalFEOperators.jl

@@ -267,46 +267,6 @@ function _generate_image_space_span(op::LocalFEOperator,
   Gridap.MultiField.MultiFieldCellField(all_febases)
 end
 
-# TO-DO: to study re-use opportunities among the following function and
-# function Gridap.FESpaces._change_argument(
-#  op::typeof(jacobian),f,trian::SkeletonTriangulation,uh::Gridap.MultiField.MultiFieldFEFunction)
-# function _generate_image_space_span(O::Gridap.FESpaces.SingleFieldFESpace,
-#                                     v::Gridap.MultiField.MultiFieldCellField,
-#                                     cell_dofs,
-#                                     basis_style,
-#                                     skeleton::SkeletonTriangulation)
-
-#   fields_on_O_triangulation=_generate_image_space_span(O,
-#                                                        v,
-#                                                        cell_dofs,
-#                                                        basis_style,
-#                                                        get_triangulation(O))
-#   nfields=length(v.single_fields)
-#   single_fields = Gridap.CellData.CellField[]
-#   for i=1:nfields
-#     du_i_b = fields_on_O_triangulation.single_fields[i]
-#     vi     = v.single_fields[i].single_field
-#     if (_is_on_skeleton(vi))
-#       m = get_background_model(get_triangulation(vi))
-#       cell_wise_facets = _get_cell_wise_facets(m)
-#       cell_vector=Gridap.CellData.get_data(du_i_b)
-#       cell_vector_split=_cell_vector_facets_split(vi)
-#       trian = skeleton
-#       sglue = Gridap.Geometry.get_glue(trian,Val(num_cell_dims(m)))
-#       cell_field=SkeletonExpandedVectorFromSplitCellVector(trian.glue,
-#                                                            cell_wise_facets,
-#                                                            cell_vector,
-#                                                            cell_vector_split)
-#       cell_field=lazy_map(Broadcasting(∘),cell_field,sglue.tcell_lface_mface_map)
-#       cf = Gridap.CellData.GenericCellField(cell_field, trian, DomainStyle(du_i_b))
-#     else
-#       cf = v.single_fields[i]
-#     end
-#     push!(single_fields, cf)
-#   end
-#   Gridap.MultiField.MultiFieldCellField(single_fields)
-# end
-
 function _compute_cell_dofs_field_offsets(U::Gridap.MultiField.MultiFieldFESpace)
   nfields = length(U.spaces)
   cell_dofs_field_offsets=Vector{Int}(undef,nfields+1)
@@ -328,7 +288,7 @@ end
 
 function _generate_image_space_span(op::LocalFEOperator{<:SingleValued},
                                     O::Gridap.MultiField.MultiFieldFESpace,
-                                    v::Gridap.MultiField.MultiFieldCellField,
+                                    v::Gridap.CellData.CellField,
                                     cell_dofs,
                                     basis_style) where N
   cell_dofs_field_offsets=_compute_cell_dofs_field_offsets(O)
@@ -392,6 +352,7 @@ function _generate_image_space_span(op::LocalFEOperator{<:MultiValued},
         cell_field=lazy_map(linear_combination,skel_facet_dofs,skel_facet_shapefuns)
         cf = Gridap.CellData.GenericCellField(cell_field, tskel, DomainStyle(mf_basis[i]))
         push!(single_fields,cf)
+        sj=ej+1
       end
       multi_field=Gridap.MultiField.MultiFieldCellField(single_fields)
     else

test/LocalFEOperatorTests.jl

@@ -10,30 +10,6 @@ module LocalFEOperatorTests
     LocalFEOperator((m,n),UK_U∂K,VK_V∂K)
   end
 
-  function setup_l2_projection_operator(UK_U∂K,VK_V∂K,dΩ,d∂K)
-    m((uK,u∂K)  , (vK,v∂K)) = ∫(vK*uK)dΩ  + ∫(v∂K*u∂K)d∂K
-    n((uhK,uh∂K), (vK,v∂K)) = ∫(vK*uhK)dΩ + ∫(v∂K*uh∂K)d∂K
-    LocalFEOperator((m,n),UK_U∂K,VK_V∂K)
-  end
-
-  function setup_difference_operator(l2_projection_op,reconstruction_op)
-     function _op(uK_u∂K)
-        u_rec=reconstruction_op(uK_u∂K)
-        uhK,uh∂K=uK_u∂K
-        l2_projection_op((u_rec-uhK,u_rec-uh∂K))
-     end
-     return _op
-  end
-
-  function setup_skeleton_difference_operator(reduction_op,reconstruction_op)
-    function _op(uK_u∂K)
-       uK,u∂K=uK_u∂K
-       u_rec=reconstruction_op(uK_u∂K)
-       reduction_op((u_rec-uK))
-    end
-    return _op
- end
-
   function setup_reconstruction_operator(model, order, dΩ, d∂K)
     T         = Float64
     nK        = get_cell_normal_vector(d∂K.quad.trian)
@@ -65,6 +41,18 @@ module LocalFEOperatorTests
                     test_space_ds_decomp=VKR_DS_DECOMP)
   end
 
+  function setup_difference_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K)
+    m((uK,u∂K)  , (vK,v∂K)) = ∫(vK*uK)dΩ + ∫(v∂K*u∂K)d∂K
+    function n(uK_u∂K, (vK,v∂K))
+      uK_u∂K_rec     = R(uK_u∂K)
+      uK,u∂K         = uK_u∂K
+      uK_rec,u∂K_rec = uK_u∂K_rec
+      ∫(vK *uK_rec)dΩ  + ∫(vK*u∂K_rec)dΩ   - ∫(vK *uK)dΩ +
+      ∫(v∂K*uK_rec)d∂K + ∫(v∂K*u∂K_rec)d∂K - ∫(v∂K*u∂K)d∂K
+    end
+    LocalFEOperator((m,n),UK_U∂K,VK_V∂K; field_type_at_common_faces=MultiValued())
+   end
+
   model=CartesianDiscreteModel((0,1,0,1),(2,2))
 
   # Geometry
@@ -91,13 +79,13 @@ module LocalFEOperatorTests
   nK     = get_cell_normal_vector(∂K)
   d∂K    = Measure(∂K,degree)
 
-  reconstruction_op=setup_reconstruction_operator(model, order, dΩ, d∂K)
+  R=setup_reconstruction_operator(model, order, dΩ, d∂K)
   uhK_uh∂K =get_trial_fe_basis(UK_U∂K)
-  reconstruction_op_image_span=reconstruction_op(uhK_uh∂K)
+  reconstruction_op_image_span=R(uhK_uh∂K)
   uh_dofs=zeros(num_free_dofs(UK_U∂K))
   uh_dofs[2]=1.0
   uhK_uh∂K=FEFunction(UK_U∂K,uh_dofs)
-  projected_uhK_uh∂K=reconstruction_op(uhK_uh∂K)
+  projected_uhK_uh∂K=R(uhK_uh∂K)
   # (uK,u∂K) = get_trial_fe_basis(UK_U∂K)
   # v_c,v_nc = get_fe_basis(VKR_DS_DECOMP)
   # u_nzm,u_zm = get_trial_fe_basis(UKR_DS_DECOMP)
@@ -111,7 +99,7 @@ module LocalFEOperatorTests
   y=reduction_op(uhk)
 
   # Check that the result of applying the reconstruction operator to the
-  # the result of applyting the reduction operator to P_K^{k+1} is P_K^{k+1}
+  # the result of applying the reduction operator to P_K^{k+1} is P_K^{k+1}
   # itself
   reffeᵤ     = ReferenceFE(lagrangian,Float64,order+1)
   VH1        = TestFESpace(Ω, reffeᵤ; conformity=:H1)
@@ -120,7 +108,89 @@ module LocalFEOperatorTests
   uh         = FEFunction(UH1, free_dofs)
   uh_reduced = reduction_op(uh)
 
-  uh_reconstructed = reconstruction_op(uh_reduced)
+  uh_reconstructed = R(uh_reduced)
   eh = uh_reconstructed-uh
   @test sum(∫(eh*eh)dΩ) < 1.0e-12
+
+  # Check that the mean value of the reconstructed cell FE space functions
+  # match the mean value of the original cell FE space functions
+  uhK_uh∂K =get_trial_fe_basis(UK_U∂K)
+  vhK_vh∂K =get_fe_basis(VK_V∂K)
+  vK,v∂K=vhK_vh∂K
+  uK,u∂K=uhK_uh∂K
+
+  R_vhK_vh∂K=R(vhK_vh∂K)
+  R_vhK,_=R_vhK_vh∂K
+
+  dc1=∫(vK)dΩ
+  dc2=∫(R_vhK)dΩ
+
+  @test all(get_array(dc1) .≈ get_array(dc2))
+
+  diff_op=setup_difference_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K)
+  v=get_fe_basis(VK_V∂K)
+  ub=get_trial_fe_basis(UK_U∂K)
+
+  ub_rec=R(ub)
+  ub_rec1,ub_rec2=ub_rec
+  xΩ=Gridap.CellData.get_cell_points(dΩ.quad)
+  x∂K=Gridap.CellData.get_cell_points(d∂K.quad)
+
+  v_rec=R(v)
+  v_rec1,v_rec2=v_rec
+  v_rec1_d∂K=Gridap.CellData.change_domain(v_rec1,d∂K.quad.trian,Gridap.CellData.ReferenceDomain())
+  v_rec1_d∂K(x∂K)[1][4][1][1]
+  v_rec2_d∂K=Gridap.CellData.change_domain(v_rec2,d∂K.quad.trian,Gridap.CellData.ReferenceDomain())
+  v_rec2_d∂K(x∂K)[1][3][2][1]
+
+  uK_u∂K_rec     = R(ub)
+  uK,u∂K         = ub
+  uK_rec,u∂K_rec = uK_u∂K_rec
+  dc=∫(vK *uK_rec)dΩ  + ∫(vK*u∂K_rec)dΩ   - ∫(vK *uK)dΩ +
+     ∫(v∂K*uK_rec)d∂K + ∫(v∂K*u∂K_rec)d∂K - ∫(v∂K*u∂K)d∂K
+
+
+  δvK,δv∂K=diff_op(v)
+  δuK,δu∂K=diff_op(ub)
+
+  δv∂K_K,δv∂K_∂K=δv∂K
+  δu∂K_K,δu∂K_∂K=δu∂K
+
+  δv∂K_K(x∂K)[1][1][1]
+  δv∂K_K(x∂K)[1][2][1]
+  δv∂K_K(x∂K)[1][3][1]
+  δv∂K_K(x∂K)[1][4][1]
+  δv∂K_∂K(x∂K)[1][1][2]
+  δv∂K_∂K(x∂K)[1][2][2]
+  δv∂K_∂K(x∂K)[1][3][2]
+  δv∂K_∂K(x∂K)[1][4][2]
+
+  (δv∂K_K*δu∂K_K)(x∂K)[1][4][1,1]
+  (δv∂K_K*δu∂K_∂K)(x∂K)[1][4][1,2]
+  (δv∂K_∂K*δu∂K_K)(x∂K)[1][4][2,1]
+  (δv∂K_∂K*δu∂K_∂K)(x∂K)[1][4][2,2]
+
+  dc=∫(δv∂K_K*δu∂K_K +
+      δv∂K_K*δu∂K_∂K+
+       δv∂K_∂K*δu∂K_K+
+         δv∂K_∂K*δu∂K_∂K)d∂K
+
+  get_array(dc)[1][1,1]
+  get_array(dc)[1][2,1]
+  get_array(dc)[1][1,2]
+  get_array(dc)[1][2,2]
+  δvK_K,δvK_∂K=δvK
+  δuK_K,δuK_∂K=δuK
+
+  dc=∫(δvK_K*δuK_K)dΩ+
+       ∫(δvK_K*δuK_∂K)dΩ+
+         ∫(δvK_∂K*δuK_K)dΩ+
+           ∫(δvK_∂K*δuK_∂K)dΩ
+
+  get_array(dc)[1][1,1]
+  get_array(dc)[1][1,2]
+  get_array(dc)[1][2,1]
+  get_array(dc)[1][2,2]
+
+
 end

test/P_m.jl

@@ -6,11 +6,11 @@ function Pₘ(uh, uhΓ, μ, d∂K)
   B = ∫(μ⋅uh)d∂K
   # [c][f][bμ,buhΓ==bμ][f,f]
   A_array = _remove_sum_facets(_remove_densify(Gridap.CellData.get_array(A)))
-  # [c][f][bμ,buh][f] (if uh Trial Basis)
-  # [c][f][bμ]    [f] (if uh FE Function)
+  # [c][f][bμ,buh][f,1] (if uh Trial Basis)
+  # [c][f][bμ]    [f]   (if uh FE Function)
   B_array = _remove_sum_facets(_remove_densify(Gridap.CellData.get_array(B)))
-  # [c][f][1,buh][1] (if uh Trial Basis)
-  # [c][f]           (if uh FE Function)
+  # [c][f][1,buh][1]   (if uh Trial Basis)
+  # [c][f]             (if uh FE Function)
   pm_uh_dofs = lazy_map(GridapHybrid.compute_bulk_to_skeleton_l2_projection_dofs, A_array, B_array)
   # [c][f][1,buhΓ][1,f]
   uhΓ_d∂K = Gridap.CellData.change_domain(uhΓ, d∂K.quad.trian, ReferenceDomain())