FE_LagP.v

(**
This file is part of the Coq Numerical Analysis library

Copyright (C) Boldo, Clément, Martin, Mayero, Mouhcine

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)

(**

 * Bibliography

 #<DIV><A NAME="ErnGuermond"></A></DIV>#
 [[ErnGuermond]]
 Alexandre Ern and Jean-Luc Guermond,
 Finite Elements I. Approximation and Interpolation,
 Texts in Applied Mathematics, vol. 72, Springer Cham, 2021,
 #<A HREF="https://doi.org/10.1007/978-3-030-56341-7">#
 https://doi.org/10.1007/978-3-030-56341-7#</A>#,#<BR>#
 #<A HREF="https://inria.hal.science/hal-03226049">#
 https://inria.hal.science/hal-03226049#</A># (with different page numbers).

 #<DIV><A NAME="RR9557v1"></A></DIV>#
 [[RR9557v1]]
 François Clément and Vincent Martin,
 Finite element method. Detailed proofs to be formalized in Coq,
 RR-9557, first version, Inria Paris, 2024,
 #<A HREF="https://inria.hal.science/hal-04713897v1">#
 https://inria.hal.science/hal-04713897v1#</A>#.
*)

From Requisite Require Import stdlib_wR ssr_wMC.

From Coquelicot Require Import Hierarchy.
Close Scope R_scope.

From Algebra Require Import Algebra_wDep.
From FEM Require Import multi_index poly_Pdk poly_LagP1k poly_LagPd1_ref.
From FEM Require Import geom_transf_affine poly_LagPd1.
From FEM Require Import LagP_node_ref LagP_node FE FE_transf.

Local Open Scope Monoid_scope.
Local Open Scope Group_scope.
Local Open Scope Ring_scope.

(* FIXME: COMMENTS TO BE REMOVED FOR PUBLICATION!
  => hypothèse k > 0 OU k.+1
  (à voir selon les énoncés pour éviter les castF et ne pas faire .+1 de façon inopinée)
  Bref : on peut mettre .+1 dans un aux mais on préfère hypothèse > 0 pour le final
*)


Section LagP_Linear_forms.

(** Lagrange linear forms. *)

(* TODO FC (07/02/2025: nnode_LagPdk ou ndof_LagPdk ?
 Ou même mettre partout (pbinom d k).+1... *)
Definition nnode_LagPdk (d k : nat) : nat := (pbinom d k).+1.
(* "nombre de dof = nombre de noeuds" *)

Lemma nnode_LagPdk_pos : forall d k, (0 < nnode_LagPdk d k)%coq_nat.
Proof. apply pbinomS_gt_0. Qed.

Lemma Pdk_has_dim : forall d k, has_dim (Pdk d k) (nnode_LagPdk d k).
Proof. intros; apply Pdk_has_dim. Qed.

Context {d : nat}.

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Def 1608, p. 102.#<BR>#
 For k=0. *)
Definition Sigma_LagPd0 vtx : '(FRd d -> R)^(nnode_LagPdk d 0) :=
  fun _ p => p (node_d0 vtx).

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Def 1608, p. 102.#<BR>#
 For k>0. *)
Definition Sigma_LagPSdSk k vtx : '(FRd d -> R)^(nnode_LagPdk d k) :=
  fun idk p => p (node vtx idk).

Lemma Sigma_lm_LagPd0 :
  forall vtx id0, lin_map (Sigma_LagPd0 vtx id0).
Proof. intros; apply pt_eval_lm. Qed.

Lemma Sigma_lm_LagPSdSk :
  forall k vtx idk, lin_map (Sigma_LagPSdSk k vtx idk).
Proof. intros; apply pt_eval_lm. Qed.

End LagP_Linear_forms.


Section Unisolvence_d0.

(** Degenerate case: unisolvence result for k=0. *)

Let PI0 d :=
  forall vtx, KerS0 (Pdk d 0) (gather (Sigma_LagPd0 vtx)).

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1615, p. 103. *)
Lemma unisolvence_inj_LagPd0 : forall {d}, PI0 d.
Proof.
unfold PI0, Sigma_LagPd0; intros;
    apply KerS0_gather_equiv; move=> p /Pd0_eq -> Hp; rewrite (Hp ord0); easy.
Qed.

End Unisolvence_d0.


Section Unisolvence_SdSk.

Let PI d k :=
  forall vtx, aff_indep_ms vtx ->
    KerS0 (Pdk d k) (gather (Sigma_LagPSdSk k vtx)).

