Added new tutorials

Author: JordiManyer

Project.toml

@@ -17,6 +17,7 @@ GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapMakie = "41f30b06-6382-4b60-a5f7-79d86b35bf5d"
 GridapP4est = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9"
 GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
+GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
@@ -35,11 +36,12 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 DataStructures = "0.18.22"
-Gridap = "0.18"
+Gridap = "0.19"
 GridapDistributed = "0.4"
 GridapGmsh = "0.7"
 GridapP4est = "0.3"
 GridapPETSc = "0.5"
+GridapSolvers = "0.5.0"
 MPI = "0.20"
 SpecialFunctions = "1"
 julia = "1.6"

src/poisson_hdg.jl

@@ -0,0 +1,170 @@
+# # [Tutorial: Hybridizable Discontinuous Galerkin Method for the Poisson equation]
+#
+# In this tutorial, we will implement a Hybridizable Discontinuous Galerkin (HDG) method
+# for solving the Poisson equation. The HDG method is an efficient variant of DG methods
+# that introduces an auxiliary variable λ on mesh interfaces to reduce the global system size.
+#
+# ## The Poisson Problem
+#
+# We consider the Poisson equation with Dirichlet boundary conditions:
+#
+# ```math
+# \begin{aligned}
+# -\Delta u &= f \quad \text{in} \quad \Omega\\
+# u &= g \quad \text{in} \quad \partial\Omega
+# \end{aligned}
+# ```
+#
+# ## HDG Discretization
+#
+# The HDG method first rewrites the problem as a first-order system:
+#
+# ```math
+# \begin{aligned}
+# \boldsymbol{q} + \nabla u &= \boldsymbol{0} \quad \text{in} \quad \Omega\\
+# \nabla \cdot \boldsymbol{q} &= -f \quad \text{in} \quad \Omega\\
+# u &= g \quad \text{on} \quad \partial\Omega
+# \end{aligned}
+# ```
+#
+# The HDG discretization introduces three variables:
+# - ``\boldsymbol{q}_h``: the approximation to the flux ``\boldsymbol{q}``
+# - ``u_h``: the approximation to the solution ``u``
+# - ``\lambda_h``: the approximation to the trace of ``u`` on element faces
+#
+# The method is characterized by:
+# 1. Local discontinuous approximation for ``(\boldsymbol{q}_h,u_h)`` in each element
+# 2. Continuous approximation for ``\lambda_h`` on element faces
+# 3. Numerical flux ``\widehat{\boldsymbol{q}}_h = \boldsymbol{q}_h + \tau(u_h - \lambda_h)\boldsymbol{n}``
+#
+# where τ is a stabilization parameter.
+#
+# First, let's load the required Gridap packages:
+
+using Gridap
+using Gridap.Geometry
+using Gridap.FESpaces
+using Gridap.MultiField
+using Gridap.CellData
+
+# ## Manufactured Solution
+#
+# We use the method of manufactured solutions to verify our implementation.
+# We choose a solution u and derive the corresponding source term f:
+
+u(x) = sin(2*π*x[1])*sin(2*π*x[2])
+q(x) = -∇(u)(x)  # Define the flux q = -∇u
+f(x) = (∇ ⋅ q)(x) # Source term f = -Δu = -∇⋅(∇u)
+
+# ## Mesh Generation
+#
+# Create a 3D simplicial mesh from a Cartesian grid:
+
+D = 2  # Problem dimension
+nc = Tuple(fill(8, D))  # 4 cells in each direction
+domain = Tuple(repeat([0, 1], D))  # Unit cube domain
+model = simplexify(CartesianDiscreteModel(domain,nc))
+
+# ## Volume and Interface Triangulations
+#
+# HDG methods require two types of meshes:
+# 1. Volume mesh for element interiors
+# 2. Skeleton mesh for element interfaces
+#
+# We also need patch-wise triangulations for local computations:
+
+Ω = Triangulation(ReferenceFE{D}, model)  # Volume triangulation
+Γ = Triangulation(ReferenceFE{D-1}, model)  # Skeleton triangulation
+
+ptopo = Geometry.PatchTopology(model)
+Ωp = Geometry.PatchTriangulation(model,ptopo)  # Patch volume triangulation
+Γp = Geometry.PatchBoundaryTriangulation(model,ptopo)  # Patch skeleton triangulation
+
+# ## Reference Finite Elements
+#
+# HDG uses three different finite element spaces:
+# 1. Vector-valued space for the flux q (Q)
+# 2. Scalar space for the solution u (V)
+# 3. Scalar space for the interface variable λ (M)
+
+order = 1  # Polynomial order
+reffe_Q = ReferenceFE(lagrangian, VectorValue{D, Float64}, order; space=:P)
+reffe_V = ReferenceFE(lagrangian, Float64, order; space=:P)
+reffe_M = ReferenceFE(lagrangian, Float64, order; space=:P)
+
+# ## Test and Trial Spaces
+#
+# Create discontinuous cell spaces for volume variables (q,u) and 
+# a discontinous face space for the interface variable λ:
+
+# Test spaces
+V = TestFESpace(Ω, reffe_V; conformity=:L2)  # Discontinuous vector space
+Q = TestFESpace(Ω, reffe_Q; conformity=:L2)  # Discontinuous scalar space
+M = TestFESpace(Γ, reffe_M; conformity=:L2, dirichlet_tags="boundary")  # Interface space
+
+# Only the skeleton space has Dirichlet BC
+N = TrialFESpace(M, u)
+
+# ## MultiField Structure
+#
+# Group the spaces for q, u, and λ using MultiField:
+
+mfs = MultiField.BlockMultiFieldStyle(2,(2,1))  # Special blocking for efficient static condensation
+X = MultiFieldFESpace([V, Q, N];style=mfs)
+
+# ## Integration Setup
+#
+# Define measures for volume and face integrals:
+
+degree = 2*(order+1)  # Integration degree
+dΩp = Measure(Ωp,degree)  # Volume measure
+dΓp = Measure(Γp,degree)  # Surface measure
+
+# ## HDG Parameters
+#
+# The stabilization parameter τ affects the stability and accuracy of the method:
+
+τ = 1.0 # HDG stabilization parameter
+
+# ## Weak Form
+#
+# The HDG weak form consists of three equations coupling q, u, and λ.
+# We need operators to help define the weak form:
+
+n = get_normal_vector(Γp)  # Face normal vector
+Πn(u) = u⋅n  # Normal component
+Π(u) = change_domain(u,Γp,DomainStyle(u))  # Project to skeleton
+
+# The bilinear and linear forms are:
+# 1. Volume integrals for flux and primal equations
+# 2. Interface integrals for hybridization
+# 3. Stabilization terms with parameter τ
+
+a((uh,qh,sh),(vh,wh,lh)) = ∫( qh⋅wh - uh*(∇⋅wh) - qh⋅∇(vh) )dΩp + ∫(sh*Πn(wh))dΓp +
+                           ∫((Πn(qh) + τ*(Π(uh) - sh))*(Π(vh) + lh))dΓp
+l((vh,wh,lh)) = ∫( f*vh )dΩp
+
+# ## Static Condensation and Solution
+#
+# A key feature of HDG is static condensation - eliminating volume variables
+# to get a smaller system for λ only:
+
+op = MultiField.StaticCondensationOperator(ptopo,X,a,l)
+uh, qh, sh = solve(op)
+
+# ## Error Analysis
+#
+# Compute the L2 error between numerical and exact solutions:
+
+dΩ = Measure(Ω,degree)
+eh = uh - u
+l2_uh = sqrt(sum(∫(eh⋅eh)*dΩ))
+
+writevtk(Ω,"results",cellfields=["uh"=>uh,"qh"=>qh,"eh"=>eh])
+
+# ## Going Further
+#
+# By modifying the stabilisation term, HDG can also work on polytopal meshes. An driver 
+# solving the same problem on a polytopal mesh is available [in the Gridap repository](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HDGPolytopalTests.jl). 
+# A tutorial for HHO on polytopal meshes is also available. 
+#

