Finite-element error with numerical quadrature

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm57

abbrev E2 := Euclidean 2

/-- Polygonal domain (Sect. 5.2): "a domain whose boundary Γ is a polygon".
Ω is open, connected and bounded. Its boundary `frontier Ω` is a simple closed polygon.
The polygon has `m ≥ 3` pairwise distinct vertices `p 0, …, p (m-1)`, and its edges join
`p k` to `p (k+1 mod m)`. Two distinct edges meet only in common endpoints. So every edge
is nondegenerate, and the polygon has no slits or isolated points. -/
def PolygonalDomain (Ω : Set E2) : Prop :=
  IsOpen Ω ∧ IsConnected Ω ∧ Bornology.IsBounded Ω ∧
    ∃ (m : ℕ) (p : Fin m → E2), 3 ≤ m ∧ Function.Injective p ∧
      (∀ k l : Fin m, k ≠ l →
        segment ℝ (p k) (p (finRotate m k)) ∩ segment ℝ (p l) (p (finRotate m l)) ⊆
          ({p k, p (finRotate m k)} : Set E2) ∩ {p l, p (finRotate m l)}) ∧
      frontier Ω = ⋃ k : Fin m, segment ℝ (p k) (p (finRotate m k))

/-- A nondegenerate triangle with vertices `P_{l,K}`. -/
structure Triangle where
  v : Fin 3 → E2
  indep : AffineIndependent ℝ v

def Triangle.set (K : Triangle) : Set E2 := convexHull ℝ (Set.range K.v)

/-- `|K| = area(K)`. -/
def Triangle.area (K : Triangle) : ℝ := (volume K.set).toReal

/-- Barycenter `P_K = (1/3) Σ_l P_{l,K}`. -/
def Triangle.bary (K : Triangle) : E2 := (1 / 3 : ℝ) • ∑ l, K.v l

/-- Conforming triangulation of `closure Ω`. -/
structure Triangulation (Ω : Set E2) where
  n : ℕ
  K : Fin n → Triangle
  cover : closure Ω = ⋃ i, (K i).set
  conforming : ∀ i j, i ≠ j →
    Set.range (K i).v ≠ Set.range (K j).v ∧
    (K i).set ∩ (K j).set = convexHull ℝ (Set.range (K i).v ∩ Set.range (K j).v)

/-- `h = max_K diam K`. -/
def Triangulation.mesh {Ω : Set E2} (T : Triangulation Ω) : ℝ :=
  ⨆ i : Fin T.n, Metric.diam (T.K i).set

/-- All angles of all triangles are at least `θ0`. -/
def Triangulation.AnglesGE {Ω : Set E2} (T : Triangulation Ω) (θ0 : ℝ) : Prop :=
  ∀ i, ∀ l j k : Fin 3, l ≠ j → l ≠ k → j ≠ k →
    θ0 ≤ EuclideanGeometry.angle ((T.K i).v j) ((T.K i).v l) ((T.K i).v k)

/-- `S_h`: continuous on `closure Ω`, linear on each `K`, zero on `Γ`. -/
def InSh {Ω : Set E2} (T : Triangulation Ω) (v : E2 → ℝ) : Prop :=
  ContinuousOn v (closure Ω) ∧
  (∀ i, ∃ (c : ℝ) (g : E2), ∀ x ∈ (T.K i).set, v x = c + inner ℝ g x) ∧
  ∀ x ∈ frontier Ω, v x = 0

/-- The two quadrature rules (5.60) and (5.64). -/
inductive QuadRule
  | barycentric
  | nodal

def qK : QuadRule → Triangle → (E2 → ℝ) → ℝ
  | .barycentric, K, φ => K.area * φ K.bary
  | .nodal, K, φ => (1 / 3 : ℝ) * K.area * ∑ l, φ (K.v l)

/-- `∇v·∇w` on `K` for `v, w` linear on `K` (gradient taken at the interior point `P_K`). -/
def elemGradDot (K : Triangle) (v w : E2 → ℝ) : ℝ :=
  ∑ j, partialDeriv v j K.bary * partialDeriv w j K.bary

/-- `a_h(v,w)` of (5.59). -/
def ah {Ω : Set E2} (q : QuadRule) (T : Triangulation Ω) (a v w : E2 → ℝ) : ℝ :=
  ∑ i, qK q (T.K i) (fun x => a x * elemGradDot (T.K i) v w)

/-- `(f,w)_h` of (5.59). -/
def fh {Ω : Set E2} (q : QuadRule) (T : Triangulation Ω) (f w : E2 → ℝ) : ℝ :=
  ∑ i, qK q (T.K i) (fun x => f x * w x)