(** Base case: unisolvence result for d=1 and k>0. *)

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1456, p. 59. *)
Lemma unisolvence_inj_LagP1Sk : forall {k}, (0 < k)%coq_nat -> PI 1 k.
Proof.
intros [| k] Hk vtx Hvtx; [easy |]; unfold PI, Sigma_LagPSdSk.
rewrite KerS0_gather_equiv; intros p Hp1 Hp2.
rewrite (LagP1k_decomp (node_inj Hvtx) Hp1); apply: lc_zero_compat_l.
apply extF; easy.
Qed.

(** Base case: unisolvence result for d>0 and k=1. *)

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1619, p. 104. *)
Lemma unisolvence_inj_LagPSd1 : forall {d}, (0 < d)%coq_nat -> PI d 1.
Proof.
intros [| d] Hd vtx Hvtx; [easy |]; unfold PI, Sigma_LagPSdSk.
rewrite KerS0_gather_equiv; intros p Hp1 Hp2.
rewrite (LagPd1_decomp_node Hvtx Hp1); apply: lc_zero_compat_l.
apply extF; easy.
Qed.

(** Induction step: unisolvence result for d,k>1. *)

Lemma unisolvence_inj_Hface :
  forall {d k}, (0 < d)%coq_nat -> (0 < k)%coq_nat -> PI d k ->
    forall i {vtx p}, aff_indep_ms vtx -> Pdk d.+1 k p ->
      (forall (idk : [0..(pbinom d.+1 k).+1)),
        Hface vtx i (node vtx idk) -> p (node vtx idk) = 0) ->
      forall x, Hface vtx i x -> p x = 0.
Proof.
intros d k Hd Hk.
assert (HSd : (0 < d.+1)%coq_nat) by apply S_pos.
destruct k as [| k1]; [easy |].
intros H i vtx p Hvtx Hp1 Hp2 x Hx.
rewrite -(f_invS_canS_r (T_geom_Hface_bijS Hvtx i) x)//.
fold (T_geom_Hface_invS Hvtx i x).
pose (pi := fun y => p (T_geom_Hface vtx i y));
  fold (pi (T_geom_Hface_invS Hvtx i x)).
rewrite (H _ vtx_ref_aff_indep pi)//; [apply T_geom_Hface_comp; easy |].
extF idk; rewrite gather_eq; unfold Sigma_LagPSdSk, pi; rewrite -node_ref_node.
destruct (ord_eq_dec i ord0) as [-> | Hi].
(* i = ord0 *)
rewrite T_geom_Hface_0_node//; apply Hp2, node_Hface_0; [easy.. |].
rewrite Adk_inv_correct_r inj_CSdk_sum//.
(* i <> ord0 *)
rewrite T_geom_Hface_S_node//; apply Hp2, (node_Hface_S Hvtx _ Hi HSd Hk).
rewrite Adk_inv_correct_r;
    [apply insertF_correct_l; easy | rewrite inj_ASdki_sum; apply Adk_sum].
Qed.

Lemma unisolvence_inj_ind_LagPSdSk :
  forall {d k}, (0 < d)%coq_nat -> (0 < k)%coq_nat ->
    PI d.+1 k -> PI d k.+1 -> PI d.+1 k.+1.
Proof.
unfold PI, Sigma_LagPSdSk; intros d k Hd0 Hk0 IH2 IH1 vtx Hvtx.
assert (HSd0 : (0 < d.+1)%coq_nat) by apply S_pos.
assert (HSk0 : (0 < k.+1)%coq_nat) by apply S_pos.
assert (HSk1 : (1 < k.+1)%coq_nat) by now rewrite -Nat.succ_lt_mono.
apply KerS0_gather_equiv; intros p Hp1 Hp2.
(* Step 1: factorization. *)
destruct (factor_zero_on_Hface HSd0 Hvtx ord0 Hp1) as [q [Hq1 Hq2]].
apply (unisolvence_inj_Hface Hd0 HSk0); easy.
(* Step 2: cancellation. *)
specialize (IH2 _ (sub_vtx_aff_indep HSk1 Hvtx));
    rewrite KerS0_gather_equiv in IH2.
rewrite Hq2 (IH2 q)//; [apply: mult_zero_r |]; intros iSdk.
apply mult_reg_l with (LagPd1 Hvtx ord0 (sub_node k.+1 vtx iSdk)).
apply invertible_equiv_R; rewrite -Hface_equiv;
    apply sub_node_out_Hface_0; easy.
