Energy error for triangular finite elements

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm53

/-- A closed nondegenerate triangle: three affinely independent vertices. -/
structure Triangle where
  v : Fin 3 → Euclidean 2
  indep : AffineIndependent ℝ v

def Triangle.toSet (K : Triangle) : Set (Euclidean 2) :=
  convexHull ℝ (Set.range K.v)

/-- `h_K = diam(K)`. -/
def Triangle.diam (K : Triangle) : ℝ := Metric.diam K.toSet

/-- All interior angles of `K` are at least `θ₀`. -/
def Triangle.MinAngleGE (θ₀ : ℝ) (K : Triangle) : Prop :=
  ∀ i j k : Fin 3, i ≠ j → i ≠ k → j ≠ k →
    θ₀ ≤ EuclideanGeometry.angle (K.v j) (K.v i) (K.v k)

/-- Polygonal domain (Sect. 5.2: "a domain whose boundary `Γ` is a polygon"): an open
connected set whose frontier is a single simple closed polygon.

The polygon has `n ≥ 3` distinct vertices `p 0, …, p (n-1)`, with edges
`[p i, p (i+1 mod n)]`. Consecutive edges meet only in their shared vertex, and
non-adjacent edges are disjoint.

This excludes slit domains and multiply connected domains, whose frontier is not one
simple closed polygon. -/
def IsPolygonalDomain (Ω : Set (Euclidean 2)) : Prop :=
  IsOpen Ω ∧ IsConnected Ω ∧
    ∃ (n : ℕ) (p : Fin n → Euclidean 2), 3 ≤ n ∧ Function.Injective p ∧
      (∀ i : Fin n,
        segment ℝ (p i) (p (finRotate n i)) ∩
          segment ℝ (p (finRotate n i)) (p (finRotate n (finRotate n i))) =
          {p (finRotate n i)}) ∧
      (∀ i j : Fin n, i ≠ j → j ≠ finRotate n i → i ≠ finRotate n j →
        Disjoint (segment ℝ (p i) (p (finRotate n i)))
          (segment ℝ (p j) (p (finRotate n j)))) ∧
      frontier Ω = ⋃ i : Fin n, segment ℝ (p i) (p (finRotate n i))

