* Added convergence analysis to HHO driver

* Added stabilization

test/PoissonHHOTests.jl

@@ -11,6 +11,8 @@ module PoissonHHOTests
     reffe_c   = ReferenceFE(monomial_basis  , Float64, order+1; subspace=:OnlyConstant)
     reffe_nc  = ReferenceFE(monomial_basis  , Float64, order+1; subspace=:ExcludeConstant)
 
+    Ω = dΩ.quad.trian
+
     VKR     = TestFESpace(Ω  , refferecᵤ; conformity=:L2)
     UKR     = TrialFESpace(VKR)
     UKR_NZM = TrialFESpace(TestFESpace(Ω, reffe_nzm; conformity=:L2))
@@ -42,57 +44,107 @@ module PoissonHHOTests
     LocalFEOperator((m,n),UK_U∂K,VK_V∂K; field_type_at_common_faces=MultiValued())
    end
 
-  u(x)=x[1]+x[2]
-  f(x)=-Δ(u)(x)
-
-  model=CartesianDiscreteModel((0,1,0,1),(2,2))
-  D  = num_cell_dims(model)
-  Ω  = Triangulation(ReferenceFE{D},model)
-  Γ  = Triangulation(ReferenceFE{D-1},model)
-  ∂K = GridapHybrid.Skeleton(model)
-
-  order     = 1
-  reffeᵤ    = ReferenceFE(lagrangian,Float64,order  ;space=:P)
-
-  VK     = TestFESpace(Ω  , reffeᵤ; conformity=:L2)
-  V∂K    = TestFESpace(Γ  , reffeᵤ; conformity=:L2,dirichlet_tags=collect(5:8))
-  UK     = TrialFESpace(VK)
-  U∂K    = TrialFESpace(V∂K,u)
-  VK_V∂K = MultiFieldFESpace([VK,V∂K])
-  UK_U∂K = MultiFieldFESpace([UK,U∂K])
-
-  degree = 2*order+1
-  dΩ     = Measure(Ω,degree)
-  dΓ     = Measure(Γ,degree)
-  d∂K    = Measure(∂K,degree)
-
-  R=setup_reconstruction_operator(model, order, dΩ, d∂K)
-  diff_op=setup_difference_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K)
-
-  function r(u,v)
-    uK,u∂K=R(u)
-    vK,v∂K=R(v)
-    ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ
+   p = 3
+   u(x) = x[1]^p+x[2]^p
+   f(x)=-Δ(u)(x)
+
+   #u(x)=x[1]+x[2]
+   #f(x)=-Δ(u)(x)
+
+  function solve_hho(cells,order)
+      partition = (0,1,0,1)
+      model = CartesianDiscreteModel(partition, cells)
+      D  = num_cell_dims(model)
+      Ω  = Triangulation(ReferenceFE{D},model)
+      Γ  = Triangulation(ReferenceFE{D-1},model)
+      ∂K = GridapHybrid.Skeleton(model)
+
+      reffeᵤ    = ReferenceFE(lagrangian,Float64,order  ;space=:P)
+
+      VK     = TestFESpace(Ω  , reffeᵤ; conformity=:L2)
+      V∂K    = TestFESpace(Γ  , reffeᵤ; conformity=:L2,dirichlet_tags=collect(5:8))
+      UK     = TrialFESpace(VK)
+      U∂K    = TrialFESpace(V∂K,u)
+      VK_V∂K = MultiFieldFESpace([VK,V∂K])
+      UK_U∂K = MultiFieldFESpace([UK,U∂K])
+
+      degree = 2*order+1
+      dΩ     = Measure(Ω,degree)
+      dΓ     = Measure(Γ,degree)
+      d∂K    = Measure(∂K,degree)
+
+      R=setup_reconstruction_operator(model, order, dΩ, d∂K)
+      diff_op=setup_difference_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K)
+
+      function r(u,v)
+        uK,u∂K=R(u)
+        vK,v∂K=R(v)
+        ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ +
+          ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ
+      end
+
+      function s(u,v)
+        h_T=CellField(get_array(∫(1)dΩ),Ω)
+        h_T_1=1.0/h_T
+        h_T_2=1.0/(h_T*h_T)
+
+        δuK,δu∂K=diff_op(u)
+        δvK,δv∂K=diff_op(v)
+        δvK_K,δvK_∂K=δvK
+        δuK_K,δuK_∂K=δuK
+        δv∂K_K,δv∂K_∂K=δv∂K
+        δu∂K_K,δu∂K_∂K=δu∂K
+        ∫(h_T_2*δvK_K*δuK_K)dΩ+
+          ∫(h_T_2*δvK_K*δuK_∂K)dΩ+
+             ∫(h_T_2*δvK_∂K*δuK_K)dΩ+
+               ∫(h_T_2*δvK_∂K*δuK_∂K)dΩ+
+        ∫(h_T_1*δv∂K_∂K*δu∂K_∂K)d∂K
+      end
+
+      a(u,v)=r(u,v)+s(u,v)
+      l((vK,))=∫(vK*f)dΩ
+
+
+      @time op=HybridAffineFEOperator((u,v)->(a(u,v),l(v)), UK_U∂K, VK_V∂K, [1], [2])
+      @time xh=solve(op)
+
+      uhK,uh∂K=xh
+      e = u -uhK
+      # @test sqrt(sum(∫(e⋅e)dΩ)) < 1.0e-12
+      return sqrt(sum(∫(e⋅e)dΩ))
   end
 