rewrite -fct_mult_eq -Hq2 sub_node_eq// Hp2 mult_zero_r; easy.
Qed.

(** Unisolvence result for d,k>0. *)

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1625, pp. 105-107. *)
Theorem unisolvence_inj_LagPSdSk :
  forall {d k}, (0 < d)%coq_nat -> (0 < k)%coq_nat -> PI d k.
Proof.
apply nat_ind2_alt_11.
intros d; apply: unisolvence_inj_LagPSd1.
intros k; apply: unisolvence_inj_LagP1Sk.
intros d k; apply: unisolvence_inj_ind_LagPSdSk.
Qed.

End Unisolvence_SdSk.


Section FE_LagPSdSk.

(** Build Lagrange FE for d,k>0. *)

(** We need a dimension and an approximation degree which are structurally nonzero.
 We take d = d1.+1 and k = k1.+1 with d1 k1 : nat. *)
Context {d1 : nat}.
Variable k1 : nat.
Let d := d1.+1.
Let k := k1.+1.
Context {vtx : 'R^{d.+1,d}}.
Hypothesis Hvtx : aff_indep_ms vtx.

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Th 1629, p. 108.#<BR>#
 For k>0. *)
Definition FE_LagPSdSk : FE d :=
  mk_FE Simplex vtx (nnode_LagPdk d k) (Pdk_has_dim d k)
    (Sigma_lm_LagPSdSk k vtx) (unisolvence_inj_LagPSdSk S_pos S_pos _ Hvtx).

End FE_LagPSdSk.


Section FE_LagPSdk.

(** Build Lagrange FE for d>0. *)

(** We need a dimension which is structurally nonzero.
 We take d = d1.+1 with d1 : nat. *)
Context {d1 : nat}.
Let d := d1.+1.

(** For k=0. *)

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Th 1629, p. 108.#<BR>#
 For k=0. *)
