Library FEM.FE

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

[Ern] Alexandre Ern, Éléments finis, Aide-mémoire, Dunod/L'Usine Nouvelle, 2005, https://www.dunod.com/sciences-techniques/aide-memoire-elements-finis.
[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 geometry.

Local Open Scope R_scope.

Section FE_Facts1.

Inductive shape_type := Simplex | Quad.

Definition nvtx_of_shape (d : nat) (shape : shape_type) : nat :=
  match shape with
    | Simplexd.+1
    | QuadNat.pow 2 d
  end.

Lemma nos_eq :
   {d1 d2 s1 s2},
    d1 = d2 s1 = s2 nvtx_of_shape d1 s1 = nvtx_of_shape d2 s2.
Proof. move=>>; apply f_equal2. Qed.


[Ern] Def 4.1, p. 75.
[ErnGuermond] Def 5.2, p. 51.
[RR9557v1] Def 1424, p. 53.
Record FE : Type := mk_FE {
  d : nat ;
  ndof : nat ;
  d_pos : (0 < d)%coq_nat ;
  ndof_pos : (0 < ndof)%coq_nat ;
  shape : shape_type ;
  nvtx : nat := nvtx_of_shape d shape ;
  vtx : '('R^d)^nvtx ;
  K_geom : 'R^d Prop := convex_envelop vtx ;
  P_approx : FRd d Prop ;
  P_approx_has_dim : has_dim P_approx ndof ;
  Sigma : '(FRd d R)^ndof ;
  Sigma_lm : i, lin_map (Sigma i) ;
  unisolvence_inj : KerS0 P_approx (gather Sigma) ;
}.

Lemma unisolvence : (f:FE),
     bijS (P_approx f) fullset (gather (Sigma f)).
Proof.
intros f.
apply lmS_bijS_KerS0 with (ndof f) (ndof f); try easy.
apply (P_approx_has_dim f).
apply has_dim_Rn.
apply gather_lm_compat, (Sigma_lm f).
apply (unisolvence_inj f).
Qed.

Lemma FE_ext :
   {fe1 fe2 : FE} (Hd : d fe2 = d fe1)
      (Hn : ndof fe2 = ndof fe1) (Hs : shape fe2 = shape fe1),
    vtx fe1 = castT (nos_eq Hd Hs) Hd (vtx fe2)
    P_approx fe1 = cast2F_fun Hd (P_approx fe2)
    Sigma fe1 = castF Hn (mapF (cast2F_fun Hd) (Sigma fe2))
    fe1 = fe2.
Proof.
intros fe1 fe2 Hd Hn Hs; destruct fe1, fe2; simpl in *; subst.
rewrite castT_id !cast2F_fun_id castF_id; intros; subst.
f_equal; apply proof_irrel.
Qed.

Variable fe: FE.

Lemma nvtx_pos : (0 < nvtx fe)%coq_nat.
Proof.
unfold nvtx, nvtx_of_shape; case shape.
auto with arith.
apply pow_pos; auto.
Qed.

Lemma shape_dec :
  {shape fe = Simplex} + {shape fe = Quad}.
Proof.
case (shape fe); [left | right]; easy.
Qed.

Lemma nvtx_dec :
  {nvtx fe = S (d fe)} + {nvtx fe = (Nat.pow 2 (d fe))%nat }.
Proof.
unfold nvtx; case (shape fe); [left|right]; easy.
Qed.

Lemma nvtx_Simplex :
  shape fe = Simplex nvtx fe = S (d fe).
Proof.
unfold nvtx; case (shape fe); auto.
intros H; discriminate.
Qed.

Lemma nvtx_Quad :
  shape fe = Quad nvtx fe = (Nat.pow 2 (d fe))%nat.
Proof.
unfold nvtx; case (shape fe); auto.
intros H; discriminate.
Qed.

