Coercivity and consistency of finite-element quadrature

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Lem51

/-- Closed triangle with vertices v 0, v 1, v 2. -/
def triangle (v : Fin 3 → Euclidean 2) : Set (Euclidean 2) :=
  convexHull ℝ (Set.range v)

/-- Polygonal domain: a bounded connected open set whose boundary Γ is a finite union of segments. -/
def IsPolygonalDomain (Ω : Set (Euclidean 2)) : Prop :=
  IsOpen Ω ∧ IsConnected Ω ∧ Bornology.IsBounded Ω ∧
    ∃ (m : ℕ) (A B : Fin m → Euclidean 2), frontier Ω = ⋃ k, segment ℝ (A k) (B k)

/-- Triangulation (Sect. 5.2): closed triangles covering the closure of Ω. Distinct triangles meet
in the empty set, a common node or a common edge. No node lies in the interior of an edge. -/
structure Triangulation (Ω : Set (Euclidean 2)) where
  n : ℕ
  vert : Fin n → Fin 3 → Euclidean 2
  indep : ∀ K, AffineIndependent ℝ (vert K)
  cover : closure Ω = ⋃ K, triangle (vert K)
  conforming : ∀ K L, K ≠ L →
    (triangle (vert K) ∩ triangle (vert L) = ∅ ∨
     (∃ p q, vert K p = vert L q ∧ triangle (vert K) ∩ triangle (vert L) = {vert K p}) ∨
     (∃ p p' q q', p ≠ p' ∧ vert K p = vert L q ∧ vert K p' = vert L q' ∧
        triangle (vert K) ∩ triangle (vert L) = segment ℝ (vert K p) (vert K p')))
  noHangingNode : ∀ K L (p p' q : Fin 3), p ≠ p' →
    vert L q ∉ openSegment ℝ (vert K p) (vert K p')

/-- h = max_K h_K, h_K = diam K. -/
def meshSize {Ω : Set (Euclidean 2)} (T : Triangulation Ω) : ℝ :=
  ⨆ K : Fin T.n, Metric.diam (triangle (T.vert K))

/-- Angles of all triangles of all members of the family are bounded below independently of h. -/
def MinAngleBounded {Ω : Set (Euclidean 2)} (𝒯 : ℝ → Triangulation Ω) : Prop :=
  ∃ θ₀ : ℝ, 0 < θ₀ ∧ ∀ h : ℝ, 0 < h → h < 1 → ∀ K : Fin (𝒯 h).n, ∀ p q r : Fin 3,
    q ≠ p → r ≠ p → q ≠ r →
      θ₀ ≤ EuclideanGeometry.angle ((𝒯 h).vert K q) ((𝒯 h).vert K p) ((𝒯 h).vert K r)

/-- S_h: continuous on the closure of Ω, linear on each K, vanishing on Γ. -/
def InSh {Ω : Set (Euclidean 2)} (T : Triangulation Ω) (v : Euclidean 2 → ℝ) : Prop :=
  ContinuousOn v (closure Ω) ∧
  (∀ K, ∃ c : ℝ, ∃ ℓ : Euclidean 2 →L[ℝ] ℝ, ∀ x ∈ triangle (T.vert K), v x = c + ℓ x) ∧
  ∀ x ∈ frontier Ω, v x = 0

/-- |K| = area(K). -/
def area (v : Fin 3 → Euclidean 2) : ℝ := (volume (triangle v)).toReal

/-- Barycenter P_K. -/
def barycenter (v : Fin 3 → Euclidean 2) : Euclidean 2 := (1 / 3 : ℝ) • ∑ l, v l

inductive QuadRule
  | barycentric
  | nodal

/-- q_K from (5.60) (barycentric) and (5.64) (nodal). -/
def quadK : QuadRule → (Fin 3 → Euclidean 2) → (Euclidean 2 → ℝ) → ℝ
  | .barycentric, v, φ => area v * φ (barycenter v)
  | .nodal, v, φ => (1 / 3 : ℝ) * area v * ∑ l, φ (v l)

/-- i-th component of the gradient of v restricted to K. This gradient is constant on K and
is taken at the interior point P_K. -/
def elemGrad {Ω : Set (Euclidean 2)} (T : Triangulation Ω) (K : Fin T.n)
    (v : Euclidean 2 → ℝ) (i : Fin 2) : ℝ :=
  partialDeriv v i (barycenter (T.vert K))

/-- a_h(v,w) = Σ_K q_K(a ∇v·∇w), (5.59). -/
def ah {Ω : Set (Euclidean 2)} (rule : QuadRule) (T : Triangulation Ω)
    (a v w : Euclidean 2 → ℝ) : ℝ :=
  ∑ K, quadK rule (T.vert K) (fun x => a x * ∑ i, elemGrad T K v i * elemGrad T K w i)

/-- (f,w)_h = Σ_K q_K(f w), (5.59). -/
def lh {Ω : Set (Euclidean 2)} (rule : QuadRule) (T : Triangulation Ω)
    (f w : Euclidean 2 → ℝ) : ℝ :=
  ∑ K, quadK rule (T.vert K) (fun x => f x * w x)

/-- a(v,w) = ∫_Ω a ∇v·∇w, with weak gradients Gv and Gw. -/
def aForm (Ω : Set (Euclidean 2)) (a : Euclidean 2 → ℝ) (Gv Gw : Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  ∫ x in Ω, a x * ∑ i, Gv i x * Gw i x

/-- (f,w) = ∫_Ω f w. -/
def l2Form (Ω : Set (Euclidean 2)) (f w : Euclidean 2 → ℝ) : ℝ := ∫ x in Ω, f x * w x

/-- |v|_1 from the weak gradient G. -/
def H1Semi (Ω : Set (Euclidean 2)) (G : Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (∑ i, ∫ x in Ω, (G i x) ^ 2)

/-- Weak first derivatives G i = D_i f and second derivatives H i j = D_j D_i f. -/
def HasH2Derivs (Ω : Set (Euclidean 2)) (f : Euclidean 2 → ℝ)
    (G : Fin 2 → Euclidean 2 → ℝ) (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ) : Prop :=
  HasH1Gradient Ω f G ∧ ∀ i, HasH1Gradient Ω (G i) (H i)

/-- ‖f‖_2 = (Σ_{|α|≤2} ‖D^α f‖²)^{1/2}. The multi-indices of order 2 are the pairs i ≤ j. -/
def H2Norm (Ω : Set (Euclidean 2)) (f : Euclidean 2 → ℝ)
    (G : Fin 2 → Euclidean 2 → ℝ) (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt ((∫ x in Ω, (f x) ^ 2) + (∑ i, ∫ x in Ω, (G i x) ^ 2) +
    ∑ i, ∑ j, if i ≤ j then ∫ x in Ω, (H i j x) ^ 2 else 0)

/-- ‖a‖_{C²}: sup over the closure of Ω of |D^α a| for |α| ≤ 2. -/
def C2Norm (Ω : Set (Euclidean 2)) (a : Euclidean 2 → ℝ) : ℝ :=
  ⨆ x : closure Ω, max |a x| (max (⨆ i : Fin 2, |partialDeriv a i x|)
    (⨆ i : Fin 2, ⨆ j : Fin 2, |partialDeriv (fun y => partialDeriv a j y) i x|))

/-- Lemma 5.1 (Larsson–Thomée). The quadrature in (5.59) is (5.60) or (5.64), selected by
`rule`. The coefficient a is smooth with a ≥ a₀ > 0 on the closure of Ω (Sect. 5.2 convention).
All constants are independent of h in the family {𝒯_h}_{0<h<1}, whose angles are bounded below.
* (5.65): the upper constant may depend on a, because the source proof uses a(x) ≤ C.
* (5.66): a single constant serves every admissible coefficient a. The bound shows its
  dependence on a only through the explicit factor ‖a‖_{C²}.
* (5.67): a single constant serves every f. -/
theorem lem_5_1 (Ω : Set (Euclidean 2)) (hΩ : IsPolygonalDomain Ω)
    (𝒯 : ℝ → Triangulation Ω) (h𝒯 : ∀ h : ℝ, 0 < h → h < 1 → meshSize (𝒯 h) = h)
    (hangle : MinAngleBounded 𝒯) (rule : QuadRule) :
    (∀ (a : Euclidean 2 → ℝ), ContDiff ℝ (⊤ : ℕ∞) a → ∀ a₀ : ℝ, 0 < a₀ →
      (∀ x ∈ closure Ω, a₀ ≤ a x) →
      ∃ C : ℝ, ∀ h : ℝ, 0 < h → h < 1 → ∀ (χ : Euclidean 2 → ℝ) (Gχ : Fin 2 → Euclidean 2 → ℝ),
        InSh (𝒯 h) χ → HasH1Gradient Ω χ Gχ →
          a₀ * H1Semi Ω Gχ ^ 2 ≤ ah rule (𝒯 h) a χ χ ∧
          ah rule (𝒯 h) a χ χ ≤ C * H1Semi Ω Gχ ^ 2) ∧
    (∃ C : ℝ, ∀ (a : Euclidean 2 → ℝ), ContDiff ℝ (⊤ : ℕ∞) a → ∀ a₀ : ℝ, 0 < a₀ →
      (∀ x ∈ closure Ω, a₀ ≤ a x) →
      ∀ h : ℝ, 0 < h → h < 1 → ∀ (ψ χ : Euclidean 2 → ℝ) (Gψ Gχ : Fin 2 → Euclidean 2 → ℝ),
        InSh (𝒯 h) ψ → HasH1Gradient Ω ψ Gψ → InSh (𝒯 h) χ → HasH1Gradient Ω χ Gχ →
          |ah rule (𝒯 h) a ψ χ - aForm Ω a Gψ Gχ| ≤
            C * h ^ 2 * C2Norm Ω a * H1Semi Ω Gψ * H1Semi Ω Gχ) ∧
    (∃ C : ℝ, ∀ (f : Euclidean 2 → ℝ) (G : Fin 2 → Euclidean 2 → ℝ)
        (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ),
        ContinuousOn f (closure Ω) → HasH2Derivs Ω f G H →
        ∀ h : ℝ, 0 < h → h < 1 → ∀ (χ : Euclidean 2 → ℝ) (Gχ : Fin 2 → Euclidean 2 → ℝ),
          InSh (𝒯 h) χ → HasH1Gradient Ω χ Gχ →
            |lh rule (𝒯 h) f χ - l2Form Ω f χ| ≤ C * h ^ 2 * H2Norm Ω f G H * H1Semi Ω Gχ) := by
  sorry

end PDETextbook.LarssonThomee.Lem51