Improved new tutorials

Author: JordiManyer

Project.toml

@@ -41,7 +41,7 @@ GridapDistributed = "0.4"
 GridapGmsh = "0.7"
 GridapP4est = "0.3"
 GridapPETSc = "0.5"
-GridapSolvers = "0.5.0"
+GridapSolvers = "0.6"
 MPI = "0.20"
 SpecialFunctions = "1"
 julia = "1.6"

deps/build.jl

@@ -26,6 +26,9 @@ files = [
   "Topology optimization"=>"TopOptEMFocus.jl",
   "Poisson on unfitted meshes"=>"poisson_unfitted.jl",
   "Poisson with AMR"=>"poisson_amr.jl",
+  "Poisson with HDG"=>"poisson_hdg.jl",
+  "Poisson with HHO on polytopal meshes"=>"poisson_hho.jl",
+  "Block assembly and solvers: Incompressible Stokes example"=>"stokes_blocks.jl",
   "Lagrange multipliers" => "lagrange_multipliers.jl",
   "Low-level API - Poisson equation"=>"poisson_dev_fe.jl",
   "Low-level API - Geometry" => "geometry_dev.jl",

docs/make.jl

@@ -60,15 +60,15 @@ for (i,(title,filename)) in enumerate(Tutorials.files)
 end
 
 makedocs(
-    sitename = "Gridap tutorials",
-    format = Documenter.HTML(),
-    pages = pages
+  sitename = "Gridap tutorials",
+  format = Documenter.HTML(),
+  pages = pages
 )
 
 # Documenter can also automatically deploy documentation to gh-pages.
 # See "Hosting Documentation" and deploydocs() in the Documenter manual
 # for more information.
 
 deploydocs(
-    repo = "github.com/gridap/Tutorials.git"
+  repo = "github.com/gridap/Tutorials.git"
 )

src/poisson_hdg.jl

@@ -1,10 +1,8 @@
-# # [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.
+# that introduces an auxiliary variable m on mesh interfaces to reduce the global system size.
 #
-# ## The Poisson Problem
+# ## HDG Discretization
 #
 # We consider the Poisson equation with Dirichlet boundary conditions:
 #
@@ -15,29 +13,28 @@
 # \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\\
+# \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
+# - ``m_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}``
+# Numerical fluxes are defindes as 
 #
-# where τ is a stabilization parameter.
+# ```math
+# \widehat{\boldsymbol{q}}_h = \boldsymbol{q}_h + \tau(u_h - m_h)\boldsymbol{n}
+# ```
+#
+# where $\tau$ is a stabilization parameter.
 #
 # First, let's load the required Gridap packages:
 
@@ -54,108 +51,105 @@ using Gridap.CellData
 
 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)
+f(x) = (∇ ⋅ q)(x) # Source term f = -Δu = -∇⋅(∇u)$
 
-# ## Mesh Generation
+# ## Geometry
 #
-# Create a 3D simplicial mesh from a Cartesian grid:
+# We generate a D-dimensional 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
+# From this mesh, we will require two triangulations where to define our HDG spaces:
 #
-# HDG methods require two types of meshes:
-# 1. Volume mesh for element interiors
-# 2. Skeleton mesh for element interfaces
+# 1. A cell triangulation $\Omega$, for the volume variables
+# 2. A face triangulation $\Gamma$, for the skeleton variables
 #
-# We also need patch-wise triangulations for local computations:
+# These are given by
 
 Ω = Triangulation(ReferenceFE{D}, model)  # Volume triangulation
 Γ = Triangulation(ReferenceFE{D-1}, model)  # Skeleton triangulation
 
+# ## PatchTopology and PatchTriangulation
+#
+# A key aspect of hybrid methods is the use of static condensation, which is the 
+# elimination of cell unknowns to reduce the size of the global system.
+# To achieve this, we need to be able to assemble and solve local problems on each cell, that 
+# involve 
+#   - contributions from the cell itself
+#   - contributions from the cell faces
+# To this end, Gridap provides a general framework for patch-assembly and solves. The idea 
+# is to define a patch decomposition of the mesh (in this case, a patch is a cell and its sourrounding 
+# faces). We can then gather contributions for each patch, solve the local problems, and 
+# assemble the results into the global system.
+#
+# The following code creates the required `PatchTopology` for the problem at hand. We then 
+# take d-dimensional slices of it by the means of `PatchTriangulation` and `PatchBoundaryTriangulation`.
+# These are the `Triangulation`s we will integrate our weakform over.
+
 ptopo = Geometry.PatchTopology(model)
 Ωp = Geometry.PatchTriangulation(model,ptopo)  # Patch volume triangulation
 Γp = Geometry.PatchBoundaryTriangulation(model,ptopo)  # Patch skeleton triangulation
 
-# ## Reference Finite Elements
+# ## FESpaces
 #
 # 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)
+# 1. A vector-valued space for the flux q (Q)
+# 2. A scalar space for the solution u (V)
+# 3. A scalar space for the interface variable m (M)
+#
+# We then define discontinuous finite element spaces of the approppriate order, locally $\mathbb{P}^k$.
+# Note that only the skeletal space has Dirichlet boundary conditions.
 
 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:
+# Since we are doing static condensation, we need assemble by blocks. In particular, the 
+# `StaticCondensationOperator` expects the variables to be groupped in two blocks:
+#   - The eliminated variables (in this case, the volume variables q and u)
+#   - The retained variables (in this case, the interface variable m)
+# We will assemble by blocks using the `BlockMultiFieldStyle` API. 
 
-mfs = MultiField.BlockMultiFieldStyle(2,(2,1))  # Special blocking for efficient static condensation
-X = MultiFieldFESpace([V, Q, N];style=mfs)
+mfs = BlockMultiFieldStyle(2,(2,1))  # Special blocking for efficient static condensation
+X = MultiFieldFESpace([V, Q, N]; style=mfs)
 
-# ## Integration Setup
+# ## Weak Form and integration
 #
-# 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:
+dΩp = Measure(Ωp,degree)  # Volume measure, on the patch triangulation
+dΓp = Measure(Γp,degree)  # Surface measure, on the patch boundary triangulation
 
 τ = 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:
+# With all these ingredients, we can now build our statically-condensed operator ans 
+# solve the problem. Note that we are solving a scatically-condensed system. We can 
+# retrieve the internal `AffineFEOperator` from `op.sc_op`.
 
 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Ω))