Energy minimization for an elliptic boundary problem

import PDETextbook.Base

open MeasureTheory
open scoped BigOperators ENNReal

noncomputable section

namespace PDETextbook.LarssonThomee.Thm37

/-- Smooth boundary: near each boundary point, Ω is the sublevel set `{ρ < 0}` of a
C^∞ defining function `ρ` with nonvanishing differential at that point. -/
def SmoothBoundary {d : ℕ} (Ω : Set (Euclidean d)) : Prop :=
  ∀ p ∈ frontier Ω, ∃ U ∈ nhds p, ∃ ρ : Euclidean d → ℝ,
    ContDiff ℝ (⊤ : ℕ∞) ρ ∧ fderiv ℝ ρ p ≠ 0 ∧ Ω ∩ U = {x ∈ U | ρ x < 0}

/-- Polygonal (d = 2) / polyhedral (d = 3) domain: the closure of Ω is a finite union of
convex polytopes (convex hulls of finitely many points), and Ω lies on one side of its
boundary, i.e. `Ω = interior (closure Ω)`. The second condition excludes punctured or
slit domains, whose boundary is not a polygon/polyhedron surface. -/
def PolyhedralBoundary {d : ℕ} (Ω : Set (Euclidean d)) : Prop :=
  interior (closure Ω) = Ω ∧
    ∃ P : Finset (Finset (Euclidean d)),
      closure Ω = ⋃ S ∈ P, convexHull ℝ (S : Set (Euclidean d))

/-- The book's standing domain convention: Ω ⊂ ℝ^d, d ∈ {1,2,3}, is a bounded domain
(open and connected) whose boundary Γ is smooth, or a polygon (d = 2) / polyhedron (d = 3). -/
def StandingDomain {d : ℕ} (Ω : Set (Euclidean d)) : Prop :=
  (d = 1 ∨ d = 2 ∨ d = 3) ∧ IsOpen Ω ∧ IsConnected Ω ∧ Bornology.IsBounded Ω ∧
    (SmoothBoundary Ω ∨ ((d = 2 ∨ d = 3) ∧ PolyhedralBoundary Ω))

/-- Condition (3.2) together with the smoothness of the coefficients `a, b, c`
(smooth on the closure of Ω): `a(x) ≥ a₀ > 0`, `c(x) ≥ 0` for all `x ∈ Ω`. -/
def CoefficientCondition {d : ℕ} (Ω : Set (Euclidean d))
    (a : Euclidean d → ℝ) (b : Euclidean d → Euclidean d) (c : Euclidean d → ℝ) : Prop :=
  ContDiffOn ℝ (⊤ : ℕ∞) a (closure Ω) ∧ ContDiffOn ℝ (⊤ : ℕ∞) b (closure Ω) ∧
  ContDiffOn ℝ (⊤ : ℕ∞) c (closure Ω) ∧
  ∃ a₀ : ℝ, 0 < a₀ ∧ ∀ x ∈ Ω, a₀ ≤ a x ∧ 0 ≤ c x

/-- `u ∈ H₀¹` with weak gradient `Gu` is the (weak/variational) solution of (3.23):
`∫_Ω (a ∇u·∇v + b·∇u v + c u v) dx = ∫_Ω f v dx` for all `v ∈ H₀¹`. -/
def IsVariationalSolution {d : ℕ} (Ω : Set (Euclidean d))
    (a : Euclidean d → ℝ) (b : Euclidean d → Euclidean d) (c f : Euclidean d → ℝ)
    (u : Euclidean d → ℝ) (Gu : Fin d → Euclidean d → ℝ) : Prop :=
  InH10Closure Ω u ∧ HasH1Gradient Ω u Gu ∧
  ∀ (v : Euclidean d → ℝ) (Gv : Fin d → Euclidean d → ℝ),
    InH10Closure Ω v → HasH1Gradient Ω v Gv →
      (∫ x in Ω, (a x * ∑ i, Gu i x * Gv i x + (∑ i, b x i * Gu i x) * v x
          + c x * u x * v x)) = ∫ x in Ω, f x * v x

/-- The functional (3.26):
`F(v) = ½ ∫_Ω (a|∇v|² + c v²) dx − ∫_Ω f v dx`, with `∇v = G` the weak gradient. -/
def dirichletFunctional {d : ℕ} (Ω : Set (Euclidean d))
    (a c f : Euclidean d → ℝ) (v : Euclidean d → ℝ) (G : Fin d → Euclidean d → ℝ) : ℝ :=
  (1 / 2 : ℝ) * (∫ x in Ω, (a x * ∑ i, G i x ^ 2 + c x * v x ^ 2))
    - ∫ x in Ω, f x * v x

/-- Theorem 3.7 (Dirichlet's principle). Assume (3.2) and `b = 0`. Let `f ∈ L₂` and
`u ∈ H₀¹` be the solution of (3.23). Then `F(u) ≤ F(v)` for all `v ∈ H₀¹`, with equality
only for `v = u` (equality in H₀¹, i.e. a.e. on Ω). -/
theorem dirichlet_principle {d : ℕ} (Ω : Set (Euclidean d)) (hΩ : StandingDomain Ω)
    (a : Euclidean d → ℝ) (b : Euclidean d → Euclidean d) (c : Euclidean d → ℝ)
    (hcoef : CoefficientCondition Ω a b c) (hb : b = 0)
    (f : Euclidean d → ℝ) (hf : InL2 Ω f)
    (u : Euclidean d → ℝ) (Gu : Fin d → Euclidean d → ℝ)
    (hu : IsVariationalSolution Ω a b c f u Gu) :
    ∀ (v : Euclidean d → ℝ) (Gv : Fin d → Euclidean d → ℝ),
      InH10Closure Ω v → HasH1Gradient Ω v Gv →
        dirichletFunctional Ω a c f u Gu ≤ dirichletFunctional Ω a c f v Gv ∧
        (dirichletFunctional Ω a c f u Gu = dirichletFunctional Ω a c f v Gv →
          v =ᵐ[domainMeasure Ω] u) := by
  sorry

end PDETextbook.LarssonThomee.Thm37