Lagrange multipliers tutorial

.github/workflows/ci.yml

@@ -25,7 +25,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
+      - 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 }}
@@ -46,7 +46,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 p_laplacian.jl hyperelasticity.jl dg_discretization.jl fsi_tutorial.jl
+      - 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 }}

Project.toml

@@ -7,6 +7,7 @@ version = "0.17.0"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
@@ -33,6 +34,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+DataStructures = "0.18.22"
 Gridap = "0.18"
 GridapDistributed = "0.4"
 GridapGmsh = "0.7"

deps/build.jl

@@ -26,6 +26,7 @@ files = [
   "Topology optimization"=>"TopOptEMFocus.jl",
   "Poisson on unfitted meshes"=>"poisson_unfitted.jl",
   "Poisson with AMR"=>"poisson_amr.jl",
+  "Lagrange multipliers" => "lagrange_multipliers.jl",
   "Low-level API - Poisson equation"=>"poisson_dev_fe.jl",
   "Low-level API - Geometry" => "geometry_dev.jl",
 ]

src/lagrange_multipliers.jl

@@ -0,0 +1,127 @@
+# In this tutorial, we will learn
+#
+#    - How to enforce constraints using Lagrange multipliers
+#    - How to work with `ConstantFESpace`
+#
+# ## Problem statement
+#
+# In this tutorial, we solve the Poisson equation with pure Neumann boundary conditions. 
+# This problem is well-known to be singular since the solution is defined up to a constant, 
+# which we have to fix to obtain a unique solution.
+# Here, we will use a Lagrange multiplier to enforce that the mean value of the solution 
+# equals a given constant.
+#
+# The problem reads: find $u$ and $λ$ such that
+#
+# ```math
+# \left\lbrace
+# \begin{aligned}
+# -\Delta u = f  \ &\text{in} \ \Omega,\\
+# \nabla u\cdot n = g \ &\text{on}\ \Gamma,\\
+# \int_{\Omega} u \ {\rm d}\Omega = \bar{u},\\
+# \end{aligned}
+# \right.
+# ```
+#
+# where $\Omega$ is our domain, $\Gamma$ is its boundary, $n$ is the outward unit normal vector,
+# and $\bar{u}$ is a given constant that fixes the mean value of the solution.
+#
+# ## Numerical scheme
+#
+# The weak form of this problem using Lagrange multipliers reads:
+# find $(u,λ) \in V \times \Lambda$ such that
+#
+# ```math
+# \begin{aligned}
+# \int_{\Omega} \nabla u \cdot \nabla v \ {\rm d}\Omega + 
+# \int_{\Omega} λv \ {\rm d}\Omega + 
+# \int_{\Omega} uμ \ {\rm d}\Omega = 
+# \int_{\Omega} fv \ {\rm d}\Omega + 
+# \int_{\Gamma} v(g\cdot n) \ {\rm d}\Gamma + 
+# \int_{\Omega} μ\bar{u} \ {\rm d}\Omega
+# \end{aligned}
+# ```
+#
+# for all $(v,μ) \in V \times \Lambda$, where $V = H^1(\Omega)$ and $\Lambda = \mathbb{R}$.
+#
+# ## Implementation
+#
+# First, we load the Gridap package and define the exact solution that we will use to 
+# manufacture the source term and boundary condition:
+
+using Gridap
+
+u_exact(x) = sin(x[1]) * cos(x[2])
+
+# Now we can create a simple Cartesian mesh of the unit square:
+
+model = CartesianDiscreteModel((0,1,0,1),(8,8))
+
+# We will use first order Lagrangian finite elements for the primal variable u.
+
+order = 1
+reffe = ReferenceFE(lagrangian, Float64, order)
+V = FESpace(model, reffe)
+
+# For the Lagrange multiplier λ, we need a space of constant functions, since λ ∈ ℝ.
+# In Gridap, we can create such a space using `ConstantFESpace`:
+
+Λ = ConstantFESpace(model)
+
+# Conceptually, a `ConstantFESpace` is a space defined on the whole domain with a 
+# single degree of freedom, which is what we need for the Lagrange multiplier λ.
+# We finally bundle both spaces into a multi-field space:
+
+X = MultiFieldFESpace([V, Λ])
+
+# ## Integration
+#
+# We need to create the triangulation and measures for both domain and boundary 
+# integration:
+
+Ω = Triangulation(model)
+Γ = BoundaryTriangulation(model)
+dΩ = Measure(Ω, 2*order)
+dΓ = Measure(Γ, 2*order)
+
+# Next, we manufacture the source term f and Neumann boundary condition g 
+# from the exact solution. We also compute the mean value ū that we want 
+# to enforce:
+
+f(x) = -Δ(u_exact)(x)
+g(x) = ∇(u_exact)(x)
+ū = sum(∫(u_exact)dΩ)
+nΓ = get_normal_vector(Γ)
+
+# ## Weak Form
+#
+# We can now define the bilinear and linear forms of our problem.
+# Note how the forms take tuples as arguments, representing the 
+# multi-field nature of our solution:
+
+a((u,λ),(v,μ)) = ∫(∇(u)⋅∇(v) + λ*v + u*μ)dΩ
+l((v,μ)) = ∫(f*v + μ*ū)dΩ + ∫(v*(g⋅nΓ))*dΓ
+
+# ## Solution
+#
+# We can now create the FE operator and solve the system:
+
+op = AffineFEOperator(a, l, X, X)
+uh, λh = solve(op)
+
+# Note how we get two values from solve: the primal solution uh and 
+# the Lagrange multiplier λh. Finally, we compute the L2 error and 
+# verify that the mean value constraint is satisfied:
+
+eh = uh - u_exact
+l2_error = sqrt(sum(∫(eh⋅eh)*dΩ))
+ūh = sum(∫(uh)*dΩ)
+
+# The L2 error should be small (of order h²) and ūh should be very close to ū,
+# showing that both the equation and the constraint are well satisfied.
+
+# ## Visualization
+#
+# We can visualize the solution and error by writing them to a VTK file:
+
+writevtk(Ω, "results", cellfields=["uh"=>uh, "error"=>eh])

src/poisson_amr.jl

@@ -79,7 +79,11 @@ function LShapedModel(n)
   cell_coords = map(mean,get_cell_coordinates(model))
   l_shape_filter(x) = (x[1] < 0.5) || (x[2] < 0.5)
   mask = map(l_shape_filter,cell_coords)
-  return simplexify(DiscreteModelPortion(model,mask))
+  model = simplexify(DiscreteModelPortion(model,mask))
+
+  grid = get_grid(model)
+  topo = get_grid_topology(model)
+  return UnstructuredDiscreteModel(grid, topo, FaceLabeling(topo))
 end
 
 # Define the L2 norm for error estimation.
@@ -179,8 +183,8 @@ for i in 1:nsteps
   )
 
   println("Error: $error, Error η: $(sum(η))")
-  last_error = error
-  model = fmodel
+  global last_error = error
+  global model = fmodel
 end
 
 # The final mesh gives the following result:

test/runtests.jl

@@ -17,23 +17,17 @@ else
 end
 
 for (title,filename) in files
-    # Create temporal modules to isolate and protect test scopes
-    tmpdir=mktempdir(;cleanup=true)
-    filename_wo_extension=split(filename,".")[1]
-    tmpmod = filename_wo_extension
-    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
+  # 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
-
-# module fsi_tutorial
-# using Test
-# @time @testset "fsi_tutorial" begin include("../src/fsi_tutorial.jl") end
-# end # module