Galerkin error for stationary transport with reaction

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm13_5

open PDETextbook

/-- Conforming triangulation of `closure Ω ⊂ ℝ²`. Triangles are nondegenerate. Any two meet in the hull of common vertices. -/
structure Triangulation (Ω : Set (Euclidean 2)) where
  tri : Finset (Fin 3 → Euclidean 2)
  nondeg : ∀ t ∈ tri, AffineIndependent ℝ t
  cover : closure Ω = ⋃ t ∈ tri, convexHull ℝ (Set.range t)
  conforming : ∀ t ∈ tri, ∀ t' ∈ tri, ∃ S : Set (Euclidean 2),
    S ⊆ Set.range t ∩ Set.range t' ∧
    convexHull ℝ (Set.range t) ∩ convexHull ℝ (Set.range t') = convexHull ℝ S

variable {Ω : Set (Euclidean 2)}

def Triangulation.IsVertex (T : Triangulation Ω) (p : Euclidean 2) : Prop :=
  ∃ t ∈ T.tri, ∃ i : Fin 3, t i = p

/-- Mesh parameter `h`: the maximal triangle diameter. -/
def Triangulation.meshSize (T : Triangulation Ω) : ℝ :=
  ⨆ t ∈ T.tri, Metric.diam (convexHull ℝ (Set.range t))

/-- `S_h`: continuous piecewise linear functions, with no boundary condition. -/
def Triangulation.InSh (T : Triangulation Ω) (χ : Euclidean 2 → ℝ) : Prop :=
  ContinuousOn χ (closure Ω) ∧
  ∀ t ∈ T.tri, ∃ A : Euclidean 2 →ᵃ[ℝ] ℝ, ∀ x ∈ convexHull ℝ (Set.range t), χ x = A x

/-- `Γ_- = {x ∈ Γ : a·n < 0}`, where `n` is an outward unit normal. -/
def inflowBoundary (Ω : Set (Euclidean 2)) (a : Euclidean 2) : Set (Euclidean 2) :=
  {x | x ∈ frontier Ω ∧ ∃ n : Euclidean 2, ‖n‖ = 1 ∧
    (∀ y ∈ Ω, inner ℝ (y - x) n ≤ 0) ∧ inner ℝ a n < 0}

/-- `S_h^- = {χ ∈ S_h : χ = 0 on Γ_-}`. -/
def Triangulation.InShMinus (T : Triangulation Ω) (a : Euclidean 2) (χ : Euclidean 2 → ℝ) : Prop :=
  T.InSh χ ∧ ∀ x ∈ inflowBoundary Ω a, χ x = 0

/-- The inflow boundary is exactly a union of triangle edges. -/
def Triangulation.MatchedInflow (T : Triangulation Ω) (a : Euclidean 2) : Prop :=
  ∃ E : Set (Euclidean 2 × Euclidean 2),
    (∀ e ∈ E, ∃ t ∈ T.tri, ∃ i j : Fin 3, i ≠ j ∧ e = (t i, t j)) ∧
    inflowBoundary Ω a = ⋃ e ∈ E, segment ℝ e.1 e.2

/-- `‖v‖`: the L² norm over Ω. -/
def l2n (Ω : Set (Euclidean 2)) (v : Euclidean 2 → ℝ) : ℝ :=
  (eLpNorm v 2 (domainMeasure Ω)).toReal

/-- `G i = ∂_i u` and `H i j = ∂_j ∂_i u`, both weak. -/
def HasH2Derivs (Ω : Set (Euclidean 2)) (u : Euclidean 2 → ℝ) (G : Fin 2 → Euclidean 2 → ℝ)
    (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ) : Prop :=
  HasH1Gradient Ω u G ∧ ∀ i, HasH1Gradient Ω (G i) (H i)

/-- `‖u‖_2 = (∑_{|α|≤2} ‖D^α u‖²)^{1/2}`, one term per multi-index. -/
def h2Norm (Ω : Set (Euclidean 2)) (u : Euclidean 2 → ℝ) (G : Fin 2 → Euclidean 2 → ℝ)
    (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (l2n Ω u ^ 2 + ∑ i, l2n Ω (G i) ^ 2 +
    (l2n Ω (H 0 0) ^ 2 + l2n Ω (H 0 1) ^ 2 + l2n Ω (H 1 1) ^ 2))

def h1Seminorm (Ω : Set (Euclidean 2)) (G : Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (∑ i, l2n Ω (G i) ^ 2)

/-- The recalled estimates (13.13), with `C` independent of `h` and `v`. Here `w = I_h v` is nodal interpolation. -/
def InterpolationEstimates (𝒯 : Set (Triangulation Ω)) : Prop :=
  ∃ C : ℝ, ∀ T ∈ 𝒯, ∀ (v : Euclidean 2 → ℝ) (G : Fin 2 → Euclidean 2 → ℝ)
    (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ) (w : Euclidean 2 → ℝ) (Gw : Fin 2 → Euclidean 2 → ℝ),
    HasH2Derivs Ω v G H → ContinuousOn v (closure Ω) → T.InSh w →
    (∀ p, T.IsVertex p → w p = v p) → HasH1Gradient Ω w Gw →
    l2n Ω (fun x => w x - v x) ≤ C * T.meshSize ^ 2 * h2Norm Ω v G H ∧
    h1Seminorm Ω (fun i x => Gw i x - G i x) ≤ C * T.meshSize * h2Norm Ω v G H

/-- (13.12): `a·∇u + a₀u = f` in Ω, and `u = g` on Γ_-. -/
def IsSolution1312 (Ω : Set (Euclidean 2)) (a : Euclidean 2) (a0 : ℝ)
    (f g u : Euclidean 2 → ℝ) (G : Fin 2 → Euclidean 2 → ℝ) : Prop :=
  HasH1Gradient Ω u G ∧ ContinuousOn u (closure Ω) ∧
  (∀ᵐ x ∂(domainMeasure Ω), (∑ i, a i * G i x) + a0 * u x = f x) ∧
  ∀ x ∈ inflowBoundary Ω a, u x = g x

/-- (13.14): `u_h ∈ S_h` and `u_h = I_h g` on Γ_- (nodal values). The Galerkin identity holds for all `χ ∈ S_h^-`. -/
def IsGalerkin1314 (T : Triangulation Ω) (a : Euclidean 2) (a0 : ℝ)
    (f g uh : Euclidean 2 → ℝ) : Prop :=
  T.InSh uh ∧ (∀ p, T.IsVertex p → p ∈ inflowBoundary Ω a → uh p = g p) ∧
  ∃ Gh : Fin 2 → Euclidean 2 → ℝ, HasH1Gradient Ω uh Gh ∧
    ∀ χ, T.InShMinus a χ →
      (∫ x in Ω, (∑ i, a i * Gh i x) * χ x) + a0 * (∫ x in Ω, uh x * χ x) =
        ∫ x in Ω, f x * χ x

/-- Theorem 13.5: `‖u_h - u‖ ≤ C h ‖u‖_2`. Ω is a bounded convex polygonal domain in ℝ², and `a`, `a₀ > 0` are constant. The family of triangulations is matched to Γ_- and satisfies (13.13). `C` is independent of `h` and of the data. -/
theorem thm_13_5 (Ω : Set (Euclidean 2)) (hopen : IsOpen Ω) (hconv : Convex ℝ Ω)
    (hbdd : Bornology.IsBounded Ω)
    (hpoly : ∃ V : Finset (Euclidean 2), closure Ω = convexHull ℝ (V : Set (Euclidean 2)))
    (a : Euclidean 2) (a0 : ℝ) (ha0 : 0 < a0)
    (𝒯 : Set (Triangulation Ω)) (hmatch : ∀ T ∈ 𝒯, T.MatchedInflow a)
    (hI : InterpolationEstimates 𝒯) :
    ∃ C : ℝ, ∀ T ∈ 𝒯, ∀ (f g u uh : Euclidean 2 → ℝ) (G : Fin 2 → Euclidean 2 → ℝ)
      (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ),
      HasH2Derivs Ω u G H → IsSolution1312 Ω a a0 f g u G → IsGalerkin1314 T a a0 f g uh →
      l2n Ω (fun x => uh x - u x) ≤ C * T.meshSize * h2Norm Ω u G H := by
  sorry

end PDETextbook.LarssonThomee.Thm13_5