import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm51

/-- The coordinate of a point of the one-dimensional Euclidean space. -/
def coord (x : Euclidean 1) : ℝ := x 0

/-- `Ω = (0,1)`, realized inside `Euclidean 1`. -/
def Omega : Set (Euclidean 1) := {x | 0 < coord x ∧ coord x < 1}

/-- A function on `ℝ` (only its values on `Ω̄ = [0,1]` matter) viewed on `Euclidean 1`. -/
def liftFn (v : ℝ → ℝ) : Euclidean 1 → ℝ := fun x => v (coord x)

/-- A partition `0 = x_0 < x_1 < ⋯ < x_M = 1` of `Ω`. -/
structure Partition where
  M : ℕ
  x : Fin (M + 1) → ℝ
  strictMono : StrictMono x
  first : x 0 = 0
  last : x (Fin.last M) = 1

/-- `h = max_j h_j`, where `h_j = x_j - x_{j-1}`. -/
def meshSize (P : Partition) : ℝ :=
  ⨆ j : Fin P.M, (P.x j.succ - P.x j.castSucc)

/-- `S_h = { v ∈ C(Ω̄) : v linear on each K_j = [x_{j-1}, x_j], v(0) = v(1) = 0 }`. -/
def InSh (P : Partition) (v : ℝ → ℝ) : Prop :=
  ContinuousOn v (Set.Icc 0 1) ∧
  (∀ j : Fin P.M, ∃ α β : ℝ, ∀ t ∈ Set.Icc (P.x j.castSucc) (P.x j.succ), v t = α * t + β) ∧
  v 0 = 0 ∧ v 1 = 0

/-- `a(v,w) = ∫_Ω (a v' w' + c v w) dx`, with `v', w'` the given weak derivatives. -/
def bilinA (a c : ℝ → ℝ) (v v' w w' : Euclidean 1 → ℝ) : ℝ :=
  ∫ x in Omega, (a (coord x) * v' x * w' x + c (coord x) * v x * w x)

/-- Energy norm (5.11): `‖v‖_a = a(v,v)^{1/2}`. -/
def energyNorm (a c : ℝ → ℝ) (v v' : Euclidean 1 → ℝ) : ℝ :=
  Real.sqrt (bilinA a c v v' v v')

/-- The `L_2(Ω)` norm `‖·‖`. -/
def l2Norm (v : Euclidean 1 → ℝ) : ℝ :=
  (eLpNorm v 2 (domainMeasure Omega)).toReal

/-- `‖v‖_2 = (‖v‖² + ‖v'‖² + ‖v''‖²)^{1/2}` (the `H²(Ω)` norm in one dimension). -/
def h2Norm (v v' v'' : Euclidean 1 → ℝ) : ℝ :=
  Real.sqrt (l2Norm v ^ 2 + l2Norm v' ^ 2 + l2Norm v'' ^ 2)

/-- (5.2): `u ∈ H_0^1` with weak derivative `u1`, and `a(u,φ) = (f,φ)` for all `φ ∈ H_0^1`. -/
def IsVariationalSolution (a c : ℝ → ℝ) (f u u1 : Euclidean 1 → ℝ) : Prop :=
  InH10Closure Omega u ∧ HasH1Gradient Omega u (fun _ => u1) ∧
  ∀ φ φ1 : Euclidean 1 → ℝ, InH10Closure Omega φ → HasH1Gradient Omega φ (fun _ => φ1) →
    bilinA a c u u1 φ φ1 = ∫ x in Omega, f x * φ x

/-- (5.3): `u_h ∈ S_h` with weak derivative `uh1`, and `a(u_h,χ) = (f,χ)` for all `χ ∈ S_h`. -/
def IsGalerkinSolution (a c : ℝ → ℝ) (P : Partition) (f : Euclidean 1 → ℝ)
    (uh : ℝ → ℝ) (uh1 : Euclidean 1 → ℝ) : Prop :=
  InSh P uh ∧ HasH1Gradient Omega (liftFn uh) (fun _ => uh1) ∧
  ∀ (χ : ℝ → ℝ) (χ1 : Euclidean 1 → ℝ), InSh P χ →
    HasH1Gradient Omega (liftFn χ) (fun _ => χ1) →
    bilinA a c (liftFn uh) uh1 (liftFn χ) χ1 = ∫ x in Omega, f x * liftFn χ x

/-- Theorem 5.1 (Larsson–Thomée). For `-(au')' + cu = f` in `Ω = (0,1)`, `u(0)=u(1)=0`,
with `a, c` smooth on `Ω̄`, `a ≥ a₀ > 0`, `c ≥ 0`, `f ∈ L₂`: if `u` solves (5.2)
(`u ∈ H²`, `u1 = u'`, `u2 = u''`) and `u_h` solves (5.3), then
(5.12) `‖u_h - u‖_a = min_{χ ∈ S_h} ‖χ - u‖_a` and (5.13) `‖u_h' - u'‖ ≤ C h ‖u‖_2`,
with `C` independent of the partition, `h`, `f` and `u`. -/
theorem thm_5_1 (a c : ℝ → ℝ)
    (ha_smooth : ContDiffOn ℝ (⊤ : ℕ∞) a (Set.Icc 0 1))
    (hc_smooth : ContDiffOn ℝ (⊤ : ℕ∞) c (Set.Icc 0 1))
    (ha_pos : ∃ a₀ : ℝ, 0 < a₀ ∧ ∀ t ∈ Set.Icc (0 : ℝ) 1, a₀ ≤ a t)
    (hc_nonneg : ∀ t ∈ Set.Icc (0 : ℝ) 1, 0 ≤ c t) :
    (∀ (P : Partition) (f u u1 u2 : Euclidean 1 → ℝ) (uh : ℝ → ℝ) (uh1 : Euclidean 1 → ℝ),
      InL2 Omega f →
      IsVariationalSolution a c f u u1 →
      HasH1Gradient Omega u1 (fun _ => u2) →
      IsGalerkinSolution a c P f uh uh1 →
      IsLeast
        {r : ℝ | ∃ (χ : ℝ → ℝ) (χ1 : Euclidean 1 → ℝ), InSh P χ ∧
            HasH1Gradient Omega (liftFn χ) (fun _ => χ1) ∧
            r = energyNorm a c (fun x => liftFn χ x - u x) (fun x => χ1 x - u1 x)}
        (energyNorm a c (fun x => liftFn uh x - u x) (fun x => uh1 x - u1 x))) ∧
    (∃ C : ℝ, ∀ (P : Partition) (f u u1 u2 : Euclidean 1 → ℝ) (uh : ℝ → ℝ)
        (uh1 : Euclidean 1 → ℝ),
      InL2 Omega f →
      IsVariationalSolution a c f u u1 →
      HasH1Gradient Omega u1 (fun _ => u2) →
      IsGalerkinSolution a c P f uh uh1 →
      l2Norm (fun x => uh1 x - u1 x) ≤ C * meshSize P * h2Norm u u1 u2) := by
  sorry

end PDETextbook.LarssonThomee.Thm51