From 385ab165045b1f40182d0ba3d9f2751992698326 Mon Sep 17 00:00:00 2001
From: "G. D. McBain" <gdmcbain@protonmail.com>
Date: Wed, 26 May 2021 16:35:45 +1000
Subject: [PATCH] constrained formulation of the Stokes problem

---
 docs/examples/ex30_minimize.py | 83 ++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)
 create mode 100644 docs/examples/ex30_minimize.py

diff --git a/docs/examples/ex30_minimize.py b/docs/examples/ex30_minimize.py
new file mode 100644
index 000000000..4af06aa74
--- /dev/null
+++ b/docs/examples/ex30_minimize.py
@@ -0,0 +1,83 @@
+r"""Stokes problem by minimization
+
+This redoes Example 30 using scipy.optimize.minimize, as in ./ex01_minimize.py.
+
+Based on the 'constraint formulation of the Stokes problem' (Ern & Guermond 2004, 
+&sect;4.1.3), it suffices to minimize the viscous dissipation subject to the
+constraint that the velocity field have no divergence.
+
+.. [ERN-GUERMOND] 
+
+"""
+from skfem import *
+from skfem.models.poisson import vector_laplace, laplace, mass
+from skfem.models.general import divergence, rot
+
+import numpy as np
+from scipy.optimize import LinearConstraint, minimize
+from scipy.sparse import csr_matrix
+from scipy.sparse.linalg import LinearOperator, minres
+
+mesh = MeshQuad.init_tensor(*(np.linspace(-.5, .5, 2**6),)*2)
+element = {'u': ElementVectorH1(ElementQuad2()),
+           'p': ElementQuad1()}
+basis = {variable: InteriorBasis(mesh, e, intorder=3)
+         for variable, e in element.items()}
+
+
+@LinearForm
+def body_force(v, w):
+    return w.x[0] * v.value[1]
+
+
+A = asm(vector_laplace, basis['u'])
+B = -asm(divergence, basis['u'], basis['p'])
+f = asm(body_force, basis['u'])
+D = basis['u'].find_dofs()['all'].all()
+
+
+def functional(u: np.ndarray) -> float:
+    return u @ A @ u / 2 - f @ u
+
+
+def jacobian(u: np.ndarray) -> np.ndarray:
+    return u @ A - f
+
+
+def hessian(u: np.ndarray) -> np.ndarray:
+    return A
+
+
+incompressibility = LinearConstraint(B, 0, 0)
+slip = csr_matrix((np.ones(D.size), (np.arange(D.size), D)), (D.size, basis["u"].N))
+noslip = LinearConstraint(slip, 0, 0)
+
+velocity = minimize(
+    functional,
+    basis["u"].zeros(),
+    jac=jacobian,
+    hess=hessian,
+    method="trust-constr",
+    constraints=[incompressibility, noslip],
+    options={"verbose": 3},
+).x
+
+basis['psi'] = basis['u'].with_element(ElementQuad2())
+vorticity = asm(rot, basis['psi'], w=basis['u'].interpolate(velocity))
+psi = solve(*condense(asm(laplace, basis['psi']), vorticity, D=basis['psi'].find_dofs()))
+psi0 = (basis['psi'].probes(np.zeros((2, 1))) @ psi)[0]
+
+
+if __name__ == '__main__':
+
+    from pathlib import Path
+
+    from matplotlib.tri import Triangulation
+    from skfem.visuals.matplotlib import draw
+
+    print(psi0)
+
+    ax = draw(mesh)
+    ax.tricontour(Triangulation(*mesh.p, mesh.to_meshtri().t.T),
+                  psi[basis['psi'].nodal_dofs.flatten()])
+    ax.get_figure().savefig(Path(__file__).with_suffix(".png"))