-  function s(u,v)
-    δuK,δu∂K=diff_op(u)
-    δvK,δv∂K=diff_op(v)
-    δvK_K,δvK_∂K=δvK
-    δuK_K,δuK_∂K=δuK
-    δv∂K_K,δv∂K_∂K=δv∂K
-    δu∂K_K,δu∂K_∂K=δu∂K
-    ∫(δvK_K*δuK_K)dΩ+∫(δvK_K*δuK_∂K)dΩ+∫(δvK_∂K*δuK_K)dΩ+∫(δvK_∂K*δuK_∂K)dΩ+
-    ∫(δv∂K_∂K*δu∂K_∂K)d∂K
+  function conv_test(ns,order)
+    el2 = Float64[]
+    hs = Float64[]
+    for n in ns
+      l2 = solve_hho((n,n),order)
+      println(l2)
+      h = 1.0/n
+      push!(el2,l2)
+      push!(hs,h)
+    end
+    println(el2)
+    el2
   end
 
-  a(u,v)=r(u,v)+s(u,v)
-  l((vK,))=∫(vK*f)dΩ
-
-  op=HybridAffineFEOperator((u,v)->(a(u,v),l(v)), UK_U∂K, VK_V∂K, [1], [2])
-  xh=solve(op)
+  function slope(hs,errors)
+    x = log10.(hs)
+    y = log10.(errors)
+    linreg = hcat(fill!(similar(x), 1), x) \ y
+    linreg[2]
+  end
 
-  uhK,uh∂K=xh
-  e = u -uhK
-  @test sqrt(sum(∫(e⋅e)dΩ)) < 1.0e-12
+  ns=[8,16,32,64,128]
+  order=1
+  el, hs = conv_test(ns,order)
+  println("Slope L2-norm u: $(slope(hs,el))")
+  slopek  =[Float64(ni)^(-(order)) for ni in ns]
+  slopekp1=[Float64(ni)^(-(order+1)) for ni in ns]
+  slopekp2=[Float64(ni)^(-(order+2)) for ni in ns]
+  plot(hs,[el slopek slopekp1 slopekp2],
+    xaxis=:log, yaxis=:log,
+    label=["L2u (measured)" "slope k" "slope k+1" "slope k+2"],
+    shape=:auto,
+    xlabel="h",ylabel="L2 error",legend=:bottomright)
 end

test/StressAssistedHenckyHDGTests.jl

