Residual-based a posteriori control of L2 error

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm56

open PDETextbook
open Classical

abbrev E2 := Euclidean 2

/-- Closed triangle K: convex hull of the vertices of a nondegenerate 2-simplex. -/
def tri (K : Affine.Simplex ℝ E2 2) : Set E2 := convexHull ℝ (Set.range K.points)

/-- Edge of K opposite vertex i. -/
def edge (K : Affine.Simplex ℝ E2 2) (i : Fin 3) : Set E2 :=
  segment ℝ (K.points (i + 1)) (K.points (i + 2))

/-- h_K = diam(K). -/
def hK (K : Affine.Simplex ℝ E2 2) : ℝ := Metric.diam (tri K)

/-- Convex polygonal domain in the plane: open, nonempty, convex, with polygonal boundary
(closure is the convex hull of finitely many points). -/
def ConvexPolygonalDomain (Ω : Set E2) : Prop :=
  IsOpen Ω ∧ Ω.Nonempty ∧ Convex ℝ Ω ∧
    ∃ P : Finset E2, closure Ω = convexHull ℝ (P : Set E2)

/-- Triangulation of Ω (Sect. 5.2): closed triangles covering closure Ω; two distinct triangles
meet in the empty set, a common node, or a common edge; no node lies in the interior of an edge. -/
def IsTriangulation (Ω : Set E2) (T : Finset (Affine.Simplex ℝ E2 2)) : Prop :=
  T.Nonempty ∧ closure Ω = (⋃ K ∈ T, tri K) ∧
  (∀ K ∈ T, ∀ K' ∈ T, K ≠ K' →
    (tri K ∩ tri K' = ∅ ∨
     (∃ P, P ∈ Set.range K.points ∧ P ∈ Set.range K'.points ∧ tri K ∩ tri K' = {P}) ∨
     (∃ P Q, P ≠ Q ∧ P ∈ Set.range K.points ∧ P ∈ Set.range K'.points ∧
        Q ∈ Set.range K.points ∧ Q ∈ Set.range K'.points ∧
        tri K ∩ tri K' = segment ℝ P Q))) ∧
  (∀ K ∈ T, ∀ K' ∈ T, ∀ i j : Fin 3,
    K'.points j ∉ openSegment ℝ (K.points (i + 1)) (K.points (i + 2)))

/-- h = max_K h_K. -/
def IsMeshsize (T : Finset (Affine.Simplex ℝ E2 2)) (h : ℝ) : Prop :=
  (∀ K ∈ T, hK K ≤ h) ∧ ∃ K ∈ T, hK K = h

/-- All angles of K are at least θ0. -/
def MinAngleGE (K : Affine.Simplex ℝ E2 2) (θ0 : ℝ) : Prop :=
  ∀ i j k : Fin 3, i ≠ j → i ≠ k → j ≠ k →
    θ0 ≤ EuclideanGeometry.angle (K.points j) (K.points i) (K.points k)

/-- S_h: continuous on closure Ω, linear (affine) on each K, vanishing on Γ. -/
def InSh (Ω : Set E2) (T : Finset (Affine.Simplex ℝ E2 2)) (v : E2 → ℝ) : Prop :=
  ContinuousOn v (closure Ω) ∧
  (∀ K ∈ T, ∃ L : E2 →ᵃ[ℝ] ℝ, Set.EqOn v (fun x => L x) (tri K)) ∧
  ∀ x ∈ frontier Ω, v x = 0

/-- a(v,w) = ∫_Ω a ∇v·∇w dx, with weak gradients G, H. -/
def aForm (Ω : Set E2) (a : E2 → ℝ) (G H : Fin 2 → E2 → ℝ) : ℝ :=
  ∫ x in Ω, a x * ∑ i, G i x * H i x

/-- u solves (5.24): u ∈ H₀¹ and a(u,v) = (f,v) for all v ∈ H₀¹. -/
def IsWeakSolution (Ω : Set E2) (a f u : E2 → ℝ) : Prop :=
  InH10Closure Ω u ∧ ∃ G, HasH1Gradient Ω u G ∧
    ∀ (v : E2 → ℝ) (Gv : Fin 2 → E2 → ℝ), InH10Closure Ω v → HasH1Gradient Ω v Gv →
      aForm Ω a G Gv = ∫ x in Ω, f x * v x

/-- u_h solves (5.26): u_h ∈ S_h and a(u_h,χ) = (f,χ) for all χ ∈ S_h. -/
def IsFESolution (Ω : Set E2) (T : Finset (Affine.Simplex ℝ E2 2)) (a f uh : E2 → ℝ) : Prop :=
  InSh Ω T uh ∧ ∃ G, HasH1Gradient Ω uh G ∧
    ∀ (χ : E2 → ℝ) (Gχ : Fin 2 → E2 → ℝ), InSh Ω T χ → HasH1Gradient Ω χ Gχ →
      aForm Ω a G Gχ = ∫ x in Ω, f x * χ x

/-- 𝒜 u_h = -∇·(a∇u_h) on K, computed classically from the affine piece L of u_h on K. -/
def elemOp (a : E2 → ℝ) (L : E2 →ᵃ[ℝ] ℝ) (x : E2) : ℝ :=
  -∑ i : Fin 2, partialDeriv (fun y => a y * partialDeriv (fun z => L z) i y) i x

/-- Outward unit normal of K on the edge opposite vertex i. -/
def edgeNormal (K : Affine.Simplex ℝ E2 2) (i : Fin 3) : E2 :=
  let p := K.points (i + 1)
  let d := K.points (i + 2) - p
  let w := p - K.points i
  let m := w - (inner ℝ w d / inner ℝ d d) • d
  ‖m‖⁻¹ • m

/-- Jump [n·∇u_h] across the edge opposite vertex i of K (neighbour K' sharing that edge). -/
def normalJump (T : Finset (Affine.Simplex ℝ E2 2))
    (L : Affine.Simplex ℝ E2 2 → (E2 →ᵃ[ℝ] ℝ)) (K : Affine.Simplex ℝ E2 2) (i : Fin 3) : ℝ :=
  ∑ K' ∈ T, if K' ≠ K ∧ edge K i ⊆ tri K' then
    (L K).linear (edgeNormal K i) - (L K').linear (edgeNormal K i) else 0

/-- ‖a[n·∇u_h]‖_{∂K∖Γ}, L² w.r.t. arc length (1-dim Hausdorff measure). -/
def jumpNorm (Ω : Set E2) (T : Finset (Affine.Simplex ℝ E2 2)) (a : E2 → ℝ)
    (L : Affine.Simplex ℝ E2 2 → (E2 →ᵃ[ℝ] ℝ)) (K : Affine.Simplex ℝ E2 2) : ℝ :=
  Real.sqrt (∑ i : Fin 3,
    ((eLpNorm (fun x => a x * normalJump T L K i) 2
      ((Measure.hausdorffMeasure (1 : ℝ)).restrict (edge K i \ frontier Ω))).toReal) ^ 2)

/-- R_K = h_K² ‖𝒜u_h − f‖_K + h_K^{3/2} ‖a[n·∇u_h]‖_{∂K∖Γ}. -/
def RK (Ω : Set E2) (T : Finset (Affine.Simplex ℝ E2 2)) (a f : E2 → ℝ)
    (L : Affine.Simplex ℝ E2 2 → (E2 →ᵃ[ℝ] ℝ)) (K : Affine.Simplex ℝ E2 2) : ℝ :=
  hK K ^ 2 * (eLpNorm (fun x => elemOp a (L K) x - f x) 2 (volume.restrict (tri K))).toReal
  + hK K ^ (3 / 2 : ℝ) * jumpNorm Ω T a L K

/-- Theorem 5.6 (a posteriori estimate).

Standing assumptions of Sect. 5.2: a is smooth on closure Ω with a ≥ a0 > 0, and f ∈ L₂.

The family convention of Sect. 5.3 says the triangulations {T_h}_{0<h<1} have angles bounded
below independently of h. It is realised by fixing a lower angle bound θ0 > 0. The constant C
is then uniform over every conforming triangulation T of Ω with meshsize h = max_K h_K ∈ (0,1)
and all angles ≥ θ0, and over all f.

This quantification is equivalent to requiring a uniform C for each family with an arbitrary
index set: the class of all such T is itself a family. It does not require a triangulation of
meshsize exactly h for every real h, so no size restriction on Ω is imposed.

L K is the linear piece of u_h on K. -/
theorem thm_5_6 (Ω : Set E2) (hΩ : ConvexPolygonalDomain Ω)
    (a : E2 → ℝ) (ha : ContDiffOn ℝ (⊤ : ℕ∞) a (closure Ω))
    (ha0 : ∃ a0 : ℝ, 0 < a0 ∧ ∀ x ∈ closure Ω, a0 ≤ a x) :
    ∀ θ0 : ℝ, 0 < θ0 → ∃ C : ℝ,
      ∀ T : Finset (Affine.Simplex ℝ E2 2), IsTriangulation Ω T →
      (∃ h : ℝ, IsMeshsize T h ∧ 0 < h ∧ h < 1) →
      (∀ K ∈ T, MinAngleGE K θ0) →
      ∀ f : E2 → ℝ, InL2 Ω f →
      ∀ u uh : E2 → ℝ, IsWeakSolution Ω a f u → IsFESolution Ω T a f uh →
      ∀ L : Affine.Simplex ℝ E2 2 → (E2 →ᵃ[ℝ] ℝ),
        (∀ K ∈ T, Set.EqOn uh (fun x => L K x) (tri K)) →
        (eLpNorm (fun x => uh x - u x) 2 (domainMeasure Ω)).toReal ≤
          C * Real.sqrt (∑ K ∈ T, (RK Ω T a f L K) ^ 2) := by
  sorry

end PDETextbook.LarssonThomee.Thm56