Library FEM.FE_LagP

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

[ErnGuermond] Alexandre Ern and Jean-Luc Guermond, Finite Elements I. Approximation and Interpolation, Texts in Applied Mathematics, vol. 72, Springer Cham, 2021, https://doi.org/10.1007/978-3-030-56341-7,
https://inria.hal.science/hal-03226049 (with different page numbers).
[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, https://inria.hal.science/hal-04713897v1.

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 geom_simplex geom_transf_affine.
From FEM Require Import poly_Pdk poly_LagPd1_ref poly_LagPd1_cur poly_LagP1k.
From FEM Require Import FE FE_transf.

Local Open Scope R_scope.

Section Unisolvence_dSk.

FE LAGRANGE current

Definition nnode_LagPdk dL kL : nat := (pbinom dL kL).+1.

Lemma nnode_LagPdk_pos : dL kL, (0 < nnode_LagPdk dL kL)%coq_nat.
Proof.
intros; apply pbinomS_gt_0.
Qed.

Lemma Pdk_has_dim_ndof_LagPdk : dL kL,
  has_dim (Pdk dL kL) (nnode_LagPdk dL kL).
Proof.
intros; apply Pdk_has_dim.
Qed.


Definition Sigma_LagPdSk_cur {dL} kL vtx_cur : '(FRd dL R)^(nnode_LagPdk dL kL) :=
  fun i pp (node_cur vtx_cur i).

[RR9557v1] Lem 1619, p. 104.
Lemma unisolvence_inj_LagPd1_cur : {dL} {vtx_cur:'R^{dL.+1,dL}},
  (0 < dL)%coq_nat
  affine_independent vtx_cur
  KerS0 (Pdk dL 1) (gather (Sigma_LagPdSk_cur 1 vtx_cur)).
Proof.
destruct dL; try easy.
intros vtx_cur dL_pos Hvtx.
rewrite KerS0_gather_equiv.
intros p Hp1 Hp2.
rewrite (LagPd1_cur_decomp_Pd1 Hvtx p Hp1).
apply: lc_zero_compat_l.
apply extF; intros i; rewrite fct_zero_eq.
rewrite (vtx_cur_node_cur_d1 vtx_cur).
apply Hp2.
Qed.

[RR9557v1] Lem 1456, p. 59.
Lemma unisolvence_inj_LagP1k_cur : {kL vtx_cur},
  (0 < kL)%coq_nat
  affine_independent vtx_cur
   KerS0 (Pdk 1 kL) (gather (Sigma_LagPdSk_cur kL vtx_cur)).
Proof.
intros kL vtx_cur HkL Hvtx.
rewrite KerS0_gather_equiv.
destruct kL; try easy.
intros p Hp1 Hp2.
destruct (basis_ex_decomp (LagP1k_cur_basis Hvtx) p Hp1) as [L HL].
rewrite HL in Hp2.
rewrite HL.
apply: lc_zero_compat_l; apply extF; intros i.
specialize (Hp2 i).
rewrite fct_zero_eq.
rewrite -Hp2.
unfold Sigma_LagPdSk_cur.
rewrite fct_lc_r_eq.
rewrite -{1}(lc_kron_r L).
rewrite fct_lc_r_eq.
apply lc_ext_r; intros j.
rewrite LagP1k_cur_kron_node_alt; try easy.
apply kron_sym.
Qed.

[RR9557v1] Lem 1621, pp. 104-105.
Lemma factorize_poly_zero_on_face_ref_d :
   {dL kL : nat}, (0 < dL)%coq_nat
   p_ref : FRd dL, Pdk dL kL.+1 p_ref
  ( x_ref : 'R^dL, face_ref ord_max x_ref p_ref x_ref = 0)
     q_Ref : FRd dL, Pdk dL kL q_Ref
         p_ref = ((LagPd1_ref ord_max) × q_Ref)%K.
Proof.
intros dL kL Hd p_ref Hp1; split.
destruct dL; try easy.
intros Hp2.
destruct (PdSk_split kL p_ref Hp1) as [p0 [p1 [H1 [H2 H3]]]].
p1.
split; try easy.
rewrite H3.
apply fun_ext; intros x.
rewrite fct_mult_eq.
assert (H : y : 'R^dL, p0 y = 0).
intros y; rewrite -(Rplus_0_r (p0 _)).
rewrite -{1}(Rmult_0_l (p1 (insertF y 0 ord_max))).
replace (_ + _) with (p_ref (insertF y 0 ord_max)).
2 : try rewrite H3 widenF_S_insertF_max insertF_correct_l; easy.
apply Hp2.
rewrite face_ref_equiv.
rewrite (LagPd1_ref_S ord_max_not_0).
rewrite insertF_correct_l; try easy.
apply ord_inj; rewrite lower_S_correct; easy.
rewrite H Rplus_0_l.
rewrite (LagPd1_ref_S ord_max_not_0).
unfold mult; simpl.
do 2!f_equal.
apply ord_inj; rewrite lower_S_correct; easy.
intros [q [Hq1 Hq2]] x Hx.
rewrite Hq2 fct_mult_eq.
rewrite face_ref_equiv in Hx.
rewrite Hx mult_zero_l; easy.
Qed.

[RR9557v1] Lem 1623, p. 105
(for i=0 only, and with the transposition instead of the circular permutation, see Rem 1622, p. 105).
Lemma factorize_poly_zero_on_face0 :
   {dL kL : nat} {vtx_cur : 'R^{dL.+1,dL}}
   (Hvtx: affine_independent vtx_cur),
  (0 < dL)%coq_nat
   p : FRd dL, Pdk dL kL.+1 p
  ( x : 'R^dL, face Hvtx ord0 x p x = 0)
     q : FRd dL, Pdk dL kL q
        p = ((LagPd1_cur Hvtx ord0) × q)%K.
Proof.
intros dL kL vtx_cur Hvtx Hd p Hp1; split.
intros Hp2.
pose (p_ref := p \o T_geom_d_0 vtx_cur).
assert (H : Pdk dL kL.+1 p_ref).
apply T_geom_d_0_comp; easy.
destruct (proj1 (factorize_poly_zero_on_face_ref_d Hd p_ref H)) as [q_ref [Y1 Y2]].
intros x_ref Hx_ref.
unfold p_ref, comp.
apply Hp2.
apply (T_geom_d_0_in_face0 Hvtx x_ref); easy.
(q_ref \o T_geom_d_0_inv Hvtx).
split.
apply T_geom_d_0_inv_comp; try easy.
assert (H0 : p = p_ref \o T_geom_d_0_inv Hvtx).
unfold p_ref, comp.
apply fun_ext; intros x_cur.
rewrite f_inv_can_r; easy.
rewrite H0 Y2.
unfold comp.
apply fun_ext; intros x.
rewrite 2!fct_mult_eq.
f_equal.
rewrite LagPd1_cur_T_geom_d_0_inv.
rewrite transpF_correct_l0; easy.
intros [q [Hq1 Hq2]] x Hx.
rewrite Hq2 fct_mult_eq Hx mult_zero_l; easy.
Qed.

[RR9557v1] Lem 1625, pp. 105-107.
Theorem unisolvence_inj_LagPdSk_cur : {dL kL}
  (HdL : (0 < dL)%coq_nat) (HkL : (0 < kL)%coq_nat) {vtx_cur :'R^{dL.+1,dL}},
  affine_independent vtx_cur
   KerS0 (Pdk dL kL) (gather (Sigma_LagPdSk_cur kL vtx_cur)).
Proof.
apply: nat_ind2_alt_11.
intros kL HkL vtx_cur Hvtx.
apply unisolvence_inj_LagPd1_cur; easy.
intros dL HdL vtx_cur Hvtx.
apply unisolvence_inj_LagP1k_cur; easy.
intros dL kL HdL HkL IH2 IH1 vtx_cur Hvtx.
apply KerS0_gather_equiv.
intros p Hp1 Hp2.
specialize (IH1 vtx_ref vtx_ref_affine_independent).
rewrite KerS0_gather_equiv in IH1.
pose (p0 := p \o T_geom_face vtx_cur ord0).
assert (Hp0 : Pdk dL kL.+1 p0) by now apply T_geom_face_comp.
specialize (IH1 p0 Hp0).
unfold Sigma_LagPdSk_cur in IH1.
rewrite -node_ref_cur in IH1; try now apply (Nat.lt_0_succ kL).
assert (Y1 : x : 'R^dL.+1,
  face Hvtx ord0 x p x = 0).
intros x Hx.
assert (Y2 : y : 'R^dL.+1, face Hvtx ord0 y
  p y = p0 (T_geom_face_inv Hvtx ord0 y)).
unfold p0.
intros y Hy.
unfold comp.
rewrite (f_invS_canS_r (T_geom_face_bijS Hvtx ord0)); try easy.
rewrite Y2; try easy.
rewrite IH1; try easy.
intros i; unfold p0, comp.
rewrite T_geom_face0_map_node; try easy.
apply Hp2.
apply Nat.lt_0_succ.
destruct (proj1 (factorize_poly_zero_on_face0 Hvtx (Nat.lt_0_succ dL) p Hp1) Y1) as [q [Hq1 Hq2]].
rewrite Hq2.
assert (H1 : (1 < kL.+1)%coq_nat); auto with arith.
specialize (IH2 (sub_vtx (kL.+1) vtx_cur) (sub_vtx_affine_independent H1 Hvtx)).
rewrite KerS0_gather_equiv in IH2.
unfold Sigma_LagPdSk_cur in IH2.
rewrite (IH2 q); try easy.
apply fun_ext; intros x; simpl.
now rewrite mult_zero_r.
intros i.
apply Rmult_eq_reg_l with (LagPd1_cur Hvtx ord0
  (node_cur (sub_vtx kL.+1 vtx_cur) i)).
2: generalize (sub_node_out_face0 Hvtx i); intros H4.
2: apply H4; try easy; auto with arith.
rewrite Rmult_0_r.
apply trans_eq with (p (node_cur (sub_vtx kL.+1 vtx_cur) i)).
rewrite Hq2; easy.
apply trans_eq with (p (sub_node kL.+1 vtx_cur i)); try easy.
rewrite sub_node_cur_eq; try easy.
apply Hp2.
Qed.

End Unisolvence_dSk.

Section Face_unisolvence_dSk.

Context {dL kL : nat}.
Hypothesis dL_pos : (1 < dL)%coq_nat.
Hypothesis kL_pos : (0 < kL)%coq_nat.

Context {vtx_cur :'R^{dL.+1,dL}}.
Hypothesis Hvtx : affine_independent vtx_cur.

[RR9557v1] Lem 1628, p. 107.
Lemma face0_unisolvence : p,
  Pdk dL kL p
  ( ipk:'I_(pbinom dL kL).+1, ((pbinom dL kL.-1).+1 ipk)%coq_nat
      p (node_cur vtx_cur ipk) = 0)
   
   ( x, face Hvtx ord0 x p x = 0).
Proof.
intros p Hp; split; intros H1.
intros x Hx.
case_eq dL.
intros H; contradict dL_pos; subst; easy.
intros pdL HpdL; subst.
replace x with (T_geom_face vtx_cur ord0 (T_geom_face_inv Hvtx ord0 x)).
2: rewrite f_invS_canS_r; easy.
pose (p0:= fun yp (T_geom_face vtx_cur ord0 y));
  fold (p0 (T_geom_face_inv Hvtx ord0 x)).
assert (HpdL : (0 < pdL)%coq_nat) by auto with arith.
replace p0 with (@fct_zero 'R^pdL R_AbelianMonoid); try easy.
apply sym_eq, (unisolvence_inj_LagPdSk_cur HpdL kL_pos
   vtx_ref_affine_independent).
apply Pdk_comp_am; try easy.
apply T_geom_face_am; easy.
apply extF; intros i.
unfold gather, swap; simpl.
unfold Sigma_LagPdSk_cur.
unfold fct_zero, p0.
replace (node_cur vtx_ref i) with (node_ref i).
rewrite T_geom_face0_map_node; try easy.
apply H1; try easy.
apply Adk_last_layer; try auto with arith.
rewrite Adk_inv_correct_r.
apply T_node_face_0_sum.
rewrite T_node_face_0_sum; easy.
rewrite node_ref_cur; easy.
intros ipk Hipk.
apply H1.
apply node_face0_in_Cdk; auto with arith.
Qed.

Lemma faceS_unisolvence : p (i:'I_dL),
  Pdk dL kL p
  ( ipk, Adk dL kL ipk i = O
      p (node_cur vtx_cur ipk) = 0)
   
   ( x, face Hvtx (lift_S i) x p x = 0).
Proof.
intros p i Hp; split.
2: intros H ipk Hipk.
2: rewrite (node_faceSj_in_Adkj Hvtx) in Hipk; auto with zarith.
destruct dL.
contradict dL_pos; subst; easy.
clear dL; generalize n dL_pos p i Hp vtx_cur Hvtx; clear n dL_pos p i Hp vtx_cur Hvtx.
intros pdL HpdL p i Hp vtx_cur Hvtx.
intros H x Hx.
replace x with (T_geom_face vtx_cur (lift_S i) (T_geom_face_inv Hvtx (lift_S i) x)).
2: rewrite f_invS_canS_r; easy.
pose (p0:= fun yp (T_geom_face vtx_cur (lift_S i) y));
  fold (p0 (T_geom_face_inv Hvtx (lift_S i) x)).
assert (HpdL2 : (0 < pdL)%coq_nat) by auto with arith.
replace p0 with (@fct_zero 'R^pdL R_AbelianMonoid); try easy.
apply sym_eq, (unisolvence_inj_LagPdSk_cur HpdL2 kL_pos
   vtx_ref_affine_independent).
apply Pdk_comp_am; try easy.
apply T_geom_face_am; easy.
apply extF; intros j.
unfold gather, swap; simpl.
unfold Sigma_LagPdSk_cur.
unfold fct_zero, p0.
replace (node_cur vtx_ref j) with (node_ref j).
rewrite T_geom_faceS_map_node; try easy.
apply H; try easy.
rewrite Adk_inv_correct_r.
unfold T_node_face_S; rewrite insertF_correct_l; easy.
rewrite T_node_face_S_sum; apply Adk_sum.
rewrite node_ref_cur; easy.
Qed.

End Face_unisolvence_dSk.

Section Unisolvence_d0.


Context {dL : nat}.

Variable vtx_cur :'R^{dL.+1,dL}.

[RR9557v1] Def 1588, Eq (9.80), p. 97.
Let node_iso : 'R^dL := isobarycenter (vtx_cur : 'fct_ModuleSpace^dL.+1).

Definition Sigma_0 : '(FRd dL R)^(nnode_LagPdk dL 0) :=
   fun _ pp node_iso.

[RR9557v1] Lem 1615, p. 103.
Lemma unisolvence_inj_LagPd0 :
  KerS0 (Pdk dL 0) (gather (Sigma_0)).
Proof.
apply KerS0_gather_equiv.
intros p Hp1 Hp2.
specialize (Hp2 ord0).
rewrite Pd0_eq in Hp1.
rewrite Hp1 in Hp2.
unfold Sigma_0 in Hp2; rewrite Hp2 in Hp1; easy.
Qed.

End Unisolvence_d0.

Section FE_LagPdk_cur.

FE LAGRANGE current

Context {dL : nat}.
Hypothesis dL_pos : (0 < dL)%coq_nat.
Variable kL : nat.

Variable vtx_cur :'R^{dL.+1,dL}.

[RR9557v1] Def 1608, p. 102.
Definition Sigma_LagPdk_cur : '(FRd dL R)^(nnode_LagPdk dL kL)
 := match kL with
   | OSigma_0 vtx_cur
   | S nSigma_LagPdSk_cur (S n) vtx_cur
end.

Lemma Sigma_LagPdk_cur_S : (0 < kL)%coq_nat
  Sigma_LagPdk_cur = Sigma_LagPdSk_cur kL vtx_cur.
Proof.
intros H; unfold Sigma_LagPdk_cur; case_eq kL; try easy.
intros H1; subst; easy.
Qed.

Lemma Sigma_LagPdk_cur_lm :
   i,
   lin_map (Sigma_LagPdk_cur i).
Proof.
intros i; unfold Sigma_LagPdk_cur; destruct kL; apply lm_pt_eval.
Qed.

Hypothesis Hvtx : affine_independent vtx_cur.

[RR9557v1] Th 1629, p. 108.
[RR9557v1] Lem 1604, p. 101.
Lemma Sigma_LagPdk_cur_ref : (vtx_cur :'R^{dL.+1,dL}) i p,
 Sigma_LagPdk_cur kL vtx_cur i p =
 Sigma_LagPdk_cur kL vtx_ref i
  (transf (T_geom vtx_cur) p).
Proof.
intros vtx_cur; unfold Sigma_LagPdk_cur; case kL.
intros i p; unfold Sigma_0; unfold transf.
f_equal; unfold isobarycenter.
rewrite T_geom_am.
f_equal; unfold mapF, mapiF; apply extF; intros j.
now rewrite T_geom_transports_vtx.
apply invertible_equiv_R; unfold plusn1, plusn.
rewrite sum_constF sum_ones_R; apply: scal_not_zero_R.
apply not_0_INR; discriminate.
unfold one, zero; simpl; apply R1_neq_R0.
intros k i; unfold Sigma_LagPdSk_cur, transf.
rewrite -node_ref_cur.
rewrite T_geom_transports_node; easy.
Qed.

Lemma FE_ref_to_cur_LagPdk_ref_eq : {vtx_cur :'R^{dL.+1,dL}}
  (Hvtx : affine_independent vtx_cur),
  FE_LagPdk_cur dL_pos kL vtx_cur Hvtx =
    FE_ref_to_cur FE_LagPdk_ref vtx_cur Hvtx.
Proof.
intros vtx_cur Hvtx.
pose (fe1 := FE_LagPdk_cur dL_pos kL vtx_cur Hvtx); fold fe1.
pose (fe2 := FE_ref_to_cur FE_LagPdk_ref vtx_cur Hvtx); fold fe2.
assert (Hd : d fe2 = d fe1) by easy.
assert (Hn : ndof fe2 = ndof fe1) by easy.
assert (Hs : shape fe2 = shape fe1) by easy.
apply (FE_ext Hd Hn Hs); simpl.
rewrite castT_id.
apply extF; intros i; rewrite mapF_correct.
rewrite T_geom_transports_vtx; easy.
rewrite cast2F_fun_id; apply subset_ext_equiv; split; intros p Hp.
simpl; apply Pdk_comp_am; try easy.
apply T_geom_am.
simpl in Hp.
apply Pdk_comp_am_rev with (T_geom vtx_cur); try easy.
apply T_geom_am.
apply T_geom_bij; easy.
rewrite castF_id cast2F_fun_id; apply extF; intros i.
rewrite mapF_correct; apply fun_ext; intros p; simpl.
apply Sigma_LagPdk_cur_ref.
Qed.

Lemma shape_fun_LagPdk_cur_eq : {vtx_cur :'R^{dL.+1,dL}}
  (Hvtx : affine_independent vtx_cur),
  shape_fun (FE_LagPdk_cur dL_pos kL vtx_cur Hvtx) =
    shape_fun (FE_ref_to_cur FE_LagPdk_ref vtx_cur Hvtx).
Proof.
intros vtx_cur Hvtx.
rewrite (shape_fun_ext (FE_ref_to_cur_LagPdk_ref_eq Hvtx)) castF_id castF_fun_id; easy.
Qed.

[ErnGuermond] Prop 9.3, p. 103.
Ip (v o T) = (Ip v) o T.
Lemma local_interp_LagPdk_comp :
   {vtx_cur :'R^{dL.+1,dL}}
  (Hvtx : affine_independent vtx_cur) (v : FRd dL),
  local_interp FE_LagPdk_ref (fun x_ref : 'R^dLv (T_geom vtx_cur x_ref)) =
   (fun x_ref : 'R^dLlocal_interp (FE_LagPdk_cur dL_pos kL vtx_cur Hvtx)
         v (T_geom vtx_cur x_ref)).
Proof.
intros v vtx_cur Hvtx.
unfold local_interp.
apply fun_ext; intros x.
rewrite 2!fct_lc_r_eq.
apply lc_scalF_compat.
rewrite shape_fun_LagPdk_cur_eq.
simpl; apply extF; intros i; unfold scalF, map2F.
rewrite shape_fun_transf.
unfold transf_inv; rewrite f_inv_can_l.
f_equal; now rewrite -Sigma_LagPdk_cur_ref.
Qed.

End FE_LagPdk_ref.