From 9db6a5afc2609def684c0ad3aecde12ec6d2955d Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Fri, 13 Jun 2025 00:11:03 +1000 Subject: [PATCH 01/10] Added new tutorials --- Project.toml | 4 +- src/poisson_hdg.jl | 170 +++++++++++++++++++++++++++++++++++ src/poisson_hho.jl | 199 +++++++++++++++++++++++++++++++++++++++++ src/stokes_blocks.jl | 207 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 579 insertions(+), 1 deletion(-) create mode 100644 src/poisson_hdg.jl create mode 100644 src/poisson_hho.jl create mode 100644 src/stokes_blocks.jl diff --git a/Project.toml b/Project.toml index 506480e..b5b01e3 100644 --- a/Project.toml +++ b/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" diff --git a/src/poisson_hdg.jl b/src/poisson_hdg.jl new file mode 100644 index 0000000..937c712 --- /dev/null +++ b/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. +# diff --git a/src/poisson_hho.jl b/src/poisson_hho.jl new file mode 100644 index 0000000..6f468c2 --- /dev/null +++ b/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. diff --git a/src/stokes_blocks.jl b/src/stokes_blocks.jl new file mode 100644 index 0000000..f434673 --- /dev/null +++ b/src/stokes_blocks.jl @@ -0,0 +1,207 @@ +# # Incompressible Stokes equations in a 2D/3D cavity +# +# This example solves the incompressible Stokes equations, given by +# +# ```math +# \begin{align*} +# -\Delta u - \nabla p &= f \quad \text{in} \quad \Omega, \\ +# \nabla \cdot u &= 0 \quad \text{in} \quad \Omega, \\ +# u &= \hat{x} \quad \text{in} \quad \Gamma_\text{top} \subset \partial \Omega, \\ +# u &= 0 \quad \text{in} \quad \partial \Omega \backslash \Gamma_\text{top} \\ +# \end{align*} +# ``` +# +# where $\Omega = [0,1]^d$. +# +# We use a mixed finite-element scheme, with $Q_k \times P_{k-1}^{-}$ elements for the velocity-pressure pair. +# +# To solve the linear system, we use a FGMRES solver preconditioned by a block-diagonal or +# block-triangular Shur-complement-based preconditioner. +# +# +# ## Block structure +# +# The discretized system has a natural 2×2 block structure: +# +# ```math +# \begin{bmatrix} +# A & B^T \\ +# B & 0 +# \end{bmatrix} +# \begin{bmatrix} +# u \\ +# p +# \end{bmatrix} = +# \begin{bmatrix} +# f \\ +# 0 +# \end{bmatrix} +# ``` +# +# where: +# - A: Vector Laplacian (velocity block) +# - B: Divergence operator +# - B^T: Gradient operator +# +# ## Solution strategy +# +# We use a FGMRES solver preconditioned by an upper block-triangular preconditioner: +# +# ```math +# P = \begin{bmatrix} +# A & B^T \\ +# 0 & -\hat{S} +# \end{bmatrix} +# ``` +# +# where $\hat{S}$ is an approximation of the Schur complement $S = BA^{-1}B^T$, which we +# will approximate using a pressure mass matrix. + +using LinearAlgebra +using BlockArrays + +using Gridap +using Gridap.MultiField + +using GridapSolvers +using GridapSolvers.LinearSolvers +using GridapSolvers.BlockSolvers: LinearSystemBlock, BiformBlock +using GridapSolvers.BlockSolvers: BlockDiagonalSolver, BlockTriangularSolver + +# ## Geometry and FESpaces +# +# See the basic Stokes tutorial for a detailed explanation. +# + +nc = (8,8) +Dc = length(nc) +domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) + +model = CartesianDiscreteModel(domain,nc) +labels = get_face_labeling(model) +if Dc == 2 + add_tag_from_tags!(labels,"top",[6]) + add_tag_from_tags!(labels,"walls",[1,2,3,4,5,7,8]) +else + add_tag_from_tags!(labels,"top",[22]) + add_tag_from_tags!(labels,"walls",[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25,26]) +end + +order = 2 +qdegree = 2*(order+1) +reffe_u = ReferenceFE(lagrangian,VectorValue{Dc,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +u_walls = (Dc==2) ? VectorValue(0.0,0.0) : VectorValue(0.0,0.0,0.0) +u_top = (Dc==2) ? VectorValue(1.0,0.0) : VectorValue(1.0,0.0,0.0) + +V = TestFESpace(model,reffe_u,dirichlet_tags=["walls","top"]); +U = TrialFESpace(V,[u_walls,u_top]); +Q = TestFESpace(model,reffe_p;conformity=:L2,constraint=:zeromean) + +# Our first difference will come from how we define our multi-field spaces: +# Because we want to be able to use the block-structure of the linear system, +# we have to assemble our problem by blocks. The block structure of the resulting +# linear system is determined by the `BlockMultiFieldStyle` we use to define the +# multi-field spaces. +# +# A `BlockMultiFieldStyle` takes three arguments: +# - `N`: The number of blocks we want +# - `S`: An N-tuple with the number of fields in each block +# - `P`: A permutation of the fields in the multi-field space, which determines +# how fields are grouped into blocks. +# +# By default, we create as many blocks as there are fields in the multi-field space, +# and each block contains a single field with no permutation. + +mfs = BlockMultiFieldStyle(2,(1,1),(1,2)) +X = MultiFieldFESpace([U,Q];style=mfs) +Y = MultiFieldFESpace([V,Q];style=mfs) + +# ## Weak form + +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) + +α = 1.e1 +f = (Dc==2) ? VectorValue(1.0,1.0) : VectorValue(1.0,1.0,1.0) +a((u,p),(v,q)) = ∫(∇(v)⊙∇(u))dΩ - ∫(divergence(v)*p)dΩ - ∫(divergence(u)*q)dΩ +l((v,q)) = ∫(v⋅f)dΩ + +op = AffineFEOperator(a,l,X,Y) + +# ## Block structure of the linear system +# +# As per usual, we can extract the matrix and vector of the linear system from the operator. +# Notice now that unlike in previous examples, the matrix is a `BlockMatrix` type, from +# the `BlockArrays.jl` package, which allows us to work with block-structured matrices. + +A = get_matrix(op) +b = get_vector(op) + +# ## Block solvers +# +# We will now setup two types of block preconditioners for the Stokes system. In both cases, +# we will use the preconditioners to solve the linear system using a Flexible GMRES solver. +# The idea behind these preconditioners is the well-known property that the Schur complement +# of the velocity block can be well approximated by a scaled pressure mass matrix. Moreover, +# in our pressure discretization is discontinuous which means that the pressure mass matrix +# is block-diagonal and easily invertible. +# In this example, we will use an exact LU solver for the velocity block and a CG solver +# with Jacobi preconditioner for the pressure block. + +# ### Block diagonal preconditioner +# +# The simplest block preconditioner is the block-diagonal preconditioner. +# The only ingredients required are +# - the sub-solvers for each diagonal block and +# - the diagonal blocks we want to use. +# +# The sub-solvers are defined as follows: + +u_solver = LUSolver() +p_solver = CGSolver(JacobiLinearSolver();maxiter=20,atol=1e-14,rtol=1.e-6) + +# The block structure is defined using the block API. We provide different types of blocks, +# that might have different uses depending on the problem at hand. We will here use two of +# the most common block types: +# - The `LinearSystemBlock` defines a block that is taken directly (and aliased) form +# the linear system matrix `A`. We will use this for the velocity block. +# - For the pressure block, however, the pressure mass matrix is not directly available +# in the system matrix. Instead, we will have to integrate it using the Gridap API, as usual. +# These abstract concept is implemented in the `BiformBlock` type, which allows the +# user to define a block from a bilinear form. +# All in all, we define the block structure as follows: + +u_block = LinearSystemBlock() +p_block = BiformBlock((p,q) -> ∫(-(1.0/α)*p*q)dΩ,Q,Q) + +# With these ingredients, we can now define the block diagonal preconditioner as follows: + +PD = BlockDiagonalSolver([u_block,p_block],[u_solver,p_solver]) +solver_PD = FGMRESSolver(20,PD;atol=1e-10,rtol=1.e-12,verbose=true) + +uh, ph = solve(solver_PD, op) + +# ### Block upper-triangular preconditioner +# +# A slighly more elaborate preconditioner (but also more robust) is the +# block upper-triangular preconditioner. The ingredients are similar: +# - the sub-solvers for each diagonal block +# - the blocks that define the block structure, now including the off-diagonal blocks +# - the coefficients for the off-diagonal blocks, where zero coefficients +# indicate that the block is not used. +# +# We will also represent the off-diagonal blocks using the `LinearSystemBlock` type. + +bblocks = [ u_block LinearSystemBlock(); + LinearSystemBlock() p_block ] +coeffs = [1.0 1.0; + 0.0 1.0] +PU = BlockTriangularSolver(bblocks,[u_solver,p_solver],coeffs,:upper) +solver_PU = FGMRESSolver(20,PU;atol=1e-10,rtol=1.e-12,verbose=true) + +uh, ph = solve(solver_PU, op) + +# As you can see, the block upper-triangular preconditioner is quite better +# than the block diagonal one. From 642c9692f28246ab6174a99467b808df78fbea7c Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 10:53:35 +1000 Subject: [PATCH 02/10] Improved new tutorials --- Project.toml | 2 +- deps/build.jl | 3 + docs/make.jl | 8 +- src/poisson_hdg.jl | 114 +++++++++++------------- src/poisson_hho.jl | 208 +++++++++++++++++++++++++------------------ src/stokes_blocks.jl | 31 +++++-- 6 files changed, 204 insertions(+), 162 deletions(-) diff --git a/Project.toml b/Project.toml index b5b01e3..b51c713 100644 --- a/Project.toml +++ b/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" diff --git a/deps/build.jl b/deps/build.jl index 777cbf8..3032ec4 100644 --- a/deps/build.jl +++ b/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", diff --git a/docs/make.jl b/docs/make.jl index cf642a0..8bd6553 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -60,9 +60,9 @@ 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. @@ -70,5 +70,5 @@ makedocs( # for more information. deploydocs( - repo = "github.com/gridap/Tutorials.git" + repo = "github.com/gridap/Tutorials.git" ) diff --git a/src/poisson_hdg.jl b/src/poisson_hdg.jl index 937c712..aa2b9b5 100644 --- a/src/poisson_hdg.jl +++ b/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,14 +13,12 @@ # \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} # ``` @@ -30,14 +26,15 @@ # 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Ω)) diff --git a/src/poisson_hho.jl b/src/poisson_hho.jl index 6f468c2..943e8dd 100644 --- a/src/poisson_hho.jl +++ b/src/poisson_hho.jl @@ -1,4 +1,3 @@ -# # [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 @@ -15,18 +14,18 @@ # \end{aligned} # ``` # -# where Ω is a bounded domain in R², f is a source term, and g is the prescribed boundary value. +# where $\Omega$ is a bounded domain in $\mathbb{R}^2$, 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 +# 2. Face unknowns defined on the mesh facets # # This hybrid structure allows for efficient static condensation by eliminating the cell unknowns # algebraically at the element level. # -# Load the required packages +# We start by loading the required packages using Gridap using Gridap.Geometry, Gridap.FESpaces, Gridap.MultiField @@ -34,9 +33,94 @@ using Gridap.CellData, Gridap.Fields, Gridap.Helpers using Gridap.ReferenceFEs using Gridap.Arrays -# ## Local projection operator +# ## Geometry +# +# We generate a 2-dimensional simplicial mesh from a Cartesian grid: + +u(x) = sin(2*π*x[1])*sin(2*π*x[2]) +f(x) = -Δ(u)(x) + +n = 10 +base_model = simplexify(CartesianDiscreteModel((0,1,0,1),(n,n))) +model = Geometry.voronoi(base_model) + +# From this mesh, we will require two triangulations where to define our HDG spaces: +# +# 1. A cell triangulation $\Omega$, for the volume variables +# 2. A face triangulation $\Gamma$, for the skeleton variables +# +# These are given by + +D = num_cell_dims(model) +Ω = Triangulation(ReferenceFE{D}, model) +Γ = Triangulation(ReferenceFE{D-1}, model) + +# ## FESpaces +# +# HHO uses two different finite element spaces: +# 1. A scalar space for the bulk variable uT (V) +# 2. A scalar space for the skeleton variable uF (M) +# +# We then define discontinuous finite element spaces of the approppriate order, locally $$\mathbb{P}^k$$. +# Because we are using a mixed-order scheme, the bulk space has a higher polynomial order +# than the skeleton space. +# Note that only the skeletal space has Dirichlet boundary conditions. + +order = 1 +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) + +# ## MultiField Structure +# +# 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,(1,1)) +X = MultiFieldFESpace([V, N];style=mfs) +Y = MultiFieldFESpace([V, M];style=mfs) +Xp = FESpaces.PatchFESpace(X,ptopo) + +# ## 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) +Γp = Geometry.PatchBoundaryTriangulation(model,ptopo) + +qdegree = 2*(order+1) +dΩp = Measure(Ωp,qdegree) +dΓp = Measure(Γp,qdegree) + +# ## Local operators # -# Define a projection operator to map functions onto local polynomial spaces +# A key feature of HHO is the use of local solves to define local projections of our bulk and +# skeleton variables. Just like for static condensation, we will use patch assembly to +# gather the contributions from each patch and solve the local problems. +# +# For the mixed-order Poisson problem, we require two local projections: +# - First, an L2 local projection operator onto the mesh faces. +# - Second, the so-called reconstruction operator. This operator is highly tied to the +# ellipic projector, and projects our bulk-skeleton variable pair onto a bulk +# space of higher order. +# The operators are defined as follows: function projection_operator(V, Ω, dΩ) Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) @@ -48,11 +132,6 @@ function projection_operator(V, Ω, dΩ) 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) @@ -71,59 +150,17 @@ function reconstruction_operator(ptopo,order,X,Ω,Γp,dΩp,dΓp) 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 +# ## Weakform # -# Define the bilinear form a(u,v) for the diffusion term +# We can now define: +# - The consistency term `a` +# - The stabilization term `s` +# - The rhs term `l` + +hTinv = CellField(1 ./ collect(get_array(∫(1)dΩp)), Ωp) function a(u,v) Ru_Ω, Ru_Γ = R(u) @@ -131,12 +168,6 @@ function a(u,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 @@ -145,13 +176,11 @@ function s(u,v) 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 +# ## Assembly without static condensation + +global_assem = SparseMatrixAssembler(X,Y) function weakform() u, v = get_trial_fe_basis(X), get_fe_basis(Y) @@ -162,17 +191,6 @@ function weakform() 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 @@ -181,7 +199,18 @@ eu = ui - u l2u = sqrt(sum( ∫(eu * eu)dΩp)) h1u = l2u + sqrt(sum( ∫(∇(eu) ⋅ ∇(eu))dΩp)) -# Static condensation +# ## Assembly with static condensation + +patch_assem = FESpaces.PatchAssembler(ptopo,X,Y) + +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 op = MultiField.StaticCondensationOperator(X,patch_assem,patch_weakform()) ui, ub = solve(op) @@ -190,10 +219,13 @@ 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: +# ## Going further +# +# This tutorial has introduced the basic concepts of HHO methods using the simplest their +# simplest form, e.g. mixed-order HHO for the Poisson equation. +# More advanced drivers can be found with Gridap's tests. In particular: # -# 1. Direct solution of the full system -# 2. Static condensation to eliminate cell unknowns +# - [Poisson with original HHO formulation](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HHOPolytopalTests.jl) +# - [Incompressible Stokes](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HHOMixedStokesPolytopal.jl) +# - [Linear Elasticity](https://github.com/gridap/Gridap.jl/blob/75efc9a7a7e286c27e7ca3ddef5468e591845484/test/GridapTests/HHOMixedElasticity.jl) # -# Both give the same solution but static condensation is typically more efficient -# for higher orders. diff --git a/src/stokes_blocks.jl b/src/stokes_blocks.jl index f434673..8523e03 100644 --- a/src/stokes_blocks.jl +++ b/src/stokes_blocks.jl @@ -3,12 +3,14 @@ # This example solves the incompressible Stokes equations, given by # # ```math -# \begin{align*} +# \left\lbrace +# \begin{aligned} # -\Delta u - \nabla p &= f \quad \text{in} \quad \Omega, \\ # \nabla \cdot u &= 0 \quad \text{in} \quad \Omega, \\ # u &= \hat{x} \quad \text{in} \quad \Gamma_\text{top} \subset \partial \Omega, \\ # u &= 0 \quad \text{in} \quad \partial \Omega \backslash \Gamma_\text{top} \\ -# \end{align*} +# \end{aligned} +# \right. # ``` # # where $\Omega = [0,1]^d$. @@ -39,9 +41,9 @@ # ``` # # where: -# - A: Vector Laplacian (velocity block) -# - B: Divergence operator -# - B^T: Gradient operator +# - $$A$$: Vector Laplacian (velocity block) +# - $$B$$: Divergence operator +# - $$B^T$$: Gradient operator # # ## Solution strategy # @@ -99,6 +101,8 @@ V = TestFESpace(model,reffe_u,dirichlet_tags=["walls","top"]); U = TrialFESpace(V,[u_walls,u_top]); Q = TestFESpace(model,reffe_p;conformity=:L2,constraint=:zeromean) +# ## Block multi-field spaces +# # Our first difference will come from how we define our multi-field spaces: # Because we want to be able to use the block-structure of the linear system, # we have to assemble our problem by blocks. The block structure of the resulting @@ -118,7 +122,7 @@ mfs = BlockMultiFieldStyle(2,(1,1),(1,2)) X = MultiFieldFESpace([U,Q];style=mfs) Y = MultiFieldFESpace([V,Q];style=mfs) -# ## Weak form +# ## Weak form and integration Ω = Triangulation(model) dΩ = Measure(Ω,qdegree) @@ -130,7 +134,7 @@ l((v,q)) = ∫(v⋅f)dΩ op = AffineFEOperator(a,l,X,Y) -# ## Block structure of the linear system +# ### Block structure of the linear system # # As per usual, we can extract the matrix and vector of the linear system from the operator. # Notice now that unlike in previous examples, the matrix is a `BlockMatrix` type, from @@ -194,14 +198,23 @@ uh, ph = solve(solver_PD, op) # # We will also represent the off-diagonal blocks using the `LinearSystemBlock` type. -bblocks = [ u_block LinearSystemBlock(); +sblocks = [ u_block LinearSystemBlock(); LinearSystemBlock() p_block ] coeffs = [1.0 1.0; 0.0 1.0] -PU = BlockTriangularSolver(bblocks,[u_solver,p_solver],coeffs,:upper) +PU = BlockTriangularSolver(sblocks,[u_solver,p_solver],coeffs,:upper) solver_PU = FGMRESSolver(20,PU;atol=1e-10,rtol=1.e-12,verbose=true) uh, ph = solve(solver_PU, op) # As you can see, the block upper-triangular preconditioner is quite better # than the block diagonal one. + +# ## Going further +# +# If you want to see more examples of how to use the block solvers, +# you can check the documentation in [GridapSolvers.jl](https://gridap.github.io/GridapSolvers.jl/stable/), +# as well as it's `test/Applications` folder. +# +# There you will find more complicated examples, such as using a GMG solver to +# solve the velocity block. From 2edeb228140096666e9fc270c6615ce2b1ce9c78 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 11:32:11 +1000 Subject: [PATCH 03/10] Minor --- src/poisson_hho.jl | 70 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 60 insertions(+), 10 deletions(-) diff --git a/src/poisson_hho.jl b/src/poisson_hho.jl index 943e8dd..43186b6 100644 --- a/src/poisson_hho.jl +++ b/src/poisson_hho.jl @@ -3,6 +3,12 @@ # 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. # +# HHO is a mathematically complex method. This tutorial will **NOT** cover the method nor +# its mathematical foundations. You should be familiar with both before going into +# the tutorial itself. +# We also recommend going through the HDG tutorial first, which shares many of the +# same concepts but is simpler to understand. +# # ## Problem statement # # We consider the Poisson equation with Dirichlet boundary conditions: @@ -82,7 +88,6 @@ 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) # ## PatchTopology and PatchTriangulation # @@ -114,13 +119,26 @@ dΓp = Measure(Γp,qdegree) # A key feature of HHO is the use of local solves to define local projections of our bulk and # skeleton variables. Just like for static condensation, we will use patch assembly to # gather the contributions from each patch and solve the local problems. +# The result is then an element of the space we are projecting to, given as a +# linear combination of the basis of that space. This whole abstraction is taken +# care of by the `LocalOperator` object. +# +# For the mixed-order Poisson problem, we require two local projections. +# +# ### L2 projection operator # -# For the mixed-order Poisson problem, we require two local projections: -# - First, an L2 local projection operator onto the mesh faces. -# - Second, the so-called reconstruction operator. This operator is highly tied to the -# ellipic projector, and projects our bulk-skeleton variable pair onto a bulk -# space of higher order. -# The operators are defined as follows: +# First, an L2 local projection operator onto the mesh faces. This is the simplest +# operator we can define, such that given $u$ in some undefined space, we find $Pu \in V$ st +# +# $$(Pu,q) = (u,q) \forall q \in V$$ +# +# This signature for the `LocalOperator` assumes that the rhs and lhs for each local problem +# are given by a single cell contribution (no patch assembly required). +# Note also the use of `FESpaceWithoutBCs`, which strips the boundary conditions from the +# space `V`. This is because we do not want to take into account boundary conditions +# when projecting onto the space. +# The local solve map is given by `LocalSolveMap`, which by default uses Julia's LU +# factorization to solve the local problem exactly. function projection_operator(V, Ω, dΩ) Π(u,Ω) = change_domain(u,Ω,DomainStyle(u)) @@ -132,6 +150,18 @@ function projection_operator(V, Ω, dΩ) return P end +# ### Reconstruction operator +# +# Finally, we build the so-called reconstruction operator. This operator is highly tied to the +# ellipic projector, and projects our bulk-skeleton variable pair onto a bulk space of higher order. +# +# It's definition can be found in HHO literature, and requires solving a constrained +# local problem on each cell and its faces. We therefore use patch assembly, and impose +# our constraint using an additional space `Λ` as a local Lagrange multiplier. +# Naturally, we want to eliminate the multiplier and return a solution in the reconstructed +# space `L`. This is taken care of by the `LocalPenaltySolveMap`, and the kwarg `space_out = L` +# which overrides the default behavior of returning a solution in the test space `Y`. + 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) @@ -153,7 +183,7 @@ end PΓ = projection_operator(M, Γp, dΓp) R = reconstruction_operator(ptopo,order,Y,Ωp,Γp,dΩp,dΓp) -# ## Weakform +# ## Weakform and assembly # # We can now define: # - The consistency term `a` @@ -178,7 +208,27 @@ end l((vΩ,vΓ)) = ∫(f⋅vΩ)dΩp -# ## Assembly without static condensation +# ### Patch-FESpaces +# +# An additional difficulty in HHO methods is that our reconstructed functions `R(u)` are +# hard to assemble. They are defined on the cells, but depend on skeleton degrees of freedom. +# We therefore cannot assemble them using the original FESpace `N`. Instead, we will create \ +# view of the original space that is defined on the patches, and can be used to assemble +# the local contributions depending on `R`. It can be done by the means of `PatchFESpace`: + +Xp = FESpaces.PatchFESpace(X,ptopo) + +# ### Assembly without static condensation +# +# We can now proceed to evaluate and assemble all our contributions. Note that some of +# the contributions depend on the reconstructed variables, e.g `a(u,v)`, and need to +# be assembled using the `PatchFESpace` `Xp`, while others can be assembled using the original +# FESpace `X` (e.g. `s(u,v)` and `l(v)`). +# Assembly is therefore somewhat more complex than in the standard case. We need to +# collect the contributions for every (test,trial) pair, and then merge them into a single +# data structure that can be passed to the assembler. +# In the case of mixed-order HHO, we only have two combinations, (`Xp`, `Xp`) and (`X`, `X`), +# but the cross-terms appear also in the case of the original HHO formulation. global_assem = SparseMatrixAssembler(X,Y) @@ -199,7 +249,7 @@ eu = ui - u l2u = sqrt(sum( ∫(eu * eu)dΩp)) h1u = l2u + sqrt(sum( ∫(∇(eu) ⋅ ∇(eu))dΩp)) -# ## Assembly with static condensation +# ### Assembly with static condensation patch_assem = FESpaces.PatchAssembler(ptopo,X,Y) From d8ae05d2c668998c2b3dd0d00547c6abd7824ae2 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 11:50:46 +1000 Subject: [PATCH 04/10] Better CI --- .github/workflows/ci.yml | 107 ++++++++------------------------------- 1 file changed, 20 insertions(+), 87 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a14c19f..e3fcb2c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -5,8 +5,8 @@ on: - create jobs: - tutorials_set_1: - name: Tutorials1 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + tutorials: + name: Tutorials ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - ${{ matrix.test_set.name }} runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -17,6 +17,23 @@ jobs: - ubuntu-latest arch: - x64 + test_set: + - name: Basics + files: "poisson.jl elasticity.jl lagrange_multipliers.jl dg_discretization.jl darcy.jl stokes.jl advection_diffusion.jl" + - name: Nonlinear + files: "p_laplacian.jl hyperelasticity.jl inc_navier_stokes.jl" + - name: Applications + files: "fsi_tutorial.jl emscatter.jl isotropic_damage.jl TopOptEMFocus.jl" + - name: Transient + files: "isotropic_damage.jl transient_linear.jl transient_nonlinear.jl " + - name: Hybrid + files: "poisson_hdg.jl poisson_hho.jl" + - name: Plugins + files: "stokes_blocks.jl poisson_amr.jl poisson_unfitted.jl" + - name: MPI + files: "runtests_mpi.jl" + - name: Other + files: "validation.jl validation_DrWatson.jl interpolation_fe.jl poisson_dev_fe.jl geometry_dev.jl" steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 @@ -25,91 +42,7 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl poisson.jl validation.jl elasticity.jl lagrange_multipliers.jl - tutorials_set_2: - name: Tutorials2 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl p_laplacian.jl hyperelasticity.jl dg_discretization.jl fsi_tutorial.jl poisson_amr.jl - tutorials_set_3: - name: Tutorials3 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl darcy.jl inc_navier_stokes.jl stokes.jl poisson_dev_fe.jl emscatter.jl - tutorials_set_4: - name: Tutorials4 ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl isotropic_damage.jl validation_DrWatson.jl interpolation_fe.jl transient_linear.jl transient_nonlinear.jl TopOptEMFocus.jl - tutorials_mpi: - name: TutorialsMPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - version: - - '1.10' - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests_mpi.jl + - run: cd test/; julia --project=.. --color=yes --check-bounds=yes ${{ matrix.test_set.files }} docs: name: Documentation runs-on: ubuntu-latest From aa74feb7dc422771013afa819bd94e084eb9f121 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 12:02:23 +1000 Subject: [PATCH 05/10] Minor --- .github/workflows/ci.yml | 5 ++++- test/maketests.jl | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 1 deletion(-) create mode 100644 test/maketests.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e3fcb2c..0c48555 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,7 +42,10 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes ${{ matrix.test_set.files }} + - run: | + cd test/ + julia --project=.. maketests.jl + julia --project=.. --color=yes --check-bounds=yes ${{ matrix.test_set.files }} docs: name: Documentation runs-on: ubuntu-latest diff --git a/test/maketests.jl b/test/maketests.jl new file mode 100644 index 0000000..7516618 --- /dev/null +++ b/test/maketests.jl @@ -0,0 +1,33 @@ +using Tutorials +using Test + +function get_title(files,filename) + for (title,fn) in files + if fn == filename + return title + end + end + error("File $filename not found!") +end + +if (length(ARGS) != 0) + files = [get_title(Tutorials.files,filename)=>filename for filename in ARGS] +else + files = Tutorials.files +end + +for (title,filename) in files + # Create temporal modules to isolate and protect test scopes + tmpdir = mktempdir(;cleanup=true) + tmpmod = split(filename,".")[1] + tmpfile = joinpath(tmpdir,tmpmod) + isfile(tmpfile) && error("File $tmpfile already exists!") + testpath = joinpath(@__DIR__,"../src", filename) + open(tmpfile,"w") do f + println(f, "# This file is automatically generated") + println(f, "# Do not edit") + println(f) + println(f, "module $tmpmod include(\"$testpath\") end") + end + # @time @testset "$title" begin include(tmpfile) end +end From 0fa2bd582763f32495191434e56a3532c6a9462b Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 12:15:40 +1000 Subject: [PATCH 06/10] Minor --- .github/workflows/ci.yml | 27 ++++++++++++++++++++++++--- test/maketests.jl | 33 --------------------------------- test/runtests.jl | 2 +- 3 files changed, 25 insertions(+), 37 deletions(-) delete mode 100644 test/maketests.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0c48555..9f19435 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -30,8 +30,6 @@ jobs: files: "poisson_hdg.jl poisson_hho.jl" - name: Plugins files: "stokes_blocks.jl poisson_amr.jl poisson_unfitted.jl" - - name: MPI - files: "runtests_mpi.jl" - name: Other files: "validation.jl validation_DrWatson.jl interpolation_fe.jl poisson_dev_fe.jl geometry_dev.jl" steps: @@ -45,7 +43,30 @@ jobs: - run: | cd test/ julia --project=.. maketests.jl - julia --project=.. --color=yes --check-bounds=yes ${{ matrix.test_set.files }} + julia --project=.. --color=yes --check-bounds=yes runtests.jl ${{ matrix.test_set.files }} + tutorials_mpi: + name: Tutorials MPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.10' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - run: | + cd test/ + julia --project=.. --color=yes --check-bounds=yes runtests_mpi.jl docs: name: Documentation runs-on: ubuntu-latest diff --git a/test/maketests.jl b/test/maketests.jl deleted file mode 100644 index 7516618..0000000 --- a/test/maketests.jl +++ /dev/null @@ -1,33 +0,0 @@ -using Tutorials -using Test - -function get_title(files,filename) - for (title,fn) in files - if fn == filename - return title - end - end - error("File $filename not found!") -end - -if (length(ARGS) != 0) - files = [get_title(Tutorials.files,filename)=>filename for filename in ARGS] -else - files = Tutorials.files -end - -for (title,filename) in files - # Create temporal modules to isolate and protect test scopes - tmpdir = mktempdir(;cleanup=true) - tmpmod = split(filename,".")[1] - tmpfile = joinpath(tmpdir,tmpmod) - isfile(tmpfile) && error("File $tmpfile already exists!") - testpath = joinpath(@__DIR__,"../src", filename) - open(tmpfile,"w") do f - println(f, "# This file is automatically generated") - println(f, "# Do not edit") - println(f) - println(f, "module $tmpmod include(\"$testpath\") end") - end - # @time @testset "$title" begin include(tmpfile) end -end diff --git a/test/runtests.jl b/test/runtests.jl index 7456f93..1ad6452 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,7 @@ function get_title(files,filename) error("File $filename not found!") end -if (length(ARGS) != 0) +if !iszero(length(ARGS)) files = [get_title(Tutorials.files,filename)=>filename for filename in ARGS] else files = Tutorials.files From 7fdb8d954003ccfe20c0173360e7509a84770a8c Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 12:32:31 +1000 Subject: [PATCH 07/10] Minor --- .github/workflows/ci.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9f19435..8832b6d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,7 +42,6 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - run: | cd test/ - julia --project=.. maketests.jl julia --project=.. --color=yes --check-bounds=yes runtests.jl ${{ matrix.test_set.files }} tutorials_mpi: name: Tutorials MPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} From 56168255ab80e84501f982310bb60bfe4e6f2ca1 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 12:37:44 +1000 Subject: [PATCH 08/10] Minor --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index b51c713..e2b01a6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Santiago Badia ", "Francesc Verdugo Date: Mon, 16 Jun 2025 12:54:51 +1000 Subject: [PATCH 09/10] Fix paths in advection-diffusion --- src/advection_diffusion.jl | 21 +++++++++++++-------- src/stokes_blocks.jl | 6 +++--- 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/src/advection_diffusion.jl b/src/advection_diffusion.jl index af58d8e..6900ef3 100644 --- a/src/advection_diffusion.jl +++ b/src/advection_diffusion.jl @@ -82,20 +82,25 @@ # The meshed model is like following: -# ![Pentagon_2D_Mesh](Pentagon_2D_Mesh.png) +# ![Pentagon_2D_Mesh](../assets/advection_diffusion/Pentagon_2D_Mesh.png) # ## Numerical Scheme -# The weak form of a PDE is a reformulation that allows the problem to be solved in a broader function space. Instead of requiring the solution to satisfy the PDE at every point (as in the strong form), the weak form requires the solution to satisfy an integral equation. This makes it particularly suitable for numerical methods such as the Finite Element Method. +# The weak form of a PDE is a reformulation that allows the problem to be solved in a broader +# function space. Instead of requiring the solution to satisfy the PDE at every point (as in the strong form), +# the weak form requires the solution to satisfy an integral equation. This makes it particularly suitable for +# numerical methods such as the Finite Element Method. -# Since we already hcave the original PDE (strong form), we can multiply each side by a test function $v \in H^1(\Omega)$(functions with square-integrable first derivatives). $v$ satisfies the Dirichlet boundary condition of $0$. Make integral on both sides. +# Since we already hcave the original PDE (strong form), we can multiply each side by a test +# function $v \in H^1(\Omega)$(functions with square-integrable first derivatives). $v$ satisfies +# the Dirichlet boundary condition of $0$. Make integral on both sides. # The weak form associated with this formulation is, find $T \in H^1(\Omega)$, such that: -# $$ +# ```math # \int_\Omega v (\mathbf{u} \cdot \nabla T) \, \mathrm{d}\Omega # + \int_\Omega D \nabla T \cdot \nabla v \, \mathrm{d}\Omega # = 0 -# $$ +# ``` # ## FE spaces @@ -106,7 +111,7 @@ using GridapGmsh # Import the model from file using GmshDiscreteModel: -model_pentagon = GmshDiscreteModel("pentagon_mesh.msh") +model_pentagon = GmshDiscreteModel("../models/pentagon_mesh.msh") # Set up the test FE space $V_0$, which conforms the zero value boundary conditions. @@ -165,8 +170,8 @@ writevtk(Ω,"results_non_zero",cellfields=["uh_non_zero"=>uh_non_zero]) # We can use the ParaView to preview the results clearly. Here is the temperature distribution without any flow velocity: -# ![Result_zero](Result_zero.png) +# ![Result_zero](../assets/advection_diffusion/Result_zero.png) # Here is the temperature distribution with a flow velocity of 2 going up: -# ![Result_zero](Result_non_zero.png) +# ![Result_zero](../assets/advection_diffusion/Result_non_zero.png) diff --git a/src/stokes_blocks.jl b/src/stokes_blocks.jl index 8523e03..9327e57 100644 --- a/src/stokes_blocks.jl +++ b/src/stokes_blocks.jl @@ -41,9 +41,9 @@ # ``` # # where: -# - $$A$$: Vector Laplacian (velocity block) -# - $$B$$: Divergence operator -# - $$B^T$$: Gradient operator +# - $A$: Vector Laplacian (velocity block) +# - $B$: Divergence operator +# - $B^T$: Gradient operator # # ## Solution strategy # From 4bc721ec36f8866bb661aad2c8a8711f0c7aa47d Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 16 Jun 2025 13:03:19 +1000 Subject: [PATCH 10/10] Minor --- src/poisson_hho.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/poisson_hho.jl b/src/poisson_hho.jl index 43186b6..e1967e4 100644 --- a/src/poisson_hho.jl +++ b/src/poisson_hho.jl @@ -41,7 +41,8 @@ using Gridap.Arrays # ## Geometry # -# We generate a 2-dimensional simplicial mesh from a Cartesian grid: +# We generate a 2-dimensional simplicial mesh from a Cartesian grid, then taking its Voronoi +# dual to create a polytopal mesh. u(x) = sin(2*π*x[1])*sin(2*π*x[2]) f(x) = -Δ(u)(x)