Gridap 0.19

.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,21 @@ 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: 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,72 +40,11 @@ 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
+      - run: | 
+          cd test/
+          julia --project=.. --color=yes --check-bounds=yes runtests.jl ${{ matrix.test_set.files }}
   tutorials_mpi:
-    name: TutorialsMPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Tutorials MPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -109,7 +63,9 @@ 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_mpi.jl
+      - run: | 
+          cd test/
+          julia --project=.. --color=yes --check-bounds=yes runtests_mpi.jl
   docs:
     name: Documentation
     runs-on: ubuntu-latest

Project.toml

@@ -4,6 +4,7 @@ authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <fver
 version = "0.17.0"
 
 [deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -17,6 +18,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"
@@ -34,12 +36,14 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+BlockArrays = "1.6.3"
 DataStructures = "0.18.22"
-Gridap = "0.18"
+Gridap = "0.19"
 GridapDistributed = "0.4"
 GridapGmsh = "0.7"
 GridapP4est = "0.3"
 GridapPETSc = "0.5"
+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/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)

src/poisson_hdg.jl

@@ -0,0 +1,164 @@
+# 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 m on mesh interfaces to reduce the global system size.
+#
+# ## HDG Discretization
+#
+# 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}
+# ```
+#
+# 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``
+# - ``m_h``: the approximation to the trace of ``u`` on element faces
+#
+# Numerical fluxes are defindes as 
+#
+# ```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:
+
+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)$
+
+# ## Geometry
+#
+# 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))
+
+# 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
+
+Ω = 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
+
+# ## FESpaces
+#
+# HDG uses three different finite element spaces:
+# 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)
+
+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
+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 = BlockMultiFieldStyle(2,(2,1))  # Special blocking for efficient static condensation
+X = MultiFieldFESpace([V, Q, N]; style=mfs)
+
+# ## Weak Form and integration
+#
+
+degree = 2*(order+1)  # Integration degree
+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
+
+n = get_normal_vector(Γp)  # Face normal vector
+Πn(u) = u⋅n  # Normal component
+Π(u) = change_domain(u,Γp,DomainStyle(u))  # Project to skeleton
+
+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
+#
+# 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)
+
+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. 
+#