|
| 1 | +# In this tutorial, we will learn |
| 2 | +# |
| 3 | +# - How to enforce constraints using Lagrange multipliers |
| 4 | +# - How to work with `ConstantFESpace` |
| 5 | +# |
| 6 | +# ## Problem statement |
| 7 | +# |
| 8 | +# In this tutorial, we solve the Poisson equation with pure Neumann boundary conditions. |
| 9 | +# This problem is well-known to be singular since the solution is defined up to a constant, |
| 10 | +# which we have to fix to obtain a unique solution. |
| 11 | +# Here, we will use a Lagrange multiplier to enforce that the mean value of the solution |
| 12 | +# equals a given constant. |
| 13 | +# |
| 14 | +# The problem reads: find $u$ and $λ$ such that |
| 15 | +# |
| 16 | +# ```math |
| 17 | +# \left\lbrace |
| 18 | +# \begin{aligned} |
| 19 | +# -\Delta u = f \ &\text{in} \ \Omega,\\ |
| 20 | +# \nabla u\cdot n = g \ &\text{on}\ \Gamma,\\ |
| 21 | +# \int_{\Omega} u \ {\rm d}\Omega = \bar{u},\\ |
| 22 | +# \end{aligned} |
| 23 | +# \right. |
| 24 | +# ``` |
| 25 | +# |
| 26 | +# where $\Omega$ is our domain, $\Gamma$ is its boundary, $n$ is the outward unit normal vector, |
| 27 | +# and $\bar{u}$ is a given constant that fixes the mean value of the solution. |
| 28 | +# |
| 29 | +# ## Numerical scheme |
| 30 | +# |
| 31 | +# The weak form of this problem using Lagrange multipliers reads: |
| 32 | +# find $(u,λ) \in V \times \Lambda$ such that |
| 33 | +# |
| 34 | +# ```math |
| 35 | +# \begin{aligned} |
| 36 | +# \int_{\Omega} \nabla u \cdot \nabla v \ {\rm d}\Omega + |
| 37 | +# \int_{\Omega} λv \ {\rm d}\Omega + |
| 38 | +# \int_{\Omega} uμ \ {\rm d}\Omega = |
| 39 | +# \int_{\Omega} fv \ {\rm d}\Omega + |
| 40 | +# \int_{\Gamma} v(g\cdot n) \ {\rm d}\Gamma + |
| 41 | +# \int_{\Omega} μ\bar{u} \ {\rm d}\Omega |
| 42 | +# \end{aligned} |
| 43 | +# ``` |
| 44 | +# |
| 45 | +# for all $(v,μ) \in V \times \Lambda$, where $V = H^1(\Omega)$ and $\Lambda = \mathbb{R}$. |
| 46 | +# |
| 47 | +# ## Implementation |
| 48 | +# |
| 49 | +# First, we load the Gridap package and define the exact solution that we will use to |
| 50 | +# manufacture the source term and boundary condition: |
| 51 | + |
| 52 | +using Gridap |
| 53 | + |
| 54 | +u_exact(x) = sin(x[1]) * cos(x[2]) |
| 55 | + |
| 56 | +# Now we can create a simple Cartesian mesh of the unit square: |
| 57 | + |
| 58 | +model = CartesianDiscreteModel((0,1,0,1),(8,8)) |
| 59 | + |
| 60 | +# We will use first order Lagrangian finite elements for the primal variable u. |
| 61 | + |
| 62 | +order = 1 |
| 63 | +reffe = ReferenceFE(lagrangian, Float64, order) |
| 64 | +V = FESpace(model, reffe) |
| 65 | + |
| 66 | +# For the Lagrange multiplier λ, we need a space of constant functions, since λ ∈ ℝ. |
| 67 | +# In Gridap, we can create such a space using `ConstantFESpace`: |
| 68 | + |
| 69 | +Λ = ConstantFESpace(model) |
| 70 | + |
| 71 | +# Conceptually, a `ConstantFESpace` is a space defined on the whole domain with a |
| 72 | +# single degree of freedom, which is what we need for the Lagrange multiplier λ. |
| 73 | +# We finally bundle both spaces into a multi-field space: |
| 74 | + |
| 75 | +X = MultiFieldFESpace([V, Λ]) |
| 76 | + |
| 77 | +# ## Integration |
| 78 | +# |
| 79 | +# We need to create the triangulation and measures for both domain and boundary |
| 80 | +# integration: |
| 81 | + |
| 82 | +Ω = Triangulation(model) |
| 83 | +Γ = BoundaryTriangulation(model) |
| 84 | +dΩ = Measure(Ω, 2*order) |
| 85 | +dΓ = Measure(Γ, 2*order) |
| 86 | + |
| 87 | +# Next, we manufacture the source term f and Neumann boundary condition g |
| 88 | +# from the exact solution. We also compute the mean value ū that we want |
| 89 | +# to enforce: |
| 90 | + |
| 91 | +f(x) = -Δ(u_exact)(x) |
| 92 | +g(x) = ∇(u_exact)(x) |
| 93 | +ū = sum(∫(u_exact)dΩ) |
| 94 | +nΓ = get_normal_vector(Γ) |
| 95 | + |
| 96 | +# ## Weak Form |
| 97 | +# |
| 98 | +# We can now define the bilinear and linear forms of our problem. |
| 99 | +# Note how the forms take tuples as arguments, representing the |
| 100 | +# multi-field nature of our solution: |
| 101 | + |
| 102 | +a((u,λ),(v,μ)) = ∫(∇(u)⋅∇(v) + λ*v + u*μ)dΩ |
| 103 | +l((v,μ)) = ∫(f*v + μ*ū)dΩ + ∫(v*(g⋅nΓ))*dΓ |
| 104 | + |
| 105 | +# ## Solution |
| 106 | +# |
| 107 | +# We can now create the FE operator and solve the system: |
| 108 | + |
| 109 | +op = AffineFEOperator(a, l, X, X) |
| 110 | +uh, λh = solve(op) |
| 111 | + |
| 112 | +# Note how we get two values from solve: the primal solution uh and |
| 113 | +# the Lagrange multiplier λh. Finally, we compute the L2 error and |
| 114 | +# verify that the mean value constraint is satisfied: |
| 115 | + |
| 116 | +eh = uh - u_exact |
| 117 | +l2_error = sqrt(sum(∫(eh⋅eh)*dΩ)) |
| 118 | +ūh = sum(∫(uh)*dΩ) |
| 119 | + |
| 120 | +# The L2 error should be small (of order h²) and ūh should be very close to ū, |
| 121 | +# showing that both the equation and the constraint are well satisfied. |
| 122 | + |
| 123 | +# ## Visualization |
| 124 | +# |
| 125 | +# We can visualize the solution and error by writing them to a VTK file: |
| 126 | + |
| 127 | +writevtk(Ω, "results", cellfields=["uh"=>uh, "error"=>eh]) |
0 commit comments