Added Lagrange multiplier tutorial

Author: JordiManyer

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])