/-- Triangulation `𝒯_h` of `Ω` as in Sect. 5.2. -/
structure Triangulation (Ω : Set (Euclidean 2)) where
  tris : Finset Triangle
  cover : closure Ω = ⋃ K ∈ tris, K.toSet
  inter : ∀ K ∈ tris, ∀ K' ∈ tris, K ≠ K' →
    K.toSet ∩ K'.toSet = ∅ ∨
    (∃ i j : Fin 3, K.v i = K'.v j ∧ K.toSet ∩ K'.toSet = {K.v i}) ∨
    (∃ i₁ i₂ j₁ j₂ : Fin 3, i₁ ≠ i₂ ∧ K.v i₁ = K'.v j₁ ∧ K.v i₂ = K'.v j₂ ∧
      K.toSet ∩ K'.toSet = segment ℝ (K.v i₁) (K.v i₂))
  noHangingNode : ∀ K ∈ tris, ∀ K' ∈ tris, ∀ i j₁ j₂ : Fin 3, j₁ ≠ j₂ →
    K.v i ∉ openSegment ℝ (K'.v j₁) (K'.v j₂)

/-- `h = max_K h_K` (diameters are nonnegative, so the iSup is the maximum). -/
def Triangulation.meshSize {Ω : Set (Euclidean 2)} (T : Triangulation Ω) : ℝ :=
  ⨆ K ∈ T.tris, K.diam

/-- `S_h`: continuous on `Ω̄`, linear (affine) on each `K`, zero on `Γ`. -/
def InSh {Ω : Set (Euclidean 2)} (T : Triangulation Ω) (v : Euclidean 2 → ℝ) : Prop :=
  ContinuousOn v (closure Ω) ∧
  (∀ K ∈ T.tris, ∃ A : Euclidean 2 →ᵃ[ℝ] ℝ, ∀ x ∈ K.toSet, v x = A x) ∧
  ∀ x ∈ frontier Ω, v x = 0

/-- `a(v,w) = ∫_Ω a ∇v·∇w`, in terms of (weak) gradients `G`, `H`. -/
def aForm (Ω : Set (Euclidean 2)) (a : Euclidean 2 → ℝ)
    (G H : Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  ∫ x in Ω, a x * ∑ i, G i x * H i x

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

/-- Energy norm `‖v‖_a = a(v,v)^{1/2}` from a gradient `G` of `v`. -/
def energyNorm (Ω : Set (Euclidean 2)) (a : Euclidean 2 → ℝ)
    (G : Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (aForm Ω a G G)

/-- `|v|_1 = ‖∇v‖` from a gradient `G` of `v`. -/
def H1Seminorm (Ω : Set (Euclidean 2)) (G : Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (∫ x in Ω, ∑ i, G i x ^ 2)

/-- `‖u‖_2 = (∑_{|α|≤2} ‖D^α u‖²)^{1/2}`. Here `H i j` is the weak `∂_j` of `G i`.
Each multi-index of order two is counted once (pairs `i ≤ j`). -/
def H2Norm (Ω : Set (Euclidean 2)) (u : Euclidean 2 → ℝ)
    (G : Fin 2 → Euclidean 2 → ℝ) (H : Fin 2 → Fin 2 → Euclidean 2 → ℝ) : ℝ :=
  Real.sqrt (∫ x in Ω, (u x ^ 2 + ∑ i, G i x ^ 2 +
    ∑ i, ∑ j, if i ≤ j then H i j x ^ 2 else 0))

/-- (5.24): `u ∈ H_0^1` with `a(u,v) = (f,v)` for all `v ∈ H_0^1`. -/
def IsVariationalSolution (Ω : Set (Euclidean 2)) (a f u : Euclidean 2 → ℝ) : Prop :=
  InH10Closure Ω u ∧
  ∀ Gu, HasH1Gradient Ω u Gu → ∀ v : Euclidean 2 → ℝ, InH10Closure Ω v →
    ∀ Gv, HasH1Gradient Ω v Gv → aForm Ω a Gu Gv = l2Pair Ω f v

/-- (5.26): `u_h ∈ S_h` with `a(u_h,χ) = (f,χ)` for all `χ ∈ S_h`. -/
def IsFEMSolution {Ω : Set (Euclidean 2)} (T : Triangulation Ω)
    (a f uh : Euclidean 2 → ℝ) : Prop :=
  InSh T uh ∧
  ∀ Guh, HasH1Gradient Ω uh Guh → ∀ χ : Euclidean 2 → ℝ, InSh T χ →
    ∀ Gχ, HasH1Gradient Ω χ Gχ → aForm Ω a Guh Gχ = l2Pair Ω f χ

/-- Theorem 5.3 (Larsson–Thomée).

Setting of Sect. 5.2: `Ω ⊂ ℝ²` is a polygonal domain (its boundary is a simple closed
polygon), `a` is smooth with `a ≥ a₀ > 0` on `Ω̄`, and `f ∈ L₂`. The triangulations form a
family with mesh sizes in `(0,1)` whose angles are bounded below by `θ₀ > 0` independently
of `h`.

Conclusions:
* (5.40) `‖u_h-u‖_a = min_{χ∈S_h} ‖χ-u‖_a`;
* (5.41) `|u_h-u|_1 ≤ C h ‖u‖_2`, with `C` independent of `h`, `f`, `u`.

(5.41) is asserted whenever `u ∈ H²`; otherwise its right side is infinite. -/
theorem thm_5_3 (Ω : Set (Euclidean 2)) (hΩ : IsPolygonalDomain Ω)
    (a : Euclidean 2 → ℝ) (ha_smooth : ∀ n : ℕ, ContDiffOn ℝ n a (closure Ω))
    (a₀ : ℝ) (ha₀ : 0 < a₀) (ha_lower : ∀ x ∈ closure Ω, a₀ ≤ a x)
    (𝒯 : Set (Triangulation Ω)) (θ₀ : ℝ) (hθ₀ : 0 < θ₀)
    (hangle : ∀ T ∈ 𝒯, ∀ K ∈ T.tris, K.MinAngleGE θ₀)
    (hmesh : ∀ T ∈ 𝒯, T.meshSize ∈ Set.Ioo (0 : ℝ) 1) :
    ∃ C : ℝ, ∀ T ∈ 𝒯, ∀ f : Euclidean 2 → ℝ, InL2 Ω f →
      ∀ u : Euclidean 2 → ℝ, IsVariationalSolution Ω a f u →
      ∀ uh : Euclidean 2 → ℝ, IsFEMSolution T a f uh →
        (∀ Gu Guh, HasH1Gradient Ω u Gu → HasH1Gradient Ω uh Guh →
          IsLeast {r : ℝ | ∃ χ : Euclidean 2 → ℝ, ∃ Gχ : Fin 2 → Euclidean 2 → ℝ,
              InSh T χ ∧ HasH1Gradient Ω χ Gχ ∧
              r = energyNorm Ω a (fun i x => Gχ i x - Gu i x)}
            (energyNorm Ω a (fun i x => Guh i x - Gu i x))) ∧
        (∀ (Gu : Fin 2 → Euclidean 2 → ℝ) (Hu : Fin 2 → Fin 2 → Euclidean 2 → ℝ) Guh,
          HasH1Gradient Ω u Gu → (∀ i, HasH1Gradient Ω (Gu i) (Hu i)) →
          HasH1Gradient Ω uh Guh →
          H1Seminorm Ω (fun i x => Guh i x - Gu i x) ≤
            C * T.meshSize * H2Norm Ω u Gu Hu) := by
  sorry

end PDETextbook.LarssonThomee.Thm53