Heat finite-element error for nonsmooth initial data

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm10_3

/-- Ch. 10 domain: bounded convex domain in R^2 with smooth boundary (smooth defining function). -/
def IsSmoothBoundedConvexDomain (Ω : Set (Euclidean 2)) : Prop :=
  IsOpen Ω ∧ Bornology.IsBounded Ω ∧ Convex ℝ Ω ∧ Ω.Nonempty ∧
    ∃ ρ : Euclidean 2 → ℝ, (∀ n : ℕ, ContDiff ℝ n ρ) ∧ Ω = {x | ρ x < 0} ∧
      ∀ x, ρ x = 0 → fderiv ℝ ρ x ≠ 0

def triangleSet (K : Fin 3 → Euclidean 2) : Set (Euclidean 2) :=
  convexHull ℝ (Set.range K)

structure Triangulation where
  triangles : Finset (Fin 3 → Euclidean 2)
  affIndep : ∀ K ∈ triangles, AffineIndependent ℝ K

/-- Omega_h, the union of the triangles. -/
def Triangulation.region (T : Triangulation) : Set (Euclidean 2) :=
  ⋃ K ∈ T.triangles, triangleSet K

def Triangulation.nodes (T : Triangulation) : Set (Euclidean 2) :=
  ⋃ K ∈ T.triangles, Set.range K

/-- Sect. 5.2: distinct triangles meet in nothing, a common node or a common edge. -/
def Triangulation.IsConforming (T : Triangulation) : Prop :=
  ∀ K ∈ T.triangles, ∀ K' ∈ T.triangles, Set.range K ≠ Set.range K' →
    triangleSet K ∩ triangleSet K' = convexHull ℝ (Set.range K ∩ Set.range K')

/-- h = max_K diam K. -/
def Triangulation.meshSize (T : Triangulation) : ℝ :=
  ⨆ K ∈ T.triangles, Metric.diam (triangleSet K)

def Triangulation.MinAngleAtLeast (T : Triangulation) (θ₀ : ℝ) : Prop :=
  ∀ K ∈ T.triangles, ∀ i j k : Fin 3, i ≠ j → i ≠ k → j ≠ k →
    θ₀ ≤ EuclideanGeometry.angle (K j) (K i) (K k)

/-- Sect. 5.2: nodes in the closed domain, boundary vertices of Omega_h on Gamma. -/
def FitsDomain (Ω : Set (Euclidean 2)) (T : Triangulation) : Prop :=
  (∀ p ∈ T.nodes, p ∈ closure Ω) ∧
    ∀ p ∈ T.nodes, p ∈ frontier T.region → p ∈ frontier Ω

/-- Family {T_h}, 0<h<1: conforming, fitted, angles bounded below uniformly,
and Omega minus Omega_h of width O(h^2). -/
def IsAdmissibleFamily {ι : Type*} (Ω : Set (Euclidean 2)) (T : ι → Triangulation) : Prop :=
  (∀ i, (T i).IsConforming ∧ FitsDomain Ω (T i) ∧ 0 < (T i).meshSize ∧ (T i).meshSize < 1) ∧
  (∃ θ₀ : ℝ, 0 < θ₀ ∧ ∀ i, (T i).MinAngleAtLeast θ₀) ∧
  (∃ Cw : ℝ, ∀ i, ∀ x ∈ Ω \ (T i).region,
    Metric.infDist x (frontier Ω) ≤ Cw * (T i).meshSize ^ 2)

/-- Continuous piecewise linear on T, vanishing outside Omega_h (hence on its boundary). -/
def IsShFunction (T : Triangulation) (χ : Euclidean 2 → ℝ) : Prop :=
  Continuous χ ∧
    (∀ K ∈ T.triangles, ∃ A : Euclidean 2 →ᵃ[ℝ] ℝ, ∀ x ∈ triangleSet K, χ x = A x) ∧
    ∀ x, x ∉ T.region → χ x = 0