src/poisson_hho.jl

@@ -0,0 +1,199 @@
+# # [Mixed-order Hybrid High-Order Method for the Poisson equation](@id poisson_hho)
+#
+# In this tutorial, we will learn how to implement a mixed-order Hybrid High-Order (HHO) method
+# for solving the Poisson equation. HHO methods are a class of modern hybridizable finite element methods
+# that provide optimal convergence rates while enabling static condensation for efficient solution.
+#
+# ## Problem statement
+#
+# We consider the Poisson equation with Dirichlet boundary conditions:
+#
+# ```math
+# \begin{aligned}
+# -\Delta u &= f \quad \text{in } \Omega \\
+# u &= g \quad \text{on } \partial\Omega
+# \end{aligned}
+# ```
+#
+# where Ω is a bounded domain in R², f is a source term, and g is the prescribed boundary value.
+#
+# ## HHO discretization
+#
+# The HHO method introduces two types of unknowns:
+# 1. Cell unknowns defined in the volume of each mesh cell
+# 2. Face unknowns defined on the mesh faces/edges
+#
+# This hybrid structure allows for efficient static condensation by eliminating the cell unknowns
+# algebraically at the element level.
+#
+# Load the required packages
+
+using Gridap
+using Gridap.Geometry, Gridap.FESpaces, Gridap.MultiField
+using Gridap.CellData, Gridap.Fields, Gridap.Helpers
+using Gridap.ReferenceFEs
+using Gridap.Arrays
+
+# ## Local projection operator
+# 
+# Define a projection operator to map functions onto local polynomial spaces
+
+function projection_operator(V, Ω, dΩ)
+  Π(u,Ω) = change_domain(u,Ω,DomainStyle(u))
+  mass(u,v) = ∫(u⋅Π(v,Ω))dΩ
+  V0 = FESpaces.FESpaceWithoutBCs(V)
+  P = LocalOperator(
+    LocalSolveMap(), V0, mass, mass; trian_out = Ω
+  )
+  return P
+end
+
+# ## Reconstruction operator
+#
+# Define a reconstruction operator that maps hybrid unknowns to a higher-order polynomial space.
+# This operator is key for achieving optimal convergence rates.
+
+function reconstruction_operator(ptopo,order,X,Ω,Γp,dΩp,dΓp)
+  L = FESpaces.PolytopalFESpace(Ω, Float64, order+1; space=:P)
+  Λ = FESpaces.PolytopalFESpace(Ω, Float64, 0; space=:P)
+
+  n = get_normal_vector(Γp)
+  Πn(v) = ∇(v)⋅n
+  Π(u,Ω) = change_domain(u,Ω,DomainStyle(u))
+  lhs((u,λ),(v,μ))   = ∫( (∇(u)⋅∇(v)) + (μ*u) + (λ*v) )dΩp
+  rhs((uT,uF),(v,μ)) =  ∫( (∇(uT)⋅∇(v)) + (uT*μ) )dΩp + ∫( (uF - Π(uT,Γp))*(Πn(v)) )dΓp
+  
+  Y = FESpaces.FESpaceWithoutBCs(X)
+  W = MultiFieldFESpace([L,Λ];style=BlockMultiFieldStyle())
+  R = LocalOperator(
+    LocalPenaltySolveMap(), ptopo, W, Y, lhs, rhs; space_out = L
+  )
+  return R
+end
+
+# ## Problem setup
+#
+# Define the exact solution and forcing term
+
+u(x) = sin(2*π*x[1])*sin(2*π*x[2])
+f(x) = -Δ(u)(x)
+
+# Setup the mesh and discretization parameters
+
+n = 10
+base_model = simplexify(CartesianDiscreteModel((0,1,0,1),(n,n)))
+model = Geometry.voronoi(base_model)
+
+D = num_cell_dims(model)
+Ω = Triangulation(ReferenceFE{D}, model)
+Γ = Triangulation(ReferenceFE{D-1}, model)
+
+ptopo = Geometry.PatchTopology(model)
+Ωp = Geometry.PatchTriangulation(model,ptopo)
+Γp = Geometry.PatchBoundaryTriangulation(model,ptopo)
+
+order = 1
+qdegree = 2*(order+1)
+
+dΩp = Measure(Ωp,qdegree)
+dΓp = Measure(Γp,qdegree)
+
+# ## FE spaces and operators
+#
+# Define the finite element spaces for cell and face unknowns
+
+V = FESpaces.PolytopalFESpace(Ω, Float64, order+1; space=:P) # Bulk space
+M = FESpaces.PolytopalFESpace(Γ, Float64, order; space=:P, dirichlet_tags="boundary") # Skeleton space
+N = TrialFESpace(M,u)
+
+mfs = MultiField.BlockMultiFieldStyle(2,(1,1))
+X   = MultiFieldFESpace([V, N];style=mfs)
+Y   = MultiFieldFESpace([V, M];style=mfs) 
+Xp  = FESpaces.PatchFESpace(X,ptopo)
+
+# Setup projection and reconstruction operators
+
+PΓ = projection_operator(M, Γp, dΓp)
+R  = reconstruction_operator(ptopo,order,Y,Ωp,Γp,dΩp,dΓp)
+
+# Setup assemblers
+
+global_assem = SparseMatrixAssembler(X,Y)
+patch_assem = FESpaces.PatchAssembler(ptopo,X,Y)
+
+# ## Bilinear and linear forms
+#
+# Define the bilinear form a(u,v) for the diffusion term
+
+function a(u,v)
+  Ru_Ω, Ru_Γ = R(u)
+  Rv_Ω, Rv_Γ = R(v)
+  return ∫(∇(Ru_Ω)⋅∇(Rv_Ω) + ∇(Ru_Γ)⋅∇(Rv_Ω) + ∇(Ru_Ω)⋅∇(Rv_Γ) + ∇(Ru_Γ)⋅∇(Rv_Γ))dΩp
+end
+
+# Compute the inverse of local cell measure for stabilization
+
+hTinv = CellField(1 ./ collect(get_array(∫(1)dΩp)), Ωp)
+
+# Define the stabilization term s(u,v) to weakly enforce continuity
+
+function s(u,v)
+  function SΓ(u)
+    u_Ω, u_Γ = u
+    return PΓ(u_Ω) - u_Γ
+  end
+  return ∫(hTinv * (SΓ(u)⋅SΓ(v)))dΓp
+end
+
+# Define the linear form l(v) for the source term
+
+l((vΩ,vΓ)) = ∫(f⋅vΩ)dΩp
+
+# ## Problem solution
+#
+# Set up the weak form and solve using direct or static condensation
+
+function weakform()
+  u, v = get_trial_fe_basis(X), get_fe_basis(Y)
+  data = FESpaces.collect_and_merge_cell_matrix_and_vector(
+    (Xp, Xp, a(u,v), DomainContribution(), zero(Xp)),
+    (X, Y, s(u,v), l(v), zero(X))
+  )
+  assemble_matrix_and_vector(global_assem,data)
+end
+
+function patch_weakform()
+  u, v = get_trial_fe_basis(X), get_fe_basis(Y)
+  data = FESpaces.collect_and_merge_cell_matrix_and_vector(patch_assem,
+    (Xp, Xp, a(u,v), DomainContribution(), zero(Xp)),
+    (X, Y, s(u,v), l(v), zero(X))
+  )
+  return assemble_matrix_and_vector(patch_assem,data)
+end
+
+# Direct monolithic solve
+
+A, b = weakform()
+x = A \ b
+
+ui, ub = FEFunction(X,x)
+eu  = ui - u 
+l2u = sqrt(sum( ∫(eu * eu)dΩp))
+h1u = l2u + sqrt(sum( ∫(∇(eu) ⋅ ∇(eu))dΩp))
+
+# Static condensation
+
+op = MultiField.StaticCondensationOperator(X,patch_assem,patch_weakform())
+ui, ub = solve(op) 
+
+eu  = ui - u
+l2u = sqrt(sum( ∫(eu * eu)dΩp))
+h1u = l2u + sqrt(sum( ∫(∇(eu) ⋅ ∇(eu))dΩp))
+
+# The code above demonstrates both solution approaches:
+#
+# 1. Direct solution of the full system
+# 2. Static condensation to eliminate cell unknowns
+#
+# Both give the same solution but static condensation is typically more efficient
+# for higher orders.