From fb9c982fee0768463448955e00473205919a4a58 Mon Sep 17 00:00:00 2001
From: Fredrik Ekre <ekrefredrik@gmail.com>
Date: Thu, 29 Sep 2022 18:14:47 +0200
Subject: [PATCH] Add mean value constraint

---
 docs/src/literate/stokes-flow.jl | 28 +++++++++++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/docs/src/literate/stokes-flow.jl b/docs/src/literate/stokes-flow.jl
index 52bf910e78..69db24c068 100644
--- a/docs/src/literate/stokes-flow.jl
+++ b/docs/src/literate/stokes-flow.jl
@@ -45,7 +45,12 @@
 # \boldsymbol{u} = \boldsymbol{0} \quad \forall \boldsymbol{x}\ \in\
 # \Gamma_2 \cup \Gamma_4 := \{ \boldsymbol{x}:\ ||\boldsymbol{x}|| \in \{0.5, 1\}\}.
 # ```
-# No boundary conditions are needed for the pressure.
+# With this formulation and boundary conditions for ``\boldsymbol{u}`` the pressure will
+# only be determined up to a constant. We will add an additional constraint which fixes this
+# constant. In particular, we will enforce that the mean value of the pressure on the
+# boundary is 0, see [deal.ii
+# step-11](https://www.dealii.org/current/doxygen/deal.II/step_11.html) for some more
+# discussion around this.
 #
 # The corresponding weak form reads as follows: Find ``(\boldsymbol{u}, p) \in \mathbb{U}
 # \times \mathrm{L}_2`` s.t.
@@ -194,6 +199,15 @@ end
 # For the remaining part of the boundary we add a homogeneous Dirichlet boundary condition
 # on both components of the velocity field. This is done using the [`Dirichlet`](@ref)
 # constructor, which we have discussed in other tutorials.
+#
+# As discussed in the introduction, we will add a constraint to the mean value of the
+# pressure. We do this by constructing an [`AffineConstraint`](@ref) for the pressure dofs
+# on the boundary:
+# ```math
+# 0 = \int_{\partial\Omega} p\ \mathrm{d}\Gamma \approx \sum_{i=1}^{D} p_i \quad \Rightarrow
+# \quad p_1 = \sum_{i=2}^{D} -p_i,
+# ```
+# where ``\{p_i: i \in 1..D\}`` are the pressure dofs on the boundary.
 
 function setup_constraints(dh)
     ch = ConstraintHandler(dh)
@@ -206,6 +220,18 @@ function setup_constraints(dh)
     Γ24 = union(getfaceset(dh.grid, "Γ2"), getfaceset(dh.grid, "Γ4"))
     dbc = Dirichlet(:u, Γ24, (x, t) -> [0, 0], [1, 2])
     add!(ch, dbc)
+    ## Mean value constraint: Ensure that \sum
+    ch_mean = ConstraintHandler(dh)
+    ∂Ω = union(getfaceset(dh.grid, "Γ1"), getfaceset(dh.grid, "Γ2"),
+               getfaceset(dh.grid, "Γ3"), getfaceset(dh.grid, "Γ4"))
+    add!(ch_mean, Dirichlet(:p, ∂Ω, (x, t) -> 0))
+    close!(ch_mean)
+    mean_value_constraint = AffineConstraint(
+        pop!(ch_mean.prescribed_dofs),
+        [d => -1. for d in ch_mean.prescribed_dofs],
+        0.0
+    )
+    add!(ch, mean_value_constraint)
     ## Finalize
     close!(ch)
     update!(ch, 0)