Definition FE_LagPd0 (vtx : 'R^{d.+1,d}) : FE d :=
  mk_FE Simplex vtx (nnode_LagPdk d 0) (Pdk_has_dim d 0)
    (Sigma_lm_LagPd0 vtx) (unisolvence_inj_LagPd0 vtx).

(** For k>=0. *)

Variable k : nat.
Context {vtx : 'R^{d.+1,d}}.
Hypothesis Hvtx : aff_indep_ms vtx.

Definition FE_LagPSdk : FE d :=
  match k with
  | 0 => FE_LagPd0 vtx
  | S k1 => FE_LagPSdSk k1 Hvtx
  end.

End FE_LagPSdk.


Section FE_LagPdk.

(** Build Lagrange FE for any d,k. *)

(** For d=0 and any k. *)

Definition FE_LagP0k : FE 0 := FE_d_0 Simplex one_not_zero_R.

(** For d,k>=0. *)

Context {d : nat}.
Variable k : nat.
Context {vtx : 'R^{d.+1,d}}.
Hypothesis Hvtx : aff_indep_ms vtx.

Let Hvtx' :
  forall {d'} (Hd : d = d'),
    aff_indep_ms (castF (eq_S _ _ Hd) (mapF (castF Hd) vtx)).
Proof. intros; subst; rewrite !castF_id; easy. Qed.

Definition FE_LagPdk : FE d :=
  match d as d' return d = d' -> _ with
  | 0 => fun Hd => FE_LagP0k
  | S _ => fun Hd => FE_LagPSdk k (Hvtx' Hd)
  end erefl.

End FE_LagPdk.


Section FE_LagPdk_from_ref.

(** Check that
 [FE_LagPdk k Hvtx = FE_ref_to_cur (FE_LagPdk k vtx_ref_aff_indep)],
 ie that the generic simplicial Lagrange FE is equal to the transformation of
 the reference simplicial Lagrange FE. *)

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1604 (for k=0), p. 101.#<BR>#
 It is actually directly on the linear forms, rather than on the nodes. *)
Lemma Sigma_LagPd0_eq_from_ref :
  forall {d} (vtx : 'R^{d.+1,d}) id0 p,
    Sigma_LagPd0 vtx id0 p =
      Sigma_LagPd0 vtx_ref id0 (p \o (T_geom vtx)).
Proof.
intros; unfold Sigma_LagPd0, node_d0; rewrite comp_correct T_geom_am.
do 2 apply f_equal; extF; apply eq_sym, T_geom_vtx.
rewrite sum_ones_R; apply invertible_equiv_R, R_compl.INR_n0, S_pos.
Qed.

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1604 (for k>0), p. 101.#<BR>#
 It is actually directly on the linear forms, rather than on the nodes. *)
Lemma Sigma_LagPSdSk_eq_from_ref :
  forall {d} k (vtx : 'R^{d.+1,d}) i p,
    Sigma_LagPSdSk k vtx i p =
      Sigma_LagPSdSk k vtx_ref i (p \o (T_geom vtx)).
Proof.
intros; unfold Sigma_LagPSdSk.
rewrite comp_correct -node_ref_node T_geom_node; easy.
Qed.

Variable d k : nat.
Context {vtx : 'R^{d.+1,d}}.
Hypothesis Hvtx : aff_indep_ms vtx.

Definition FE_LagPdk_ref := FE_LagPdk k (@vtx_ref_aff_indep d).

Lemma FE_LagPdk_eq_from_ref :
  FE_LagPdk k Hvtx = FE_from_ref FE_LagPdk_ref vtx Hvtx.
Proof.
unfold FE_LagPdk_ref, FE_LagPdk; destruct d as [| d1]; [| destruct k as [| k1]].
(* d = 0 *)
apply FE_d_0_uniq with 1; [easy | apply one_not_zero_R |..];
    simpl; exists ord0; rewrite singleF_0 fct_scal_eq scal_one_l//.
(* d <> 0, k = 0 *)
apply: FE_ext; simpl; [| unfold P_transf | unfold Sigma_transf]; simpl.
rewrite !castF_id; extF; rewrite !mapF_correct T_geom_vtx; easy.
apply eq_sym, transf_am_Pdk; [apply T_geom_bij; easy | apply T_geom_am].
rewrite !castF_id; fun_ext2; apply Sigma_LagPd0_eq_from_ref.
(* d,k <> 0 *)
apply: FE_ext; simpl; [| unfold P_transf | unfold Sigma_transf]; simpl.
extF; do 3 rewrite {1}castF_id; rewrite !mapF_correct; do 2 rewrite castF_id.
rewrite T_geom_vtx; easy.
apply eq_sym, transf_am_Pdk; [apply T_geom_bij; easy | apply T_geom_am].
rewrite !castF_id; fun_ext2; apply Sigma_LagPSdSk_eq_from_ref.
Qed.

Let Hndof := f_equal ndof FE_LagPdk_eq_from_ref.

Lemma shape_fun_LagPdk_eq_from_ref :
  shape_fun (FE_LagPdk k Hvtx) =
    castF (eq_sym Hndof) (shape_fun (FE_from_ref FE_LagPdk_ref vtx Hvtx)).
Proof. apply shape_fun_ext. Qed.

(* Useless... *)
Lemma local_interp_LagPdk_comp :
  forall f,
    local_interp FE_LagPdk_ref (f \o (T_geom vtx)) =
      (local_interp (FE_LagPdk k Hvtx) f) \o (T_geom vtx).
Proof.
intros; rewrite (local_interp_transf (T_geom_bij Hvtx))
    FE_LagPdk_eq_from_ref; easy.
Qed.

End FE_LagPdk_from_ref.


Section Face_unisolvence.

(** Face unisolvence for d>1 and k>0. *)

Context {d k : nat}.
Hypothesis Hd : (1 < d)%coq_nat.
Hypothesis Hk : (0 < k)%coq_nat.
Context {vtx : 'R^{d.+1,d}}.
Hypothesis Hvtx : aff_indep_ms vtx.

(**
 #<A HREF="##RR9557v1">#[[RR9557v1]]#</A>#
 Lem 1628, p. 107. *)
Lemma Hface_unisolvence :
  forall {i} p, Pdk d k p ->
    (forall (idk : 'I_(pbinom d k).+1),
      Hface vtx i (node vtx idk) -> p (node vtx idk) = 0) <->
    (forall x, Hface vtx i x -> p x = 0).
Proof.
intros i p Hp; split; [| move=> H idk Hidk; apply H; easy].
destruct d as [| d1]; [easy |].
assert (Hd0 : (0 < d1)%coq_nat) by now apply Nat.succ_lt_mono.
apply: (unisolvence_inj_Hface _ _ (unisolvence_inj_LagPSdSk _ _)); easy.
Qed.

End Face_unisolvence.