/-- S_h as a subspace of L2(Omega) (the span of an already linear set). -/
def Sh (Ω : Set (Euclidean 2)) (T : Triangulation) : Submodule ℝ (L2 ℝ Ω) :=
  Submodule.span ℝ {f : L2 ℝ Ω | ∃ χ, IsShFunction T χ ∧
    (⇑f : Euclidean 2 → ℝ) =ᵐ[domainMeasure Ω] χ}

/-- (5.38): p = P_S v. -/
def IsL2Projection {Ω : Set (Euclidean 2)} (S : Submodule ℝ (L2 ℝ Ω)) (v p : L2 ℝ Ω) : Prop :=
  p ∈ S ∧ ∀ χ ∈ S, l2Inner (p - v) χ = 0

/-- (10.3) with f = 0: (u_h,t, chi) + a(u_h, chi) = 0 for chi in S_h, t>0; u_h(0) = v_h. -/
def IsSemidiscreteSolution (Ω : Set (Euclidean 2)) (T : Triangulation)
    (vh : L2 ℝ Ω) (uh : ℝ → L2 ℝ Ω) : Prop :=
  (∀ t, 0 ≤ t → uh t ∈ Sh Ω T) ∧ ContinuousOn uh (Set.Ici 0) ∧ uh 0 = vh ∧
  ∀ t, 0 < t → ∃ duh : L2 ℝ Ω, HasDerivAt uh duh t ∧
    ∀ χ ∈ Sh Ω T, ∀ Gu Gχ : Fin 2 → Euclidean 2 → ℝ,
      HasH1Gradient Ω (⇑(uh t)) Gu → HasH1Gradient Ω (⇑χ) Gχ →
      l2Inner duh χ + (∫ x in Ω, ∑ k, Gu k x * Gχ k x) = 0

/-- (8.11): weak Dirichlet eigenpair in H_0^1. -/
def IsDirichletEigenpair (Ω : Set (Euclidean 2)) (φ : L2 ℝ Ω) (lam : ℝ) : Prop :=
  InH10Closure Ω (⇑φ) ∧ ∀ w : L2 ℝ Ω, InH10Closure Ω (⇑w) →
    ∀ Gφ Gw : Fin 2 → Euclidean 2 → ℝ,
      HasH1Gradient Ω (⇑φ) Gφ → HasH1Gradient Ω (⇑w) Gw →
      (∫ x in Ω, ∑ k, Gφ k x * Gw k x) = lam * l2Inner φ w

/-- (8.13): u(t) = E(t)v = sum_i e^{-lambda_i t} (v, phi_i) phi_i, via coefficients. -/
def IsEigenExpansionSolution (Ω : Set (Euclidean 2)) (v : L2 ℝ Ω) (u : ℝ → L2 ℝ Ω) : Prop :=
  ∀ t, 0 ≤ t → ∀ φ lam, IsDirichletEigenpair Ω φ lam →
    l2Inner (u t) φ = Real.exp (-lam * t) * l2Inner v φ

/-- Theorem 10.3 (f = 0, v_h = P_h v): ||u_h(t) - u(t)|| <= C h^2 t^{-1} ||v||, t > 0. -/
theorem thm_10_3 (Ω : Set (Euclidean 2)) (hΩ : IsSmoothBoundedConvexDomain Ω)
    {ι : Type*} (T : ι → Triangulation) (hT : IsAdmissibleFamily Ω T) :
    ∃ C : ℝ, ∀ i, ∀ v vh : L2 ℝ Ω, ∀ u uh : ℝ → L2 ℝ Ω,
      IsEigenExpansionSolution Ω v u →
      IsL2Projection (Sh Ω (T i)) v vh →
      IsSemidiscreteSolution Ω (T i) vh uh →
      ∀ t : ℝ, 0 < t →
        ‖uh t - u t‖ ≤ C * (T i).meshSize ^ 2 * t⁻¹ * ‖v‖ := by
  sorry

end PDETextbook.LarssonThomee.Thm10_3