/-- `G`, `H` are first and second weak derivatives of `u` in `L₂`. -/
def H2Data (Ω : Set E2) (u : E2 → ℝ) (G : Fin 2 → E2 → ℝ) (H : Fin 2 → Fin 2 → E2 → ℝ) : Prop :=
  HasH1Gradient Ω u G ∧ ∀ i j, HasWeakPartial Ω (G i) j (H i j) ∧ InL2 Ω (H i j)

/-- `‖u‖_2`. -/
def H2Norm (Ω : Set E2) (u : E2 → ℝ) (G : Fin 2 → E2 → ℝ) (H : Fin 2 → Fin 2 → E2 → ℝ) : ℝ :=
  Real.sqrt ((∫ x in Ω, u x ^ 2) + (∑ i, ∫ x in Ω, G i x ^ 2) +
    ∑ i, ∑ j, ∫ x in Ω, H i j x ^ 2)

/-- `‖a‖_{C²}` over `closure Ω`. -/
def C2Norm (Ω : Set E2) (a : E2 → ℝ) : ℝ :=
  ⨆ x : closure Ω, ⨆ k : Fin 3, ‖iteratedFDerivWithin ℝ (k : ℕ) a (closure Ω) x‖

/-- `|u_h - u|_1`, computed elementwise. -/
def H1SemiErr {Ω : Set E2} (T : Triangulation Ω) (uh : E2 → ℝ) (G : Fin 2 → E2 → ℝ) : ℝ :=
  Real.sqrt (∑ i, ∫ x in (T.K i).set,
    ∑ j, (partialDeriv uh j (T.K i).bary - G j x) ^ 2)

/-- Solution of (5.24). -/
def IsWeakSolution (Ω : Set E2) (a f u : E2 → ℝ) (G : Fin 2 → E2 → ℝ) : Prop :=
  InH10Closure Ω u ∧ HasH1Gradient Ω u G ∧
  ∀ (v : E2 → ℝ) (Gv : Fin 2 → E2 → ℝ), InH10Closure Ω v → HasH1Gradient Ω v Gv →
    (∫ x in Ω, a x * ∑ i, G i x * Gv i x) = ∫ x in Ω, f x * v x

/-- Solution of (5.58). -/
def IsQuadSolution {Ω : Set E2} (q : QuadRule) (T : Triangulation Ω) (a f uh : E2 → ℝ) : Prop :=
  InSh T uh ∧ ∀ χ, InSh T χ → ah q T a uh χ = fh q T f χ

/-- Theorem 5.7 (5.69). No source warnings.
Standing assumptions of Sect. 5.2 are kept: Ω is polygonal (its boundary is a simple closed
polygon), and `a` is smooth on `closure Ω` with `a ≥ a0 > 0`. The source assumes `f ∈ L_2`.
The right-hand side of (5.69) contains `‖f‖_2`, so `f` is taken with H² data, which sharpens
that standing assumption. The quadrature (5.59) evaluates `f` at points, so we use a
representative of `f` that is continuous on `closure Ω`.
The triangulation family `{T_h}_{0<h<1}` has `h = max_K h_K` and angles bounded below
independently of `h`. The constant `C` is chosen before `h`, `f`, `u` and `u_h`. -/
theorem thm_5_7 (Ω : Set E2) (hΩ : PolygonalDomain Ω)
    (a : E2 → ℝ) (ha : ContDiffOn ℝ (⊤ : ℕ∞) a (closure Ω))
    (ha0 : ∃ a0 : ℝ, 0 < a0 ∧ ∀ x ∈ closure Ω, a0 ≤ a x)
    (q : QuadRule)
    (T : Set.Ioo (0 : ℝ) 1 → Triangulation Ω)
    (hmesh : ∀ h, (T h).mesh = h)
    (hang : ∃ θ0 : ℝ, 0 < θ0 ∧ ∀ h, (T h).AnglesGE θ0) :
    ∃ C : ℝ, ∀ h : Set.Ioo (0 : ℝ) 1,
      ∀ (f : E2 → ℝ) (Gf : Fin 2 → E2 → ℝ) (Hf : Fin 2 → Fin 2 → E2 → ℝ)
        (u : E2 → ℝ) (G : Fin 2 → E2 → ℝ) (Hu : Fin 2 → Fin 2 → E2 → ℝ) (uh : E2 → ℝ),
      ContinuousOn f (closure Ω) → H2Data Ω f Gf Hf →
      IsWeakSolution Ω a f u G → H2Data Ω u G Hu →
      IsQuadSolution q (T h) a f uh →
      H1SemiErr (T h) uh G ≤
        C * (h : ℝ) * H2Norm Ω u G Hu +
          C * (h : ℝ) ^ 2 * (C2Norm Ω a * H2Norm Ω u G Hu + H2Norm Ω f Gf Hf) := by
  sorry

end PDETextbook.LarssonThomee.Thm57