@@ -165,31 +165,16 @@ function solve_stress_assisted_diffusion_hencky_hdg(cells,order;write_results=fa
     MuΓ_TR = TrialFESpace(MuΓ,u₀)
     Y_TR   = MultiFieldFESpace([W_TR, W_TR, Ψ_TR, S_TR, S_TR, V_TR, MϕΓ_TR, MuΓ_TR])
 
-    # FE formulation params
-    degree = 2 * (order + 1) # TO-DO: To think which is the minimum degree required
-    dΩ = Measure(Ω, degree)
-    n = get_cell_normal_vector(∂K)
+    # Integration
+    degree = 2 * (order + 1)
+    dΩ  = Measure(Ω, degree)
     d∂K = Measure(∂K, degree)
 
-    # Stabilization operator (solute concentration)
-    τρΓ = 1.0 # To be defined???
-    # Stabilization operator (displacements)
-    τuΓ = Float64(cells[1]) # To be defined???
-
-    println("h=$(1.0/cells[1]) τuΓ=$(τuΓ) τρΓ=$(τρΓ)")
-
-    Y_basis    = get_fe_basis(Y)
-    Y_TR_basis = get_trial_fe_basis(Y_TR)
-    (_,_,_,_,_,_,_,uhΓ_basis) = Y_TR_basis
+    # Outer normal, cell boundary
+    n = get_cell_normal_vector(∂K)
 
-    # # Testing for correctness of Pₘ projection
-    # xh              = FEFunction(Y_TR,rand(num_free_dofs(Y_TR)))
-    # _,_,_,_,_,_,_,μ = Y_basis
-    # _, uh, _, _ = xh
-    # dc=∫(μ⋅(uh-Pₘ(uh, uhΓ_basis, μ, d∂K)))d∂K
-    # for k in dc.dict[d∂K.quad.trian]
-    #   @test all(k[8] .< 1.0e-14)
-    # end
+    # Stabilization operator (displacements)
+    τuΓ = Float64(cells[1])
 
     # Residual
     # current_iterate  => (ωh,ρh,ϕh,th,σh,uh,ϕhΓ,uhΓ)
@@ -207,16 +192,7 @@ function solve_stress_assisted_diffusion_hencky_hdg(cells,order;write_results=fa
       ∫( μh⋅(σh⋅n) )d∂K - ∫( τuΓ*μh⋅Pm_uh)d∂K + ∫( τuΓ*μh⋅uhΓ )d∂K
     end
 
-    # free_values = zeros(num_free_dofs(Y_TR))
-    # xh          = FEFunction(Y_TR,free_values)
-    # (ωh,ρh,ϕh,th,σh,uh,ϕhΓ,uhΓ) = xh
-    # (rh,nh,ψh,sh,τh,vh,ηh,μh)   = Y_basis
-    # dc=residual((ωh,ρh,ϕh,th,σh,uh,ϕhΓ,uhΓ),(rh,nh,ψh,sh,τh,vh,ηh,μh))
-    # res=assemble_vector(dc,Y)
-    # dc=jacobian(x->residual(x,Y_basis),xh)
-    # A=assemble_matrix(dc,Y_TR,Y)
     @time op=HybridFEOperator(residual,Y_TR,Y,collect(1:6),[7,8])
-    # op=FEOperator(residual,Y_TR,Y)
 
     nls = NLSolver(show_trace=true, method=:newton, ftol=1.0e-14)
     solver = FESolver(nls)
@@ -259,8 +235,6 @@ function solve_stress_assisted_diffusion_hencky_hdg(cells,order;write_results=fa
      end
 
      norm2_ω,norm2_ρ,norm2_t,norm2_ϕ,norm2_σ,norm2_u
-    #@test sqrt(sum(∫(eu ⋅ eu)dΩ)) < 1.0e-12
-    #@test sqrt(sum(∫(eσ ⊙ eσ)dΩ)) < 1.0e-12
   end