Streamline-diffusion error for stationary transport

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm13_6

abbrev E2 := Euclidean 2

/-- Bounded convex polygon Ω = {x | ⟪ν j, x⟫ < c j ∀ j}, with distinct unit outward normals ν j. -/
structure ConvexPolygon where
  m : ℕ
  ν : Fin m → E2
  c : Fin m → ℝ
  unit : ∀ j, ‖ν j‖ = 1
  inj : Function.Injective ν

def ConvexPolygon.dom (P : ConvexPolygon) : Set E2 := {x | ∀ j, inner ℝ (P.ν j) x < P.c j}

/-- Relative interior of the edge with outward normal ν j: only constraint j is active. -/
def ConvexPolygon.openFacet (P : ConvexPolygon) (j : Fin P.m) : Set E2 :=
  {x | inner ℝ (P.ν j) x = P.c j ∧ ∀ k, k ≠ j → inner ℝ (P.ν k) x < P.c k}

/-- Inflow boundary Γ₋ = {a·n < 0}, taken with its closure (a union of closed edges). -/
def inflowBdry (P : ConvexPolygon) (a : E2) : Set E2 :=
  closure (⋃ (j : Fin P.m) (_ : inner ℝ a (P.ν j) < 0), P.openFacet j)

/-- |v|²_{Γ₊} = ∫_{Γ₊} |a·n| v² ds, ds = 1-dim Hausdorff (arc length) measure. -/
def outflowNormSq (P : ConvexPolygon) (a : E2) (v : E2 → ℝ) : ℝ :=
  ∑ j, max 0 (inner ℝ a (P.ν j)) *
    ∫ x in P.openFacet j, v x ^ 2 ∂(Measure.hausdorffMeasure (1 : ℝ) : Measure E2)

def tri (t : Fin 3 → E2) : Set E2 := convexHull ℝ (Set.range t)

/-- Triangulation of Ω (Sect. 5.2): closed triangles covering Ω̄, pairwise meeting in
∅, a common node, or a common edge. -/
structure Triangulation (Ω : Set E2) where
  K : Finset (Fin 3 → E2)
  nondeg : ∀ t ∈ K, AffineIndependent ℝ t
  cover : closure Ω = ⋃ t ∈ K, tri t
  conform : ∀ t ∈ K, ∀ s ∈ K, t ≠ s →
    tri t ∩ tri s = ∅ ∨
    (∃ p ∈ Set.range t ∩ Set.range s, tri t ∩ tri s = {p}) ∨
    (∃ p ∈ Set.range t ∩ Set.range s, ∃ q ∈ Set.range t ∩ Set.range s,
      p ≠ q ∧ tri t ∩ tri s = segment ℝ p q)

variable {Ω : Set E2}

def nodes (T : Triangulation Ω) : Set E2 := ⋃ t ∈ T.K, Set.range t

/-- h = max_K diam K. -/
def meshSize (T : Triangulation Ω) : ℝ := ⨆ t : T.K, Metric.diam (tri t.1)

/-- All angles of all triangles are ≥ θ₀. -/
def AnglesBoundedBelow (T : Triangulation Ω) (θ₀ : ℝ) : Prop :=
  ∀ t ∈ T.K, ∀ i j k : Fin 3, i ≠ j → j ≠ k → i ≠ k →
    θ₀ ≤ EuclideanGeometry.angle (t j) (t i) (t k)

/-- The inflow boundary is exactly a union of triangle edges. -/
def MatchedInflow (T : Triangulation Ω) (Γm : Set E2) : Prop :=
  ∀ x ∈ Γm, ∃ t ∈ T.K, ∃ i j : Fin 3, i ≠ j ∧
    x ∈ segment ℝ (t i) (t j) ∧ segment ℝ (t i) (t j) ⊆ Γm

/-- S_h: continuous piecewise linear functions on T (no boundary conditions). -/
def InSh (T : Triangulation Ω) (v : E2 → ℝ) : Prop :=
  ContinuousOn v (closure Ω) ∧ ∀ t ∈ T.K, ∃ φ : E2 →ᵃ[ℝ] ℝ, Set.EqOn v φ (tri t)