Definition is_local_shape_fun :=
  fun (theta : '(FRd (d fe))^(ndof fe)) ⇒
  ( i : 'I_(ndof fe), P_approx fe (theta i))
  ( i j : 'I_(ndof fe), Sigma fe i (theta j) = kronR i j).

Lemma is_local_shape_fun_inj : theta1 theta2,
  is_local_shape_fun theta1 is_local_shape_fun theta2
    theta1 = theta2.
Proof.
intros theta1 theta2 [H1 H1'] [H2 H2'].
pose (PC:= (has_dim_cms _ (P_approx_has_dim fe))).
generalize (lmS_bijS_val_full_equiv _ (P_approx_has_dim fe) (ndof fe) has_dim_Rn
    (proj1 (gather_lm_compat _) (Sigma_lm fe)) eq_refl); intros T.
destruct T as [T1 _].
fold PC in T1.
specialize (T1 (unisolvence fe)).
apply fun_ext; intros j.
apply minus_zero_reg_sym.
generalize (cms_minus PC _ _ (H2 j) (H1 j)); intros H.
apply trans_eq with (val (mk_sub H)); try easy.
eapply trans_eq with (val zero); try easy.
apply val_eq, T1.
simpl; unfold gather, swap; simpl.
apply fun_ext; intros i.
unfold minus; rewrite (proj1 (Sigma_lm fe i)).
rewrite lm_opp.
2: apply Sigma_lm.
rewrite H1' H2'.
rewrite plus_opp_r; easy.
Qed.

[Ern] p. 76.
[ErnGuermond] Prop 5.5, pp. 51-52.
Definition shape_fun : '(FRd (d fe))^(ndof fe) := predual_basis (unisolvence fe).

Lemma shape_fun_correct : is_local_shape_fun shape_fun.
Proof.
split.
apply predual_basis_in_sub.
intros; apply predual_basis_dualF.
Qed.


Lemma shape_fun_basis : basis (P_approx fe) shape_fun.
Proof.
apply predual_basis_basis.
apply fe.
apply (Sigma_lm fe).
Qed.

[Ern] Def 4.4, p. 78.
Definition local_interp : FRd (d fe) FRd (d fe) :=
  fun vlin_comb (fun iSigma fe i v) shape_fun.

Lemma local_interp_is_poly : v : FRd (d fe),
      (P_approx fe) (local_interp v).
Proof.
intros v; unfold local_interp.
rewrite (proj1 shape_fun_basis); easy.
Qed.

Lemma local_interp_lm :
    lin_map local_interp.
Proof.
split.
intros v1 v2.
unfold local_interp.
rewrite -lc_plus_l.
apply lc_ext_l; intros i.
apply (Sigma_lm fe).
intros v l.
unfold local_interp.
rewrite -lc_scal_l.
apply lc_ext_l; intros i.
apply (Sigma_lm fe).
Qed.

Lemma local_interp_shape_fun : i : 'I_(ndof fe),
  local_interp (shape_fun i) = shape_fun i.
Proof.
intros j.
unfold local_interp.
apply trans_eq with
   (lin_comb (fun i : 'I_(ndof fe)kronR i j) shape_fun).
apply lc_scalF_compat, extF; intros i.
f_equal; apply extF; intro.
apply shape_fun_correct.
apply lc_kron_l_in_l.
Qed.

[Ern] Eq (4.6), p. 78.
Lemma local_interp_proj : p : FRd (d fe),
  (P_approx fe) p local_interp p = p.
Proof.
intros p Hp.
rewrite (proj1 shape_fun_basis) in Hp.
inversion Hp as [ L HL ].
rewrite f_lc_compat_lm.
apply lc_ext_r; intros i.
apply local_interp_shape_fun.
apply local_interp_lm.
Qed.

[Ern] p. 78.
Lemma local_interp_idem : v : FRd (d fe),
  local_interp (local_interp v) = local_interp v.
Proof.
intros v; rewrite local_interp_proj; try easy.
apply local_interp_is_poly; try easy.
Qed.

[Ern] Eq (4.4), p. 77.
Lemma Sigma_local_interp : (i : 'I_(ndof fe))
  (v : FRd (d fe)),
    Sigma fe i (local_interp v) = Sigma fe i v.
Proof.
intros i v.
unfold local_interp.
rewrite f_lc_compat_lm.
apply trans_eq with
  (lin_comb (fun jkronR i j)(fun jSigma fe j v)).
apply lc_scalF_compat; rewrite scalF_comm_R.
f_equal; apply extF; intro; apply shape_fun_correct.
rewrite lc_kron_l_in_r; easy.
apply (Sigma_lm fe).
Qed.

End FE_Facts1.

Section FE_Facts2.

Lemma shape_fun_ext :
   {fe1 fe2 : FE} (H : fe1 = fe2),
    shape_fun fe1 =
      castF (eq_sym (f_equal ndof H))
        (mapF (castF_fun (eq_sym (f_equal d H))) (shape_fun fe2)).
Proof. intros; subst; rewrite castF_id castF_fun_id; easy. Qed.

End FE_Facts2.