/-- w = I_h g: the element of S_h with nodal values of g. -/
def IsInterp (T : Triangulation Ω) (g w : E2 → ℝ) : Prop :=
  InSh T w ∧ ∀ p ∈ nodes T, w p = g p

def nrm (Ω : Set E2) (v : E2 → ℝ) : ℝ := (eLpNorm v 2 (domainMeasure Ω)).toReal

/-- ‖u‖₂² = Σ_{|α|≤2} ‖D^α u‖², with G = ∇u and H i j = ∂_j ∂_i u (multi-indices: i ≤ j). -/
def H2NormSq (Ω : Set E2) (u : E2 → ℝ) (G : Fin 2 → E2 → ℝ) (H : Fin 2 → Fin 2 → E2 → ℝ) : ℝ :=
  nrm Ω u ^ 2 + ∑ i, nrm Ω (G i) ^ 2 + ∑ i, ∑ j, if i ≤ j then nrm Ω (H i j) ^ 2 else 0

/-- Theorem 13.6 (streamline diffusion error estimate). For a bounded convex polygon Ω,
constant a and a₀ > 0, a family of triangulations with angles bounded below and inflow
boundary matched to edges, there is C (independent of h, u, f, g) such that for the solution
u ∈ H² of (13.12) and u_h of (13.20), e = u_h − u satisfies
|e|_{Γ₊} + a₀^{1/2}‖e‖ + h^{1/2}‖a·∇e‖ ≤ C h^{3/2} ‖u‖₂. -/
theorem thm_13_6 (P : ConvexPolygon) (hbdd : Bornology.IsBounded P.dom)
    (hne : P.dom.Nonempty) (a : E2) (a0 : ℝ) (ha0 : 0 < a0)
    {ι : Type*} (T : ι → Triangulation P.dom) (θ₀ : ℝ) (hθ : 0 < θ₀)
    (hang : ∀ i, AnglesBoundedBelow (T i) θ₀)
    (hmesh : ∀ i, 0 < meshSize (T i) ∧ meshSize (T i) < 1)
    (hmatch : ∀ i, MatchedInflow (T i) (inflowBdry P a)) :
    ∃ C : ℝ, ∀ (i : ι) (f g u uh : E2 → ℝ) (G Gh : Fin 2 → E2 → ℝ)
      (H : Fin 2 → Fin 2 → E2 → ℝ),
      -- u ∈ H²(Ω) (continuous representative on Ω̄) solves (13.12)
      InL2 P.dom f →
      HasH1Gradient P.dom u G → (∀ k, HasH1Gradient P.dom (G k) (H k)) →
      ContinuousOn u (closure P.dom) →
      (∀ᵐ x ∂(domainMeasure P.dom), (∑ k, a k * G k x) + a0 * u x = f x) →
      (∀ x ∈ inflowBdry P a, u x = g x) →
      -- u_h ∈ S_h solves (13.20) with h = meshSize
      InSh (T i) uh → HasH1Gradient P.dom uh Gh →
      (∀ χ : E2 → ℝ, InSh (T i) χ → (∀ x ∈ inflowBdry P a, χ x = 0) →
        ∀ Gχ : Fin 2 → E2 → ℝ, HasH1Gradient P.dom χ Gχ →
          (∫ x in P.dom, ((∑ k, a k * Gh k x) + a0 * uh x) *
              (χ x + meshSize (T i) * ∑ k, a k * Gχ k x)) =
          ∫ x in P.dom, f x * (χ x + meshSize (T i) * ∑ k, a k * Gχ k x)) →
      (∃ Ig, IsInterp (T i) g Ig ∧ ∀ x ∈ inflowBdry P a, uh x = Ig x) →
      Real.sqrt (outflowNormSq P a (fun x => uh x - u x))
        + Real.sqrt a0 * nrm P.dom (fun x => uh x - u x)
        + Real.sqrt (meshSize (T i)) * nrm P.dom (fun x => ∑ k, a k * (Gh k x - G k x))
        ≤ C * meshSize (T i) ^ ((3 : ℝ) / 2) * Real.sqrt (H2NormSq P.dom u G H) := by
  sorry

end PDETextbook.LarssonThomee.Thm13_6