# A Rocq Formalization of Simplicial Lagrange Finite Elements

Sylvie Boldo, François Clément, Vincent Martin, Micaela Mayero, Houda Mouhcine

► To cite this version:

Sylvie Boldo, François Clément, Vincent Martin, Micaela Mayero, Houda Mouhcine. A Rocq Formalization of Simplicial Lagrange Finite Elements. 2025. hal-05128954

HAL Id: hal-05128954

https://inria.hal.science/hal-05128954v1

Preprint submitted on 25 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A Rocq Formalization of Simplicial Lagrange
Finite Elements

Sylvie Boldo¹*, François Clément²,³, Vincent Martin⁴,
Micaela Mayero¹,⁵, Houda Mouhcine¹,²,³,⁵

¹*Université Paris-Saclay, CNRS, ENS Paris-Saclay, Inria, Laboratoire Méthodes Formelles, 91190 Gif-sur-Yvette, France.

²Centre Inria de Paris, 48 rue Barrault, CS 61534, 75647 Paris Cedex,
France.

³CERMICS, École des Ponts, 77455 Marne-la-Vallée, France.

⁴Université de technologie de Compiègne, LMAC, Compiègne, France.

⁵LIPN, Université Paris 13 - USPN, CNRS UMR 7030, 93430
Villetaneuse, France.

*Corresponding author(s). E-mail(s): sylvie.boldo@inria.fr;

Contributing authors: francois.clement@inria.fr; vincent.martin@utc.fr; mayero@lipn.univ-paris13.fr; houda.mouhcine@inria.fr;

Abstract

The finite elements method is a popular method to numerically solve partial dif-
ferential equations. In the long-term goal of proving its correctness, we focus here
on the formal definition of what is a finite element: a record in the Rocq proof
assistant with both values and proofs of validity, including the main one called
unisolvence. We then instantiate this record with the most popular and useful, the
simplicial Lagrange finite elements for evenly distributed nodes, for any dimension
and any polynomial degree. These proofs require many results (definitions, lem-
mas, canonical structures) about finite families, affine spaces, multidimensional
polynomials, in the context of finite or infinite-dimensional spaces.

**Keywords:** Formal proof, Rocq, Coq, Algebra, Affine Space, Numerical Analysis, Finite Element Method, Lagrange Finite Element
# Contents

<table><tbody><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>1.1</td><td>General Introduction</td><td>4</td></tr><tr><td>1.2</td><td>State of the Art</td><td>6</td></tr><tr><td>2</td><td>Mathematical Presentation of Finite Elements</td><td>8</td></tr><tr><td>2.1</td><td>A Few Mathematical References</td><td>8</td></tr><tr><td>2.2</td><td>Functional Spaces</td><td>9</td></tr><tr><td>2.3</td><td>The Finite Element Method on a Simple Example</td><td>9</td></tr><tr><td>2.3.1</td><td>Continuous Poisson Problem</td><td>9</td></tr><tr><td>2.3.2</td><td>Discrete Poisson Problem</td><td>10</td></tr><tr><td>2.3.3</td><td>Mesh</td><td>11</td></tr><tr><td>2.3.4</td><td>A Conforming Finite Element Space, Log P<sub>k</sub><sup>d</sup></td><td>11</td></tr><tr><td>2.3.5</td><td>Building the Linear System</td><td>12</td></tr><tr><td>2.4</td><td>The Finite Element as a Triple</td><td>13</td></tr><tr><td>2.5</td><td>A Zoology of Finite Elements</td><td>14</td></tr><tr><td>3</td><td>Rocq Support Libraries</td><td>16</td></tr><tr><td>3.1</td><td>Focus on Some Specific Syntax of Rocq</td><td>16</td></tr><tr><td>3.2</td><td>Reals and Coquelicot Libraries</td><td>17</td></tr><tr><td>3.3</td><td>MathComp Library</td><td>17</td></tr><tr><td>3.4</td><td>Classical Aspects</td><td>18</td></tr><tr><td>4</td><td>Subsets and Algebra</td><td>18</td></tr><tr><td>4.1</td><td>Functions and Restrictions</td><td>19</td></tr><tr><td>4.2</td><td>Finite Families</td><td>20</td></tr><tr><td>4.3</td><td>Algebraic Structures</td><td>21</td></tr><tr><td>4.3.1</td><td>Monoid</td><td>21</td></tr><tr><td>4.3.2</td><td>Ring</td><td>23</td></tr><tr><td>4.3.3</td><td>Module Space</td><td>24</td></tr><tr><td>4.3.4</td><td>Algebraic Substructures</td><td>24</td></tr><tr><td>4.3.5</td><td>Kernel</td><td>24</td></tr><tr><td>4.4</td><td>Finite-dimensional Linear Algebra</td><td>25</td></tr><tr><td>4.5</td><td>Binomials</td><td>26</td></tr><tr><td>5</td><td>Affine Algebra</td><td>26</td></tr><tr><td>5.1</td><td>Affine Space and Barycenter</td><td>27</td></tr><tr><td>5.2</td><td>Finite-dimensional Affine Algebra and Barycentric Coordinates</td><td>30</td></tr><tr><td>6</td><td>Formalization of Finite Elements</td><td>31</td></tr><tr><td>6.1</td><td>Rocq Definition as a Record</td><td>32</td></tr><tr><td>6.2</td><td>Degenerate Cases</td><td>33</td></tr><tr><td>7</td><td>Geometric Transformations</td><td>34</td></tr><tr><td>7.1</td><td>Reference Geometry</td><td>34</td></tr><tr><td>7.2</td><td>Lagrange Polynomials and Geometric Transformation</td><td>35</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td>7.2.1</td>
      <td>Reference Lagrange Polynomials</td>
      <td>35</td>
    </tr>
    <tr>
      <td>7.2.2</td>
      <td>Affine Geometric Transformation</td>
      <td>36</td>
    </tr>
    <tr>
      <td>7.2.3</td>
      <td>Lagrange Polynomials</td>
      <td>38</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Geometric Transformations and Face Hyperplanes</td>
      <td>38</td>
    </tr>
    <tr>
      <td>7.3.1</td>
      <td>Face Hyperplanes</td>
      <td>39</td>
    </tr>
    <tr>
      <td>7.3.2</td>
      <td>Injection onto a Face Hyperplane</td>
      <td>40</td>
    </tr>
    <tr>
      <td>7.3.3</td>
      <td>Geometric Transformation Swapping Vertices</td>
      <td>41</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Geometric Transformation of a Finite Element</td>
      <td>43</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Multi-indices and Polynomial Space P<sup>d</sup><sub>k</sub></td>
      <td>44</td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>Multi-indices</td>
      <td>45</td>
    </tr>
    <tr>
      <td>8.1.1</td>
      <td>Formal Definition of Multi-indices</td>
      <td>46</td>
    </tr>
    <tr>
      <td>8.1.2</td>
      <td>Ordering Multi-indices</td>
      <td>48</td>
    </tr>
    <tr>
      <td>8.1.3</td>
      <td>Bijection and Multi-indices</td>
      <td>50</td>
    </tr>
    <tr>
      <td>8.1.4</td>
      <td>Face Multi-indices</td>
      <td>50</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Polynomial Space P<sup>d</sup><sub>k</sub></td>
      <td>51</td>
    </tr>
    <tr>
      <td>8.2.1</td>
      <td>P<sup>d</sup><sub>k</sub> Definition and Basis</td>
      <td>51</td>
    </tr>
    <tr>
      <td>8.2.2</td>
      <td>P<sup>d</sup><sub>k</sub> for Small Values of d and k</td>
      <td>53</td>
    </tr>
    <tr>
      <td>8.2.3</td>
      <td>Factorization of Polynomials</td>
      <td>54</td>
    </tr>
    <tr>
      <td>9</td>
      <td>LogP<sup>d</sup><sub>k</sub> Simplicial Lagrange Finite Element</td>
      <td>55</td>
    </tr>
    <tr>
      <td>9.1</td>
      <td>Simplicial Lagrange Nodes</td>
      <td>56</td>
    </tr>
    <tr>
      <td>9.1.1</td>
      <td>Lagrange Reference and Generic Nodes</td>
      <td>57</td>
    </tr>
    <tr>
      <td>9.1.2</td>
      <td>Lagrange Nodes and Faces</td>
      <td>58</td>
    </tr>
    <tr>
      <td>9.1.3</td>
      <td>Incrementing the Spatial Dimension</td>
      <td>58</td>
    </tr>
    <tr>
      <td>9.1.4</td>
      <td>Incrementing the Degree of Approximation: Lagrange Subnodes</td>
      <td>59</td>
    </tr>
    <tr>
      <td>9.2</td>
      <td>Lagrange Degrees of Freedom</td>
      <td>61</td>
    </tr>
    <tr>
      <td>9.3</td>
      <td>Unisolvence: Double Induction for Nonzero d and k</td>
      <td>62</td>
    </tr>
    <tr>
      <td>9.4</td>
      <td>Degenerate Case k = 0</td>
      <td>63</td>
    </tr>
    <tr>
      <td>9.5</td>
      <td>Construction of the LogP<sup>d</sup><sub>k</sub> Simplicial Lagrange Finite Element</td>
      <td>63</td>
    </tr>
    <tr>
      <td>9.6</td>
      <td>Face Unisolvence</td>
      <td>65</td>
    </tr>
    <tr>
      <td>10</td>
      <td>Conclusion and Perspectives</td>
      <td>66</td>
    </tr>
    <tr>
      <td>10.1</td>
      <td>Contributions</td>
      <td>66</td>
    </tr>
    <tr>
      <td>10.2</td>
      <td>Difficulties</td>
      <td>67</td>
    </tr>
    <tr>
      <td>10.3</td>
      <td>Design Choices</td>
      <td>68</td>
    </tr>
    <tr>
      <td>10.4</td>
      <td>Perspectives</td>
      <td>70</td>
    </tr>
    <tr>
      <td colspan="2">Bibliography</td>
      <td>70</td>
    </tr>
    <tr>
      <td colspan="2">List of Figures and Tables</td>
      <td>76</td>
    </tr>
  </tbody>
</table>
# 1 Introduction

## 1.1 General Introduction

Our long-term goal is to provide greater guarantees for numerical programs simulating physical systems, as they are widely used, e.g. in healthcare, civil engineering, aeronautics, climatology, and astrophysics. In this article we deal with one of the most popular method for solving partial differential equations, the finite element method (FEM). We focus on the formalization of finite elements (FE), and more specifically on the construction of the most commonly used, the simplicial Lagrange FE. This is a necessary first step towards a formal proof of the whole FEM. The mathematical presentation of the FEM and the corresponding references can be found in Section 2. Informally, the FEM is used to approximate the solution of partial differential equations (PDE). These equations are used to model a wide range of problems, including fluid mechanics (Navier-Stokes equations), electromagnetism (Maxwell equations), heat equation (Fourier), mechanics (elasticity), quantum mechanics (Schrödinger and Heisenberg), weather forecasting and aeronautics, among others.

<figure><img src="image_2.png" /><figcaption>**Fig. 1** 3D mesh of a cerebral aneurysm (colors only help to distinguish zones).</figcaption></figure>

Exact solutions of these PDE are generally not computable. To be able to compute an approximate solution, the idea is to discretize, i.e. to cut the problem to be solved into small pieces. As described in Section 2.3.3, we then create a mesh consisting of “simple” elements on which we know how to perform the calculations. We are now moving from a continuous problem to a discrete one. These elements are called “finite elements”. The mesh is therefore made up of these finite elements covering the total area or volume, see Figure 1. To perform approximate calculations on these different elements, a change of coordinate frame is necessary. This involves changing the representation of space from some “current element” to the “reference element”, where calculations are simpler (see Figure 2); for example in 2D, to the right triangle, with sides of length 1. The technique of using them to approximate the solution to the overall problem is the “finite element method”. For example, Figure 1 represents an application of FEM to healthcare. This is a 3D mesh of a cerebral aneurysm used to calculate blood flow inside the artery (using Navier-Stokes equation) and wall mechanics (shell equation, not depicted here), with fluid structure interaction between wall and blood, see for instance [1].
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Reference Lagrange finite elements $\mathcal{L}^{\text{ag}}\mathbb{P}_k^d$ in dimension $d=2$ (top) and $d=3$ (bottom), with degrees of approximation $k=1$ (left) and $k=2$ (right). When $k=1$, the nodes are only the vertices, and when $k=2$, nodes are added in the middle of edges.</figcaption></figure>

As explained in Section 2.4, an FE is not just a geometric shape, it is a triple $(K, P, \Sigma)$. Simply said, $K$ represents the geometry, $P$ represents the set of approximation polynomials, and $\Sigma$ is a way to characterize the polynomials, for instance by the values taken at some nodes (the blue points in Figure 2). Moreover, a “correct” FE needs to satisfy the unisolvence property, which means that the above characterization is unique. There are several types of FE, such as Lagrange, Raviart-Thomas, the choice depending on the equations to be approximated. Representations of simplicial Lagrange FE are given in Figure 2 for space dimensions $d=2$ and 3, and for affine ($k=1$) and quadratic ($k=2$) polynomials (see Section 9 for a formal description of Lagrange FE for evenly distributed nodes).

There are many programs for performing approximation calculations based on the FEM, and many FE libraries are available, such as FreeFEM++,<fnref n="1" /> FELiScE,<fnref n="2" />, and XLiFE++.<fnref n="3" /> As in any program, errors are possible, either unexpected (bugs) or expected such as method errors, as well as errors due to calculations made by the processor, in particular due to floating-point arithmetic. To make these programs safe, our aim is to formally prove the correctness of their implementation, in terms of both method and floating-point errors. To formally prove the correctness of such a program, we need to start by formalizing and proving the necessary mathematical properties in a proof assistant. The proof tool we have chosen is Rocq,<fnref n="4" /> previously called Coq. First, Rocq provides a library of floating-point numbers [2]. Since the 2000s, this tool has enabled us to formalize and prove a 1D wave equation resolution program [3], and also a number of real analysis properties necessary for our goal of proving critical numerical programs, such as the Lax–Milgram theorem [4], the Lebesgue integral [5, 6], and the Bochner integral [7]. These formalizations and the present work are prerequisites to the formalization of the FEM.

<footnote n="1">¹https://freefem.org</footnote>

<footnote n="2">²https://gitlab.inria.fr/felisce/felisce/</footnote>

<footnote n="3">³https://xlifepp/pages.math.cnrs.fr/</footnote>

<footnote n="4">⁴https://rocq-prover.org</footnote>
In Section 2, we start by describing the main mathematical notions for the finite element method. In Section 3, we briefly introduce the Rocq libraries we use, such as Coquelicot and MathComp. In Section 4, we begin to describe our developments in Rocq, starting with subsets and algebra. Section 5 present a corpus of results dedicated to affine algebra. Section 6 presents the definition of an FE as a Rocq record. In Section 7, we present several geometric transformations (i.e. related to K), including some that allow us to transform one FE into another. In Section 8, we focus on a polynomial approximation space (P). Section 9 is dedicated to the application of our formalization to Lagrange FE (with the corresponding Σ), including the proof of unisolvence, before concluding and describing our perspectives in Section 10.

The library, now called rocq-num-analysis, is available as Opam<fnref n="5" /> packages, and the
GIT repository is available at
https://lipn.univ-paris13.fr/rocq-num-analysis/-/tree/2.0/

The Rocq code and documentation is also browsable at
https://lipn.univ-paris13.fr/rocqdoc-num-analysis/2.0/

The pen-and-paper proofs corresponding to the present work, developed for their
formalization, are detailed in [8].

Note also that throughout the paper, annotations such as in “FE” contain a click-
able link towards the corresponding point in our Rocq formalization, be it a definition
(as here), a lemma, or an entire file. Moreover, annotations in pink boxes such as
[m..n] and [n.+1] indicate first references to notations, be they mathematical or formal.

1.2 State of the Art

To our knowledge, finite elements have not yet been formalized elsewhere. However, the ongoing *SciLean* project is quite related to our concerns, and there have also been some developments in the specific areas of geometrical aspects such as mesh construction, polynomial interpolation, and approximate solution to differential equations.

SciLean<fnref n="6" /> is a scientific computing library written in Lean [9]. The goal is to mix in the same framework the implementation of scientific computing algorithms together with the underlying mathematics, and the proofs of correctness of the algorithms. However, as the authors write: “This library is in an *early stage of development*.” It provides effective computations but with many admitted proofs.

The *SciLean* library includes a *FiniteElement/Mesh* directory<fnref n="7" /> that provides a unified representation of various geometric shapes, called *prism*, such as *n*-simplices, *n*-cubes, and pyramids through induction.<fnref n="8" /> Then, the usual concepts are specified by induction on the prism geometry, such as barycentric coordinates, Lagrange basis polynomials, and Lagrange nodes. And this goes up to the construction of a triangular mesh.<fnref n="9" /> However, there does not seem to be any specification of reference element and

<footnote n="5">⁵https://rocq-prover.org/p/rocq-num-analysis/2.0.0</footnote>

<footnote n="6">⁶https://github.com/lecopivo/SciLean/blob/efbb2d0</footnote>

<footnote n="7">⁷https://github.com/lecopivo/SciLean/blob/efbb2d0/SciLean/Modules/FiniteElement/Mesh</footnote>

<footnote n="8">⁸https://github.com/lecopivo/SciLean/blob/efbb2d0/SciLean/Modules/FiniteElement/Mesh/PrismRepr.lean</footnote>

<footnote n="9">⁹https://github.com/lecopivo/SciLean/blob/efbb2d0/SciLean/Modules/FiniteElement/Mesh/TriangularMesh.lean</footnote>
geometric transformation (see Section 7), or of unisolvence and finite element (see
Sections 6 and 9).

About geometry, we can mention the following formalization conducted in
Isabelle/HOL [10]. L.I. Meikle and J.D. Fleuriot [11] focused on enhancing two-
dimensional planar geometry, such as convex hulls and Delaunay triangulations.
This includes basic geometric definitions such as points, and orientation of points
(collinearity, betweenness).

Additionally, Y. Bertot and J.-F. Dufourd [12], focused on the formal verification of an algorithm for constructing Delaunay triangulations of a mesh using the Rocq proof assistant. This algorithm was implemented to construct triangles such that none of the input points lie inside the circumcircle of any triangle. In a more recent work, Y. Bertot extended his earlier work in [13] by addressing a broader class of triangulation algorithms, that guarantees that the resulting mesh satisfies important geometric properties, such as covering the convex hull and maintaining boundary edges. In our development, we do not deal with the building of the mesh. Instead, we specifically focus on a given mesh cell as the convex hull of its affinely independent vertices, and subsequently construct a finite element on this cell.

In Isabelle/HOL, R. Thiemann and A. Yamada address polynomial interpolation, including Lagrange interpolation and the case of integer polynomials,<fnref n="10" /> and, using different HOL algebraic structures, E. Karayel focuses on the case of finite fields.<fnref n="11" /> In Lean, the mathlib library offers a thorough polished study,<fnref n="12" /> by emphasizing the definition of the Lagrange polynomial basis, along with the Lagrange nodes for the interpolation process, supported by various theorems to validate that the Lagrange polynomial exactly match the provided value at each specified point. These three developments are considering univariate polynomials only.

In our Rocq work, we also define the Lagrange polynomial bases and various theorems. The construction is different (we define them from reference Lagrange polynomials, see Section 7.2), but the main difference is that we are in a multivariate setting. We define Lagrange polynomial bases for $d = 1$ (with any degree), and for $k = 1$ (and any dimension). This suffices for our needs.

As stated in Section 1.1, the FEM is essentially used to approximate PDE solutions. The formalization of ordinary differential equations (ODE) and PDE is therefore also of interest. There have been various formal developments in these domains. In Isabelle/HOL, F. Immler and J. Hölzl [14] present a formalization of ODE and the verification of rigorous numerical algorithms. In Lean, ODE developments are also in progress<fnref n="13" />. In Mizar, S. Otsuki et al. [15] formalize a simple example for linear PDE and its solution using the method of separation of variables. In HOL, E. Deniz et al. [16] use the separation of variables method for formally verifying the solution of the PDE, which allows the modeling the heat transfer in the slab under various initial and boundary conditions. In Rocq, M. Tekriwal et al. [17] formally prove the Lax

<footnote n="10">¹⁰https://www.isa-afp.org/entries/Polynomial_Interpolation.html</footnote>

<footnote n="11">¹¹https://www.isa-afp.org/entries/Interpolation_Polynomials_HOL_Algebra.html</footnote>

<footnote n="12">¹²https://leanprover-community.github.io/mathlib4_docs/Mathlib/LinearAlgebra/Lagrange.html</footnote>

<footnote n="13">¹³https://github.com/lecopivo/SciLean/blob/efbb2d0/SciLean/Analysis/ODE</footnote>
equivalence theorem for the finite difference method, which generalizes a part of the
mathematical proof of [3].

## 2 Mathematical Presentation of Finite Elements

In order to present the mathematics of finite elements (FE), we start by introducing
the principle of the finite element method (FEM) in a simple, standard example. Note
that the present formalization work does not include the full FEM; it focuses only on
the construction of finite elements. We limit this presentation to the Poisson problem,
that is Problem (1), in order to exemplify a typical stationary problem. Note that
many extensions of the FEM (such as Discontinuous Galerkin methods) derive from
this framework. However, they are outside the scope of this study. We do not study
either the Finite Volume Methods. The goal of this section is to explain and motivate
our work and to put our formalization into perspective with respect to the full FEM
and to the applied mathematics literature.

After giving a few references in Section 2.1, we set some notations on functional spaces in Section 2.2. The FEM is described in Section 2.3: the continuous Poisson formulation is presented, then discretized with a given mesh of the domain. This done, in Section 2.3.4, we give a first introduction of the most common FE, the Lagrange finite element over a simplex, that are fully formalized in this paper. This allows to write the linear system whose solution is the approximation we are looking for. This presentation concludes with the mathematical definition of a finite element as a triple in Section 2.4 (whose formalization is the object of Section 6), and a short overview of various FE in Section 2.5.

Note that all sections but 2.3.4 and 2.4 can be skipped if the reader is less interested
by the broader picture or wants to focus on the Rocq formalization.

### 2.1 A Few Mathematical References

Among the vast literature about the FEM, we can cite a few references: P.G. Ciar-
let [18], I. Babuška and T. Strouboulis [19], A. Quarteroni and A. Valli [20], A. Ern
and J.-L. Guermond [21–24], S.C. Brenner and L.R. Scott [25], and O.C. Zienkiewicz,
R.L. Taylor and J.Z. Zhu [26]. The first part of the Handbook of Numerical Analysis
is dedicated to the finite element methods [27–31].

The present paper is mostly based on detailed pen-and-paper proofs of [8], that
were mainly derived from [22]. The parts on algebra and affine geometry are inspired
by B. Gostiaux [32, 33].

In reference works on FE, such as [22], the building of simplicial Lagrange FE can
be based on more powerful tools than strictly necessary, such as differential calculus,
e.g. to characterize affine functions as those with zero second order derivative. This is
mainly because those tools are put to good use afterwards to build other FE in more
complex settings (such as Raviart-Thomas FE), or for the computation of integrals on
the FE. But in fact, simplicial Lagrange FE are essentially a matter of affine proper-
ties, and it is possible to work with almost no tools of analysis, and to employ purely
affine arguments, only based on barycenters, for almost all steps of their construction.
And this is done without having to deal with the linear part of affine objects, or explicitly going down to the affine coordinates, except for a few very simple cases. Hence the consequent corpus on affine algebra developed in Section 5. However, derivation remains a convenient way to establish the linear independence of monomials in the space of multivariate polynomials of given maximum total degree, but this is no affine property. And obviously, in the next steps, when using these FE or extending to other finite elements, notions such as integration, differentiation, and normed spaces will become necessary.

## 2.2 Functional Spaces

Let us give some notations and definitions used throughout this article, including the functional spaces on which we look for the solution to the PDE.

In this paper, the vector quantities in mathematical expressions are written with **bold type**. For instance, this is the case of points $\mathbf{x} \stackrel{\text{def.}}{=} (x_0, \dots, x_{d-1})$ in $\mathbb{R}^d$, or of the gradient $\nabla u \stackrel{\text{def.}}{=} (\frac{\partial u}{\partial x_0}, \dots, \frac{\partial u}{\partial x_{d-1}})^T$ of a scalar function $u \in \mathcal{F}(\mathbb{R}^d, \mathbb{R})$, i.e. from $\mathbb{R}^d$ to $\mathbb{R}$.

In this section, we assume that the space dimension $d$ lies in $\{1, 2, 3\}$. Let $\Omega \subset \mathbb{R}^d$ be an open bounded connected domain, and let $\Gamma \stackrel{\text{def.}}{=} \partial\Omega$ be its regular boundary. In the sequel, for the sake of simplicity, we assume that $\Omega$ is polyhedral.

The most commonly used functional spaces are

$$L^2(\Omega) \stackrel{\text{def.}}{=} \{f : \Omega \to \mathbb{R} \mid \|f\|_0 < +\infty\}, \\ H^1(\Omega) \stackrel{\text{def.}}{=} \{f \in L^2(\Omega) \mid \nabla f \in (L^2(\Omega))^d\}, \quad \text{and} \quad H_0^1(\Omega) \stackrel{\text{def.}}{=} \{f \in H^1(\Omega) \mid f|_{\Gamma} = 0\},$$

where $L^2(\Omega)$ is the Lebesgue space of measurable, square integrable, real functions over $\Omega$, it is equipped with the $L^2$-norm defined by $\|f\|_0 \stackrel{\text{def.}}{=} \sqrt{\int_{\Omega} |f|^2}$. $H^1(\Omega)$ is the Sobolev space made up of $L^2(\Omega)$ functions that have a gradient in $(L^2(\Omega))^d$, and $H_0^1(\Omega)$ is the subspace of functions vanishing on the boundary $\Gamma$. They are equipped with the seminorms defined by $\|f\|_{1,0} \stackrel{\text{def.}}{=} \|\nabla f\|_0$, and $\|f\|_1 \stackrel{\text{def.}}{=} \sqrt{\|f\|_0^2 + \|f\|_{1,0}^2}$. The spaces $(L^2(\Omega), \|\cdot\|_0)$, $(H^1(\Omega), \|\cdot\|_1)$, and $(H_0^1(\Omega), \|\cdot\|_{1,0})$ are Hilbert spaces (i.e. complete normed vector spaces, with an inner product, e.g. see [34, Cor. IX.19, p. 174]).

## 2.3 The Finite Element Method on a Simple Example

Let us present the finite element method on a simple example, the Poisson problem.

### 2.3.1 Continuous Poisson Problem

Let $f \in L^2(\Omega)$. The Poisson problem is defined in its strong form as

$$\text{find } u, \text{ such that } \begin{cases} -\Delta u = f & \text{in } \Omega \\ u = 0 & \text{on } \Gamma, \end{cases} \qquad (1)$$
where the Laplace operator $\Delta$ is defined by $\Delta u \stackrel{\text{def.}}{=} \sum_{i=0}^{d-1} \frac{\partial^2 u}{\partial x_i^2}$. If we multiply by a regular test function $v$, integrate over $\Omega$, integrate by parts, and take the boundary condition into account, we obtain the weak form of the Poisson equation,

$$ \text{find } u \in H_0^1(\Omega), \text{ such that } \forall v \in H_0^1(\Omega), \quad \int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega} f v. \qquad (2) $$

Setting $a(u, v) \stackrel{\text{def.}}{=} \int_{\Omega} \nabla u \cdot \nabla v$, $l(v) \stackrel{\text{def.}}{=} \int_{\Omega} f v$, and $V \stackrel{\text{def.}}{=} H_0^1(\Omega)$, problem (2) takes the general algebraic form

$$ \text{find } u \in V, \text{ such that } \forall v \in V, a(u, v) = l(v). \qquad (3) $$

The Lax–Milgram theorem states that problem (3) admits a unique solution, assuming some hypotheses: continuity and coercivity of $a$, continuity of $l$, and that $(V, \|\cdot\|_V)$ is a Hilbert space. This theorem has already been formalized in Rocq [4].

Note that more general results exist, taking different, Banach spaces for the test function and the solution (Petrov–Galerkin framework and Banach–Nečas–Babuška theorem, e.g. see [21, Chap. 2]).

The Lax–Milgram hypotheses are fulfilled in the case of (2) (e.g. see [27, Chap. 1, pp. 35–36]), proving that the weak Poisson formulation admits a unique solution. It is called the continuous solution (i.e. living in a continuous space), in contrast to the discrete solution (living in a finite dimensional space) that is constructed in the sequel.

### 2.3.2 Discrete Poisson Problem

The conforming Galerkin method consists in replacing $V$ by a finite-dimensional space $V_h \subset V$. Assuming that the bilinear and linear forms remain unchanged, the discrete counterpart of (3) reads

$$ \text{find } u_h \in V_h, \text{ such that } \forall v_h \in V_h, a(u_h, v_h) = l(v_h). \qquad (4) $$

The Lax–Milgram theorem also applies to (4), and Céa’s Lemma provides an estimate on the error,

$$ \|u - u_h\|_V \le C \inf_{v_h \in V_h} \|u - v_h\|_V, $$

which stresses the importance of interpolating $u$, i.e. finding a “good” approximation space $V_h$ to be as close as possible to $u$. The finite element method is a way of constructing such $V_h$, for wide ranges of problems.

We carry on with the discrete version of the simple model (2), that reads

$$ \text{find } u_h \in V_h, \text{ such that } \forall v_h \in V_h, \quad \int_{\Omega} \nabla u_h \cdot \nabla v_h = \int_{\Omega} f v_h. \qquad (5) $$
### 2.3.3 Mesh

As in Figure 1, a mesh $\mathcal{T}_h$ of $\Omega$ is a collection of closed subsets of $\Omega$ called mesh cells or geometric elements $K$, having nonempty interior, disjoints interiors, and covering $\bar{\Omega}$,

$$ \forall K \in \mathcal{T}_h, \mathring{K} \neq \emptyset, \quad (\forall K, K' \in \mathcal{T}_h, K \neq K' \implies \mathring{K} \cap \mathring{K}' = \emptyset), \quad (6) \\ \text{and } \bar{\Omega} = \bigcup_{K \in \mathcal{T}_h} K. $$

The index $h > 0$ represents the degree of refinement of the mesh. It is defined as the maximum diameter of all the mesh cells. Often, the mesh is made up of triangles or tetrahedra, or of quadrangles or hexahedra.

In the sequel, we assume that all the mesh cells are the image of a reference polyhedron $\hat{K}$ (such as the unit rectangular simplex) by a smooth invertible function, the geometric mapping. Let $(\hat{v}_i)_{i \in [0..n_{vtx})}$ be the $n_{vtx}$ vertices of $\hat{K}$, and $(v_i)_{i \in [0..n_{vtx})}$ be the ones of a cell $K$, often called the current cell, where for all $m, n \in \mathbb{N}$, the notation $[m..n]$ represents the interval of integers $[m, n) \cap \mathbb{N}$. Then, the geometric mapping is often set as a combination using $n_{vtx}$ shape functions,

$$ T^{\nu}(\hat{\mathbf{x}}) \stackrel{\text{def.}}{=} \sum_{i \in [0..n_{\text{vtx}})} \hat{\mathcal{L}}_i(\hat{\mathbf{x}})v_i, \quad (7) $$

where $\hat{\mathcal{L}}_i$ is a polynomial in $\mathbb{R}^d$ that satisfies for all $i, j \in [0..n_{vtx})$, $\hat{\mathcal{L}}_i(\hat{v}_j) = \delta_{ij}$, with $\delta$ being the Kronecker delta function (see Section 7.2.1 in the simplicial case). Thus, $T^{\nu}(\hat{v}_i) = v_i$. Let $J^{\nu}$ be the invertible Jacobian matrix of $T^{\nu}$.

When the geometric mapping is affine, the current cell $K$ is a polyhedron. In other cases, the faces may not be planar and $K$ may not be convex. See an example of geometric mapping in Figure 4 (p. 36).

Moreover, the mesh is assumed to be conforming (or matching), i.e. the intersection in dimension ($d-1$) of two neighboring cells is a full face. For instance, when $d=3$, the intersection of two distinct mesh cells $K$ and $K'$ is either empty, a vertex, a full edge, or an entire face for $K$ and for $K'$.

### 2.3.4 A Conforming Finite Element Space, $\mathcal{L}^{\text{ag}}\mathbb{P}_k^d$

Let us take an example to set the scene, the simplicial Lagrange FE $\mathcal{L}^{\text{ag}}\mathbb{P}_k^d$. Note that this brief presentation is detailed in the rest of the article: the formal definition of $\mathcal{L}^{\text{ag}}\mathbb{P}_k^d$ is the object of Sections 7, 8, and 9. Recall also that the vector quantities in mathematical expressions are written with bold type. For instance, this is the case of multi-index $\alpha$, monomials functions $X^\alpha$, and the family of linear forms $\Sigma^{\mathcal{L}^{\text{ag}}\mathbb{P}_k^d}$ below.

Assume that the conforming mesh $\mathcal{T}_h$ is made up of simplices. Let $K \in \mathcal{T}_h$. From (6), $K$ has a nonempty interior. Let $k \in \mathbb{N}$. Let the multi-index set be

$$ \mathcal{A}_k^d \stackrel{\text{def.}}{=} \{ \alpha \in \mathbb{N}^d \mid |\alpha| \le k \}, \quad (8) $$
where $|\alpha|^{\text{def.}} = \sum_{i \in [0..d)} \alpha_i$. For $\alpha \in \mathcal{A}_k^d$, let the monomial functions be

$$ \mathbf{X}^{\alpha} \stackrel{\text{def.}}{=} \left( \mathbf{x} \in \mathbb{R}^d \longmapsto \prod_{i \in [0..d)} x_i^{\alpha_i} \right). \qquad (9) $$

The space of polynomials of degree at most $k$ of $d$ variables is defined by

$$ \mathbb{P}_k^d \stackrel{\text{def.}}{=} \operatorname{span}(\mathbf{X}^{\alpha})_{\alpha \in \mathcal{A}_k^d} \qquad (10) \\ = \left\{ \left( \mathbf{x} \longmapsto \sum_{\alpha \in \mathcal{A}_k^d} c_{\alpha} x_1^{\alpha_1} \dots x_d^{\alpha_d} \right) : \mathbb{R}^d \to \mathbb{R} \;\middle|\; (c_{\alpha})_{\alpha \in \mathcal{A}_k^d} \in \mathbb{R} \right\}. $$

Let $(\mathbf{a}_{\alpha})_{\alpha \in \mathcal{A}_k^d}$ be pairwise distinct points in $K$, called the nodes of $K$ (fully defined in Section 9.1). Let $\Sigma^{\text{Log}\mathbb{P}_k^d} \stackrel{\text{def.}}{=} (\sigma_{\alpha})_{\alpha \in \mathcal{A}_k^d}$ be the family of linear forms on $\mathcal{F}(\mathbb{R}^d, \mathbb{R})$ such that

$$ \forall \alpha \in \mathcal{A}_k^d, \forall p \in \mathcal{F}(\mathbb{R}^d, \mathbb{R}), \quad \sigma_{\alpha}(p) \stackrel{\text{def.}}{=} p(\mathbf{a}_{\alpha}) \in \mathbb{R}. $$

The linear forms represent somewhat the type of unknowns one wants to compute. In this case of *nodal FE*, $u_h$ is represented by its values at the nodes, see other FE in Section 2.5.

Then, the global finite element space $V_h$ is defined as the space of continuous functions whose restriction to any mesh cell is in $\mathbb{P}_k^d$,

$$ V_h \stackrel{\text{def.}}{=} \{ v_h \in C^0(\mathbb{R}^d, \mathbb{R}) \mid \forall K \in \mathcal{T}_h, v_{h|K} \in \mathbb{P}_k^d \}. $$

The global continuity condition ensures that $V_h \subset V$. Such a space is called *Lagrange finite element space*. It can be shown that it suffices to define $V_h$ element-wise, the continuity over $\Omega$ is recovered by applying discrete inter-element conditions (see face unisolvency in Section 9.6). For computational efficiency, it is convenient to build a basis for $V_h$ where each basis function has a small support. Typical FE basis functions are nonzero only on elements sharing a face, an edge, or a vertex.

### 2.3.5 Building the Linear System

Let $(\phi_i)_{i \in [0..N_{\text{dof}})}$ be a basis of $V_h$, where $N_{\text{dof}}$ is the total number of “degrees of freedom” (i.e. the dimension) of $V_h$. Take successively $v_h = \phi_i$ in (5) and write

$$ \text{find } u_h = \sum_{j \in [0..N_{\text{dof}})} U_j^h \phi_j \in V_h, \text{ such that} \\ \qquad \forall i \in [0..N_{\text{dof}}), \quad \sum_{j \in [0..N_{\text{dof}})} \left( \int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i \right) U_j^h = \int_{\Omega} f \phi_i, $$
which can be seen as the linear system

$$
\begin{array}{ll}
\text{find } U^h \stackrel{\text{def.}}{=} [U_0^h, \dots, U_{N_{\text{dof}}-1}^h]^\top \in \mathbb{R}^{N_{\text{dof}}}, \text{ such that } A U^h = b, & \\
\text{where } \left\{ \begin{array}{ll} \forall i, j \in [0..N_{\text{dof}}), & A_{ij} \stackrel{\text{def.}}{=} \int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i, \\ \\ \forall i \in [0..N_{\text{dof}}), & b_i \stackrel{\text{def.}}{=} \int_{\Omega} f \phi_i. \end{array} \right. & (11)
\end{array}
$$

It then suffices to build A and b, and apply a linear solver. To compute A, one usually approximates the integrals on the reference cell $\hat{K}$, using a change of variable and a quadrature formula (at $N_q$ points $(\hat{\xi}_{i_q})_{i_q \in [0..N_q]}$) with weights $(\hat{\omega}_{i_q})_{i_q \in [0..N_q]}$),

$$
\begin{align}
\int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i &= \sum_{K \in \mathcal{T}_h} \int_K \nabla \phi_j \cdot \nabla \phi_i = \sum_{K \in \mathcal{T}_h} \int_{\hat{K}} (\mathbf{J}^{\mathbf{v}})^{-T} \hat{\nabla} \hat{\phi}_j \cdot (\mathbf{J}^{\mathbf{v}})^{-T} \hat{\nabla} \hat{\phi}_i |\det \mathbf{J}^{\mathbf{v}}| \nonumber \\
&\approx \sum_{K \in \mathcal{T}_h} \sum_{i_q \in [0..N_q]} \left( (\mathbf{J}^{\mathbf{v}})^{-T} \hat{\nabla} \hat{\phi}_j \cdot (\mathbf{J}^{\mathbf{v}})^{-T} \hat{\nabla} \hat{\phi}_i |\det \mathbf{J}^{\mathbf{v}}| \right) (\hat{\xi}_{i_q}) \hat{\omega}_{i_q}. \tag{12}
\end{align}
$$

Note that for an invertible matrix $M$, $M^{-T}$ represents the transposed of its inverse. Formula (12) is used in practice to build the matrix $A$ piece by piece. For each element, taking advantage of the small support of the $\phi_i$'s, one computes small matrices containing the nonzero $\int_K \nabla \phi_j \cdot \nabla \phi_i$ terms (it is done in the reference element), then one fills the corresponding components of the matrix. This procedure is called the *assembly* of the matrix, it outputs a sparse matrix $A$. Once $A$ is built, it remains to solve the (sparse) linear system (11) using one's favorite library, to compute $U^h$, and then to reconstruct $u_h$.

We presented the finite element method on a simple but typical example, from the continuous variational formulation to the construction of the linear system that needs to be solved. This is done with numerous mesh cells, on which FE are set. It is time to present the mathematical definition of a finite element. From now on, we focus on a single element.

## 2.4 The Finite Element as a Triple

A finite element depends on three components: a geometric element $K$, a discrete approximation space $P$ where the approximated solution $u_{h|K}$ is sought (typically a finite dimensional polynomial space), and a finite set of linear forms $\Sigma$ that specify the discrete quantities that are primarily computed, and that are used to reconstruct $u_{h|K}$. The main property of a FE is that these values of the linear forms uniquely determine an approximation function in $P$.

The construction of a FE relies on several values. The integer $d$ is the spatial dimension of the domain. The unknown function $u$ can be either scalar or vector, i.e. belonging to $\mathcal{F}(\mathbb{R}^d, \mathbb{R}^q)$, the space of functions from $\mathbb{R}^d$ to $\mathbb{R}^q$, where $q$ is generally 1 or $d$. The approximation space $P$ is often characterized by a degree $k$, representing
the accuracy of the approximation space, e.g. the total degree of polynomial functions in (10). And finally, $[n_{\text{dof}}]$ is the number of linear forms in $\Sigma$.

Following [18] and [22, Sec. 5.2, pp. 50–52], a finite element is mathematically defined as a triple $(K, P, \Sigma)$, where:

• $K$ is the *geometric element*, i.e. a closed bounded subset of $\mathbb{R}^d$ with nonempty interior. Typically, $K$ is a nondegenerate $d$-polytope (the generalization of points, segments, polygons or polyhedra to dimension $d$), or the image of such a subset by a regular function.

• $P$ is the *approximation space*, i.e. a nonzero, finite-dimensional vector space of functions from $\mathbb{R}^d$ to $\mathbb{R}^q$. It is meant to be a space of polynomial functions, or the image of such a space by a regular function.

• $\Sigma \stackrel{\text{def.}}{=} (\sigma_i)_{i \in [0..n_{\text{dof}})}$ is the collection of $n_{\text{dof}}$ *local degrees of freedom*, i.e. $n_{\text{dof}}$ linear forms on $P$, such that the application $[\Phi_\Sigma] \stackrel{\text{def.}}{=} (p \in P \mapsto (\sigma_i(p))_{i \in [0..n_{\text{dof}})} \in \mathbb{R}^{n_{\text{dof}})}$ is an isomorphism (a linear and bijective application between $P$ and $\mathbb{R}^{n_{\text{dof}}}$).

The bijectivity of $\Phi_\Sigma$ is called *unisolvence*. It is easy to see that unisolvence requires that the dimension of $P$ is equal to $n_{\text{dof}}$. The unisolvence property means that an approximation function of $P$ is uniquely characterized by the $n_{\text{dof}}$ discrete values $u \in \mathbb{R}^{n_{\text{dof}}}$ such that for all $i \in [0..n_{\text{dof}})$, $\sigma_i(p) = u_i$. For instance, in the case of nodal FE, see Section 2.3.4, the values at the nodes defines the approximation polynomial.

Depending on the hypotheses on the geometric transformation from the reference to the current element $T^\nu$ (e.g. see (7)), $K$ can be a $d$-polytope, i.e. with planar hyperfaces, or not, and can be convex, or not. Not all geometric elements are possible, though. In this document, we assume that the mesh is *affine*, which means that $T^\nu$ is supposed to be invertible and affine. The non-affine meshes are more tricky to handle, notably because $T^\nu(P)$ may not be included in $P$, creating a discrepancy between the reference and current elements. It is also harder to check the conformity of the mesh and to compute (12). Note that the geometric mapping is affine if and only if each of its components is in $\mathbb{P}_1^d$. With a slight abuse of notation, the geometric element is commonly labeled $\mathbb{P}_1^d$, $\mathbb{P}_2^d$, or $\mathbb{Q}_1^d$ (see below), when the components of $T^\nu$ are in the corresponding polynomial space.

In this work, the FE triple is specified in Section 6, using the hypothesis of an affine mesh, thus $K$ is assumed to be a $d$-polytope.

## 2.5 A Zoology of Finite Elements

There exists several families of finite elements. Note that, except the simplicial Lagrange FE for evenly distributed nodes described in this article, none of these are formalized in Rocq yet, but most fit into the generic framework of Section 6. The object of this section is merely to show that there exists a large variety of FE.

Nodal FE are the most popular. Their degrees of freedom (i.e. the linear forms $\Sigma$) are the evaluation operators at given points called nodes, and the discrete solution is reconstructed from these values. Among them, the case of evenly distributed (i.e. Lagrange) nodes, leaning on the vertices of the geometric element $K$, is very common.
For a degree $k \in \mathbb{N}$, they are denoted here as $^{\text{Log}}\mathbb{P}_k^d$ for $d$-simplices, and $^{\text{Log}}\mathbb{Q}_k^d$ for $d$-dimensional hexahedra (see [22, Sec. 7.4 and 6.4]). There exist also the Hermite FE, where the linear forms are the evaluation at the nodes of both the function and its derivatives (e.g. see [35]), and nodal FE for which the nodes are not evenly distributed (e.g. see [22, Sec. 6.3.5], and Remark 7.14)].

For some (non-nodal) families of FE, the degrees of freedom are average values or moments in the geometric element $K$. For instance, the Crouzeix-Raviart FE on simplices is based on $\mathbb{P}_1^d$ polynomials, and its degrees of freedom are the mean values on the faces of $K$, inducing discontinuous functions over $\Omega$ (see [22, Sec. 7.5]).

For some families where the approximated solution $\mathbf{u}_{h|K}$ is $d$-vectorial, the normal fluxes through the faces of $K$ are computed, and in others it is the tangential integral along the edges. Examples of face-flux driven FE are the Raviart-Thomas ($\mathbb{RT}$), Brezzy-Douglas-Marini ($\mathbb{BDM}$), and Brezzy-Douglas-Fortin-Marini ($\mathbb{BDFM}$) families, both for $d$-simplices, and regular $d$-hexahedra (see [22, Chap. 14]). They aim at approximating $\mathbf{H}(\text{div})$ functions, that is vector functions $\mathbf{u}$ in $(L^2)^d$, whose divergence $\nabla \cdot \mathbf{u} \stackrel{\text{def}}{=} \sum_{i \in [0..d)} \frac{\partial u_i}{\partial x_i}$ is in $L^2$. For instance, for a tetrahedron $K$, whose faces are denoted by $(F_i)_{i \in [0..3]}$ and have a normal unit vector $\mathbf{n}_i$, the space $\mathbb{RT}_0(K)$ is represented by the quantities $\frac{1}{|F_i|} \int_{F_i} \mathbf{u}_h \cdot \mathbf{n}_i$, where $|F_i|$ is the area of $F_i$. The Nedelec family is an example of FE approximating $\mathbf{H}(\text{curl})$ functions, that is vector functions $\mathbf{u}$ in $(L^2)^d$, whose curl $\nabla \times \mathbf{u}$ is in $L^2$, where $\times$ is the cross product. Its degrees of freedom are the tangential component on the edges (both for $d=2$ or $3$).

Note that the term “Lagrange” can be used for several distinct notions. First, given $n+1$ distinct points $(\mathbf{a}_i)_{i \in [0..n]}$ in $\mathbb{R}^d$, the Lagrange polynomials $(\mathcal{L}_i)_{i \in [0..n]}$ denote the family of polynomials that satisfies $\mathcal{L}_i(\mathbf{a}_j) = \delta_{ij}$. Under some hypotheses, they constitute a basis of some finite dimensional polynomial space, called the *Lagrange basis*. Then, in the context of FE, the points $(\mathbf{a}_i)_{i \in [0..n]}$ are called *Lagrange nodes*. Finally, the *Lagrange finite element* is $^{\text{Log}}\mathbb{P}_k^d$ or $^{\text{Log}}\mathbb{Q}_k^d$, according to the shape of $K$ (see Section 2.3.4 for the simplicial case).

Note also that in the literature, the term “finite element (FE)” often refers to different notions. The genuine meaning refers to a triple $(K, P, \mathbf{\Sigma})$ that specifies how functions are approximated (in $P$) and represented (with $\mathbf{\Sigma}$) on a cell $K$. The mesh cell $K$ can be called (*geometric finite element*, or simply *element*). When $K$ and $\mathbf{\Sigma}$ are unambiguous, the approximation space $P$ is often assimilated with the FE. For instance, the Lagrange FE $^{\text{Log}}\mathbb{P}_k^d$ is often denoted in the literature as the polynomial space $\mathbb{P}_k^d$. Finally, the *finite element method (FEM)* refers to a methodology to solve some partial differential equations, as described in Section 2.3.

In this article, we try and be as clear as possible, and distinguish the different notions. Thus, the FE always refers to the triple. In this respect, $^{\text{Log}}\mathbb{P}_k^d$ denotes the $d$-simplicial Lagrange FE of degree $k$, whereas $\mathbb{P}_k^d$ denotes the space of polynomials of total degree at most $k$.
# 3 Rocq Support Libraries

As stated in Section 1, this work fits in a larger framework of formally verified real analysis results and programs, mainly with Rocq. A full description of the Rocq system is beyond the scope of this paper. We just remind that it is a proof assistant, based on the calculus of inductive constructions [36, 37], and therefore on type theory and intuitionistic logic by default.

Without going into details, we just describe in Section 3.1 a few Rocq-specific key-
words that are often used. Moreover, for readers unfamiliar with the concrete syntax of
Rocq, we explain unusual notations as and when required. Then, Sections 3.2 and 3.3
present the Rocq libraries that we use, and Section 3.4 lists the axioms we rely on.

## 3.1 Focus on Some Specific Syntax of Rocq

A salient point is that we make extensive use of implicit arguments. An implicit argu-
ment of a function is an argument which can be inferred from contextual knowledge.
For example, if a function input has type $\mathbb{R}^d$, then the integer $d$ can be implicit. In
Rocq, implicit arguments are enclosed in curly braces, as we can see for example in
the definition of `concatnF` and described at the end of Section 4.3.1.

Curly braces are not just for implicit types. They are also used to distinguish between two or and two exists. For example, ({ A } + { B }) : Set and (A ∨ B) : Prop both mean “A or B”. The difference lies in the computational content of the first one. Similarly, { x | A x } : Set and (∃ x, A x ) : Prop both mean that “there exists an x such that A x holds”. Most of the time, we do not need any computational content, except to decide on an equality or to use a witness in a definition. See, for example, Sections 3.2, 3.4 and 4.1.

In Rocq, there are several ways of declaring types and functions, depending on their intended use. In the following, we briefly outline the differences between the main keywords we are using.

* **Inductive** allows defining types inductively, with possible recursive calls and generation of an induction scheme.

* **Variant** is similar to **Inductive** but without recursive calls, and no induction scheme is generated.

* **Definition** defines a term in the global environment.

* **Fixpoint** allows defining functions by pattern matching over inductive objects using a fixed point construction.

* **Let** is similar to **Definition** but is for local definitions. When the section is closed, all persistent definitions and theorems within it that depend on the constant are wrapped with a let with the same declaration.

* **Hypothesis** and **Variable** are synonymous, they postulate in the environment of a section that a type is inhabited (it can be a property).

* **Context** declares variables in the context of the current section, like **Variable**, but also allows implicit variables.
## 3.2 Reals and Coquelicot Libraries

Since we are formalizing real analysis dedicated to problems in the field of numerical analysis for and by numerical analysts, all our statements, properties and theorems are classical (in the sense of classical logic, and not intuitionistic).

For real numbers, we use the Coquelicot [38] library. This library is an extension of the Reals [39] standard library. The Reals library was developed in the early 2000s. It was a classical axiomatization. In particular, it contains a decidability property of the total order and a multiplicative inverse operation represented as a total function,

total_order_T : ∀ r1 r2 : ℝ, { r1 < r2 } + { r1 = r2 } + { r2 < r1 } .
Rinv : ℝ → ℝ.
Rinv_r : ∀ r, r ≠ 0 → r * Rinv r = 1.

This part of the standard library provides a number of definitions and properties of analysis such as limit, univariate derivative, and Riemann integral. In the early 2020s, Reals was modified [40] to include a constructive common part, inspired by C-CoRN [41], but the backwards compatibility of the classical part was retained and is what we and Coquelicot rely on.

Coquelicot was developed in the mid-2010s and provides a hierarchy of nested algebraic structures, see [38] and Section 4.3 for more details, and [42] about Rocq canonical structures. By assuming functional extensionality, we provided canonical structures so that functions fit into these algebraic structures. This was done in our development for the proof of the Lax-Milgram theorem in LM linear_map, see [4]. We rely on the hierarchy of algebraic structures of the Coquelicot library at the base of which we contributed to add the commutative monoid structure AbelianMonoid (available from version 3.4.0). We also worked on substructures, addressed in Section 4.3.4.

## 3.3 MathComp Library

We also use parts of the Mathematical Components (MathComp) [43] library and a small part of its analysis companion (MathComp-Analysis) [44], and some SSReflect [45] syntax, although most of our formalization is “vanilla Rocq”. For example from Math-Comp, we use injective, bijective, and involutive properties from ssrfun. We use the MathComp notations `n.+1` for the Peano successor `S n` and `n.-1` for the predecessor. We also use the finite dependent type ordinal, “bigops”, and binomial coefficients (see Section 4). From MathComp-Analysis, the only part we use is boolp, for logical results and axioms related to classical logic (see Section 3.4 below). We explain other notations specific to SSReflect and MathComp as and when required.

Our formalization is aimed at both mathematicians (applied mathematics) and logicians (use of theorem provers). The use of a syntax and concepts that are simply accessible for students and new users is therefore also one of our objectives [46, 47].

We have neither embraced all of MathComp, nor all of MathComp-Analysis, and we have developed some algebra we needed and kept the analysis we already developed. This is mainly because of the incompatibility of the MathComp/MathComp-Analysis and Coquelicot canonical structures, used for the representation of algebraic structures.
As explained in Section 3.2 and in more details in Section 4.3, we rely on canonical structures for nested algebraic structures. For instance, lemmas or notations such as `+` can be applied at the same time to monoid elements, vector spaces elements, and functions on $\mathbb{R}^d$, with `Rocq` being able to instantiate the needed algebraic structure. Unfortunately, the hierarchies of `MathComp/MathComp-Analysis` and `Coquelicot`, although similar, cannot be both used in the same development. As we have already built the Hilbert spaces needed for the Lax-Milgram theorem [4] upon the `Coquelicot` hierarchy, we decided not to switch our base, even if it sometimes means similar formalized mathematical results (on Lebesgue integral for example). Our goal is to have a less abstract library, driven by our numerical analysis applications (especially the finite element method for the moment), therefore the Lebesgue integral, Hilbert spaces, and vector subspaces of functions. Moreover, in the future we will need floating-point numbers (`Flocq` library [2]), which is built on top of the `Reals` standard library.

Among the `MathComp` developments, there is also a formalization of multinomials based on `MathComp` polynomials.<fnref n="14" /> Unfortunately, there are not enough results for our purpose, in particular the dimension of the subspace. We found easier to rely on the multi-indices of Section 8.1 which are also used for the Lagrange nodes in Section 9.1 that allow us to explicit a basis of $\mathbb{P}_k^d$.

## 3.4 Classical Aspects

Let us explicit the axioms we rely on, especially classical logic. Indeed, the formalization and proof of the Lax-Milgram theorem [4] has shown us that avoiding the classical axioms can lead to unnecessarily complicated statements for this kind of mathematics. As one objective is to convince mathematicians, we assume classical logic and other useful axioms. We have also developed in `rocq-num-analysis` many easy lemmas dedicated to logical facts in this framework in the file `logic_compl`².

We require the `Reals` standard library [39], which is classical (as opposed to intuitionist). We require `boolp` from the `MathComp-Analysis` [44] library, that provides propositional and functional extensionality axioms, proof irrelevance, and constructive indefinite description. We rename the latter `ex_EX`²,

Check `ex_EX`.
: $\forall (A : \text{Type}) (P : A \to \text{Prop}), (\exists x, P x) \to \{ x : A \mid P x \}.$

Finally, we require `ClassicalChoice` from the standard library, that provides (decidable) excluded middle (`classic`² and `classic_dec`²), and a choice axiom,

Check `choice`.
: $\forall (A B : \text{Type}) (R : A \to B \to \text{Prop}), (\forall x : A, \exists y : B, R x y) \to \exists f : A \to B, (\forall x : A, R x (f x)).$

## 4 Subsets and Algebra

This section is dedicated to “algebraic” developments needed or useful for the formalization of finite elements. These complements are mainly based on the `Coquelicot`

<footnote n="14">¹⁴https://github.com/math-comp/multinomials/blob/3f697b7</footnote>
and MathComp libraries (see Section 3), and on our previous Opam package coq-num-analysis.15 The former for its hierarchy of canonical structures representing nested mathematical algebraic structures, and the latter mainly for finite types ("ordinals"), iterated operators ("bigops"), and binomials.

For the formalization of Lebesgue integration for nonnegative measurable functions in [5, 6], we provided features about subsets of a given type, including support for possibly infinite families of subsets. In the present work, we need complements about functions, including support for restrictions of functions to subsets (Section 4.1), and support for finite families of any type (Section 4.2). All related files are now provided in the dedicated Opam package rocq-num-analysis-subset.16

For the formal proof of the Lax-Milgram theorem in [4], we complemented the Coquelicot library with proofs that the type of functions to some algebraic structure is an instance of the same canonical structure, and with minimal support for linear maps. Similarly, in the present work, we need to further add support for finitely iterated operations such as sums and linear combinations, for algebraic substructures (Section 4.3), and for the specific case of finite dimensional linear algebra (Section 4.4). Just as for subsets, all related files are now provided in the dedicated Opam package rocq-num-analysis-algebra.17

Note that Sections 4.3.2, 4.3.4, and 4.4 may be skipped on first reading.

## 4.1 Functions and Restrictions

Functions are basic objects we have to deal with, be they the exact solutions, or piecewise polynomial approximations. This includes the study of special cases such as with input/output sets being empty or not. This also includes handling subsets of functions, typically polynomials, and bijections from/to theses subsets. Unisolvence (see Section 2.4) is precisely a bijection between a vector subspace of functions and a vector space.

File Function therefore provides additional support for functions. First, depending on whether the input and output types are inhabited or not. For instance, given a proof H1 that some type T1 is empty, fun_from_empty H1 provides a witness that the type T1 → T2 is inhabited, and fun_from_empty_is_unit states that this functional type is actually a single-element type (a.k.a. unit type), making the previous witness unique.

Then, based on Rocq.ssrfun (where bijective f is an inductive type taking the proofs that some function is the left and right inverse of f), support for surjective functions is added, and additional support for injective and bijective functions is provided.

Finally, we provide support for the inverse of bijective functions. For instance, the definition of the inverse f_inv takes as argument a proof of bijectivity Hf of some function f and extracts the inverse from bij_EX Hf that states the strong existence of a left-and-right inverse for f (unsurprisingly, the proof of the latter needs ex_EX to reify the inductive weak existential into a strong one in Set).

<footnote n="15">15 https://rocq-prover.org/p/coq-num-analysis/1.0.0</footnote>

<footnote n="16">16 https://rocq-prover.org/p/rocq-num-analysis-subset/2.0.0</footnote>

<footnote n="17">17 https://rocq-prover.org/p/rocq-num-analysis-algebra/2.0.0</footnote>
Function_sub[1] provides support for restrictions of functions to subsets, formalized as predicates (see [5]) that is similar to the support for regular functions. All definitions, such as `injS[1]` and `surjS[1]`, have a name with suffix “S” for “subset”. For instance, given a function `f`: `T1` → `T2` and subsets `P1`: `T1` → `Prop` and `P2`: `T2` → `Prop`, the predicate `funS[1] P1 P2 f := incl (image[1] f P1) P2` states that the image of `P1` by `f` is included in `P2`, and thus that the restriction of `f` to both `P1` (input) and `P2` (output) actually represents a function. Similarly to the predicate `cancel` defined in `Rocq.ssrfun`, given an additional function `g`: `T2` → `T1`, the predicate `canS[1] P1 f g` states that `g` is the left inverse of `f` on `P1`, i.e. “cancels” `f` on `P1`.

Note that the definition of bijectivity for a restricted function is slightly more tricky than for a regular function, as we must ensure that the restrictions of both the function and its inverse actually represent functions, and thus satisfy the predicate `funS`, as shown in the following code snippet.

Definition bijS_spec[1] : (T2 → T1) → Prop :=
    λ g ⇒ funS P1 P2 f ∧ funS P2 P1 g ∧ canS P1 f g ∧ canS P2 g f.
Variant bijS[1] : Prop := BijS_ : ∀ (g : T2 → T1), bijS_spec g → bijS.

See Section 3.1 for more details on the inductive-like `Variant` definition.

The lemma `bijS_ex_uniq_rev[1]` states that a proper restriction, i.e. satisfying the predicate `funS[1]`, and such that all elements in `P2` have a unique preimage in `P1`, also satisfies the `bijS[1]` predicate. Its proof needs the use of `choice` to build the inverse of the restriction. The definition of the inverse of a restriction is similar to that of a regular function, with the same reification of the inductive weak existential into a strong existential in `Set` using `ex_EX[1]`.

## 4.2 Finite Families

We need support for arrays of items of the same type, called here *finite families*, or *n-families* to emphasize their length, be they actual vectors or points in $\mathbb{R}^n$, or simply multi-indices in $\mathbb{N}^n$. Rather than tuples, we choose to use functions from a finite type for their simplicity of use, including direct access, and also because we already accepted to use functional extensionality.

This can seem (and sometimes is) a burden to handle finite families of a given size, with possibly casts between equal size values. Nevertheless, it is related to safety proofs and preventing out-of-array accesses. We ensure the size of each family/vector, even if handling this size means additional proofs. Moreover, we sometimes explicitly need the size of the family, as in Section 8.1.1.

Given a type `E` and some length `n`: $\mathbb{N}$, the type of *n*-families on `E` is $[0..n) \to E$, denoted `E^n`, where the type `ordinal n` (from `fintype` in the `MathComp` library), denoted here `[0..n]`, represents the interval $[0..n) \subset \mathbb{N}$. The type `[0..0]`, denoted `Ø`, being empty, `E^0` is a unit type (`hat0F_unit[1]`). For nonzero length, `[0..n.+1)` is denoted `[0..n]`, and then, `ord0` and `ord_max` respectively represent the bounds 0 and `n`.

Finite_family[1] provides support for structural operations on finite families where most operators are based on transformations on ordinals provided in `ord_compl[1]`. All definitions about finite families have a name containing “F” for “family”, usually as suffix. Let us describe some of the most useful operators. The operator `constF[1]` builds
a constant family, the binary predicate `eqAF` states equality for all indices, and `inclF` states that all items belong to a given subset.

Moreover, we are interested in subfamilies: `liftF_S` removes the first item, `skipF` : `[0..n]` → 'E^n.+1 → 'E^n removes the specified item, it serves as a projector from 'E^n.+1 onto 'E^n, and given a proof `H`: `n1 ≤ n2`, `widenF` `H`: 'E^n2 → 'E^n1 only keeps the `n1` first items from `n2`-families.

Finally, `insertF` : `[0..n]` → `E` → 'E^n → 'E^n.+1 inserts the specified value at the specified location, it serves as an injection from 'E^n to 'E^n.+1. Given a function `f`: `E` → `F`, `mapF` `f`: 'E^n → 'F^n applies `f` to all items of the argument `n`-family on `E`, `to_listF` and `of_listF` provide bijections between types `∀n`, 'E^n and `list E`, and `gather` : '(E → F)^n → 'E → 'F^n transforms a family of functions into the function that gathers all the results.

Now, let us discuss casts. Given a proof `H`: `n1 = n2`, `castF` `H` casts an `n1`-family into an `n2`-family. This may seem silly from a mathematical viewpoint but it has often proved necessary. Many operators on finite families operate through casts. For instance, given a proof `H`: `n = n1 + n2` and an `n`-family `A`, `firstF` `H` (`castF H A`) extracts the first items of `A` into an `n1`-family. And given `A1`: 'E^n1, `A2`: 'E^n2 and a proof `H`: `n1 + n2 = n`, `castF H` (`concatF` `A1 A2`) concatenates the two families into an `n`-family.

*Finite tables*, or (*m*, *n*-tables, are rectangular families of families, i.e. matrices when linear algebra is involved. *Finite table* is similar to *Finite family*, where all definitions have a name containing “T” for “table”, usually as suffix.

## 4.3 Algebraic Structures

For mathematical algebraic structures, the present work relies on the *Coquelicot* [38] hierarchy of canonical structures [42], including the *AbelianMonoid* structure that we recently added to the library. The Figure 3 depicts the part of the hierarchy needed for finite elements, where new features are in blue. Support for functions to algebraic structures, denoted `fct_<Alg_Struct>`, such as `fct_ModuleSpace` is provided in `Hierarchy_compl` and almost all other features have their own file such as `Monoid_sum` and `ModuleSpace_lin_map`. Our inspiration comes from the mathematical books [32, 33].

### 4.3.1 Monoid

We first describe some complements about the *AbelianMonoid* structure. We now have a nonempty structure with a zero, this value being needed for several functions.

For instance, the notation `insertOF` := `λi ⇒ insertFi 0` serves as an injection from 'G^n to 'G^n.+1. We define the new operator on families `itemF` : ∀(n : N), `G` → `[0..n]` → 'G^n that outputs the family made of zeros except a given value in a given position. We also prove `hatOF_eq` that states that 'G^0 is the unit type with the single value 0 (as it is a monoid).

Let us now take advantage of the operator. The introduction of the commutative monoid structure opens the way to define finite sums for both `nat` and `R`, but also
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Hierarchy of algebraic structures used in the present development.
Additions to the Coquelicot library (version 3.4.0) and to previous work [4] are in blue. The main points are: finite iterations of operators using “bigops” from MathComp (`sum`/$\prod$ and `linotherwise`), most morphisms, algebraic substructures, and the affine space structure have been added. Downward arrows specify inheritance, horizontal arrows specify parameters, and the bent double-headed arrows mention that the input structure is shown to be an instance of the output one. For instance, any `ModuleSpace` has the `AffineSpace` structure over itself. Specifics of affine algebra are presented in Section 5.</figcaption></figure>

finite products, and even finite composition of functions (dropping commutativity).
Technically, this is implemented using “bigops” from MathComp that finitely iterate
an operator that is assumed to be associative, with an identity element, and possi-
bly commutative, i.e. endowing the underlying type with a, possibly commutative,
monoid structure. Note that this is not fully automatic since Coquelicot and MathComp
hierarchies are not compatible, although they both represent the same mathemati-
cal objects. A vanilla-Rocq interface is provided through wrappers around the main
results from `bigop`, and the interaction with most operators on families is provided.
For instance, `sum_ind_1` states that the sum of a family is equal to the sum of its first
element and the sum of the rest of the family, it is proved using `big_ord_recl`, and
alike, `sum_skipF`, that extracts any item, is similar to `bigD1`.

In practice, we need both (`N`, `+`) and (`N`, `×`) to be monoids. The Coquelicot structure `AbelianMonoid` is fully abstract, but the identity element is called `zero`, the operator `plus`, and in everybody's mind, this has an additive connotation. Moreover, only one instance of a given canonical structure can be assigned to a given type. Nonetheless, we define aliases `natMul` := `N` and `RMul` := `R` in `MonoidMult`, provide two-way coercions, and endow these “new types” with a multiplicative `AbelianMonoid` structure.
Then, the specific products `prod` and `prod_R` are simply defined as the generic `sum`⁶, and all results are available. Only a few are supplied with a dedicated name.

The iterated composition of functions is similar, with the identity function as the identity element, but the underlying monoid structure is not commutative. Hence, it cannot be represented by an `AbelianMonoid` structure, and `sum` cannot be recycled a second time. Thus, we provide in `MonoidComp`⁷ the definition of `comp_n`⁸ as another “bigop” from `MathComp` and with a similar vanilla-Rocq-like interface. Results about `comp_n` are mainly about the preservation of morphisms (that transports algebraic laws). For instance, `comp_n_mm`⁹ states that the iterated composition of monoid morphisms is a monoid morphism (where “mm” stands for “monoid morphism”). It is also deployed for morphisms of other algebraic structures.

Finally, using the sum (denoted here `Σ`) for natural numbers, support for the concatenation of (not necessarily rectangular) families of families is provided in `ConcatnF`⁹. Given an *n*-family *b* of integers and given an *n*-family of families *A* such that *A* *i* has size *b* *i*, we construct the flattened family of the elements in the expected order, that is of size `Σ b`. This is obtained through the concatenation of lists (from the standard `Rocq` library) and family-to-list transformations (see Section 4.2),

**Definition** `concatnF`⁹ `{E : Type} {n} {b : 'N^n} : (∀ i, 'E^(b i)) → 'E^(Σ b) :=
λ A ⇒ castF _ (of_listF (List.concat (to_listF (λ i ⇒ to_listF (A i))))).`

where the proof that the concatenated family is actually of length `Σ b` is omitted
here for legibility, and simply represented by the unspecified argument `_-`, and the
expressions between curly braces `{ }` represents implicit arguments (see Section 3.1).
Main results about `concatnF` include `concatnF_order`⁹ which states the “continuity”
conditions so that the concatenation of sorted families is also sorted, and `sumnell`⁹
stating that the sum of a concatenation is the double sum of the family of families.

4.3.2 Ring

`Ring.compl`⁹ is for additional definitions and results specific to the `Ring` structure,
having also a one and a multiplication. For instance, the predicate `invertible`⁹ states
that an element of a ring has a left-and-right inverse for the multiplication law, and
the unary operator `inv`⁹ is a total function extracting this inverse when the argument
is indeed invertible, and outputing `0` otherwise, similar to what is done in the `Reals`
with `Rinv`, see Section 3.2.

For barycenters, we use the family operator `part1F`⁹ : `[0..n] → 'K^n → 'K^n.+1` that
inserts one minus the sum at the specified location, such that the output is of
sum 1 (thus creating somehow a partition of unity, hence its name). The notation
`part1F0`⁹ := `part1F ord0` serves as an injection from 'K^n onto the subset of 'K^n.+1 of
elements of sum 1.

Note that a singleton, usually denoted `{0}`, can be endowed with a ring structure
where identity elements for both additive and multiplicative laws are the same. Some
results only hold on specific rings. Thus, we also provide predicates such as `is_field`⁹
for nonzero commutative rings in which all nonzero elements have an inverse, a.k.a.
fields, and `is_field_not_charac_2`⁹ for fields in which 1 + 1 is distinct from 0 (e.g. to
be able to define the middle of two points).
The Kronecker delta function needs (distinct) 0 and 1, hence it is defined in rings, in `Ring_kron`, and some results only hold on nonzero rings. This should be generalized to `semirings` (in which the existence of an additive inverse for all elements is dropped), and thus would also apply to $\mathbb{N}$. But such an algebraic structure is not yet present in Coquelicot's hierarchy.

### 4.3.3 Module Space

`ModuleSpace_compl` collects additional definitions and results about module spaces, i.e. vector spaces with a ring of scalars (not necessarily a field). The scalar multiplication of a vector is called `scal`. Since there is no canonical structure for commutative rings in Coquelicot's hierarchy yet, we choose to state results needing that assumption for vector spaces on the field we use here, real numbers only.

Given a module space $E$ on some ring $K$, given an $n$-family of scalars $a : 'K^n$ and an $n$-family of vectors $u : 'E^n$, `scalF` $a u$ is the $n$-family with items `scal` ($a$ $i$) ($u$ $i$), and its sum is the linear combination `lin comb` $a u$. For instance, results about `lin comb` include `lcNil` stating that the linear combination of the empty family is 0, and `lc_filter_nOF_l` stating that the linear combination remains the same when filtering out zero scalars, or vectors (with suffix `_x`). Note the common prefix `"lc_"`, that stands for "linear combination".

More results on finite-dimensional module spaces are given in Section 4.4.

### 4.3.4 Algebraic Substructures

Given a type $E$ equipped with some algebraic structure, subsets of $E$ that are closed under the laws associated with the structure are of particular interest, such as polynomials being a vector subspace of the vector space of functions. This is declined for each algebraic structure in the dedicated file `<Alg Struct>_sub`. As a general naming rule, `<op>_closed` is a predicate stating closedness under some operator `<op>`, such as `plus_closed` and `sum_closed` in abelian monoids, and `lc_closed` in module spaces. Then, `compatible_
` expresses closedness under all operations constitutive to the algebraic structure corresponding to the shorthand `"
"`: "m", "g", "r", and "ms", respectively for `AbelianMonoid`, `AbelianGroup`, `Ring`, and `ModuleSpace`. Associated statements have `"c
"` in their name (where "c" stands for "compatible"). For instance, `compatible_m` states both `plus_closed` and `zero_closed` (i.e. contains 0), and `cms_lc` states that `compatible_ms` implies `lc_closed`.

Compatible subsets can also be equipped with the same algebraic structure. This is formalized following the representation of subsets by dependent types (see Section 4.1), but this was not needed for the present formalization work.

### 4.3.5 Kernel

A subspace of interest is the kernel of a function. The kernel (or nullspace) of functions can be defined in commutative monoids, as it only needs a zero (with an additive interpretation of it). But of course, most interesting properties about injectivity and bijectivity are only valid in groups, and much more in module spaces. The kernel is the preimage of the zero subspace $\{0\}$. Note that the kernel of the restriction of a function
to some subset is the intersection of the kernel and the subset. We also provide the predicate `KerS0` stating that such kernel of a restriction reduces to the zero subspace, corresponding to the injectivity of the restricted function (when it is a morphism). Given `G1 G2`: `AbelianMonoid`, a subset `PG1`: `G1` → `Prop` and a function `f`: `G1` → `G2`,

**Definition KerS0**: `G1` → `Prop` := preimage `f` (singleton 0).

**Definition KerS**: `G1` → `Prop` := inter `PG1` (`Ker f`).

**Definition KerS0**: `Prop` := ∀ `x1`, `PG1` `x1` → `f x1 = 0` → `x1 = 0`.

For a family of functions `f`: `(G1 → G2)ⁿ`, we provide

**Lemma KerS0_gather Equiv**:

`KerS0 PG1` (`gather f`) ↔ (∀ `x1`: `G1`, `PG1` `x1` → (∀ `i`, `f i x1 = 0` → `x1 = 0`).

where `gather` is defined in Section 4.2, that is used for the unisolvence property (see Section 6.1). When `f` is a morphism and `PG1` is compatible with the algebraic structure, `Ker f` and `KerS PG1 f` are shown to also be compatible with the same algebraic structure. Moreover, we provide the equivalence between `KerS0` and the bijectivity of the restriction of a linear mapping to module subspaces of the same dimension. Given `E1 E2`: `ModuleSpace R_Ring`, `PE1`: `E1` → `Prop` and `PE2`: `E2` → `Prop`, `n1 n2`: `nat` and `f`: `E1` → `E2`,

**Lemma lmS_bijS_KerS0**:

`has_dim PE1 n1` → `has_dim PE2 n2` → `lin_map f` → `funS PE1 PE2 f` → `n2 = n1` → `KerS0 PE1 f` ↔ `bijS PE1 PE2 f`.

where `funS` and `bijS` are defined in Section 4.1 and `has_dim` in Section 4.4 below. It is a consequence of the rank-nullity theorem (see also Section 4.4). Note that this result needs the commutativity of the ring of scalars, and is here only stated for the ring of real numbers.

## 4.4 Finite-dimensional Linear Algebra

Although FEM requires infinite-dimensional spaces for functional spaces, it is also based on a variety of finite-dimensional vector (and affine) subspaces, such as the approximation spaces with their polynomial bases and dual bases, not to mention barycentric coordinates that are at the core of simplicial geometry. Files in folder `Finite_dim` provide support for module subspaces and affine subspaces of finite dimension. Again, definitions and results needing a commutative ring of scalars are stated for the field of real numbers only.

Given a module space `E` on some ring `K`, a subset `PE`: `E` → `Prop` and a finite family of vectors `B`: `E^n`, `lin_span` `B` is the linear span of `B` (i.e. the subset of all linear combinations of `B`), the predicate `lin_gen` `PE` states that `B` spans `PE` (i.e. is a linear generator of `PE`), `lin_indep` `B` that `B` is linear independent, or free, `basis` `PE` that `B` is a basis of `PE` (i.e. is generator and linear independent), and `has_dim` `PE` that `PE` is of dimension `n` (i.e. admits a basis of length `n`). Note that a total function for the dimension with values in `N` would not be satisfactory since all values are meaningful (including 0), and hence there would be no way to discriminate infinite-dimensional subspaces, and subsets that are not subspaces. Thus, we opted for a predicate version, which we preferred to an option type.
We provide technical results from the math book [32]. For instance, linDep_S Lin_Span states that if a family is included in the span of a smaller family, then it is not linearly independent. Then, lin_gen_has_dim states that the linear span of a finite family admits bases, dim_uniq, a.k.a. the dimension theorem, states that all bases have the same cardinality, and rank_nullity_thm states that, for the restriction of any linear map $f: E \to F$ to any finite dimensional subspace $E' \subset E$, one has dim $f(E') + dim ker $f \cap E' = dim E'$. Finally, incomplete_basis_theorem, close to the Steinitz exchange lemma, states that any free family Bf in a subspace PE can be completed from a family generating PE into a basis of PE. These results need commutativity of the multiplication of scalars, they are stated for vector spaces on $\mathbb{R}$.

## 4.5 Binomials

We also need additional results about the type $\mathbb{N}$ of Peano natural numbers. For instance, we provide several variants of double induction schemes, including the following well-founded double induction,

**Lemma nat_ind2_strong:** $\forall (P : \mathbb{N} \to \mathbb{N} \to \text{Prop})$, $(\forall m n, (\forall m1 n1, m1 \le m \to n1 \le n \to (m1, n1) \ne (m, n) \to P m1 n1) \to P m n) \to \forall m n, P m n$.

In Section 8.1.1, we construct multi-indices, that are families whose size is a binomial. We therefore need additional definitions and results about binomial coefficients, based on binomial from MathComp where $\binom{m}{n}$ is denoted $^{\text{C}}(m, n)$. A useful result is `binom_rising_sum_1` stating that

$$ \forall m, n \in \mathbb{N}, \sum_{i=0}^{n} \binom{m+i}{m} = \binom{m+1+n}{m+1}. \qquad (13) $$

Moreover, binomials are at least 1 when the second (or lower) argument is at most equal to the other (e.g. see [46]), and it may be convenient to structurally grant that the value is nonzero. Hence the

**Definition pbinom:** $(m n : \mathbb{N}) : \mathbb{N} := (^{\text{C}}(m + n, m), -1)$.

that returns the natural number $\binom{m+n}{n} - 1$. This is precious in the sequel where the subspace of $d$-variable polynomials of total degree at most $k$ is shown to be of dimension $\binom{d+k}{k}$, i.e. $(\text{pbinom d k})$. +1, thus avoiding technical casts of finite families when a nonzero size is mandatory.

## 5 Affine Algebra

As claimed in Section 2.1, the construction of simplicial Lagrange finite elements is largely a matter of affine objects such as affine mappings and affine subspaces (e.g. hyperplanes), in particular for the geometric aspects (see Section 7). Actually, geometrical transformations are based on the Lagrange polynomials of degree at most 1 (in $\mathbb{P}_1^d$), that happen to be equal to the barycentric coordinates, the key ingredient
in affine algebra, and this explains the developments made for $\mathbb{P}_1^d$, see Sections 7.2.3 and 8.2.2.

In practice, affine spaces are often already equipped with a vector space structure, such as the Euclidean spaces $\mathbb{R}^d$, or more generally the module spaces $K^d$ for any ring $K$. They behave as if the zero vector had lost its special status, and that any other vector (then called *point*) could play a similar role. In this context, affine subspaces are known to be vector subspaces *translated by some vector*, and affine mappings are known to be linear mappings *shifted by a constant*. Thus, affine algebra can reduce to “shifted linear algebra”, and no new skills are needed.

However, the notion of barycenter brings symmetry in formulas, and more importantly, it is the central concept around which all affine aspects revolve. For instance, B. Gostiaux writes in [33, Chap. 1, p. 5]: “une propriété vraiment affine se traduit par des barycentres” [a truly affine property translates into barycenters]. Thus, rather than dealing with the linear part of affine objects, we decided to focus on their true affine nature, and provide the bases of a full-flavored affine algebra library. This vision has the advantage of separating linear aspects from affine ones.

Section 5.1 describes the abstract construction of affine spaces, the definition and first properties of barycenters and affine mappings, and the specific case of finite dimension is discussed in Section 5.2.

## 5.1 Affine Space and Barycenter

An affine space over a module space $V$ (the type of vectors) is any nonempty type $E$ (the type of points) equipped with an external law $\mathrm{vect}^2: E \to E \to V$, denoted with the infix symbol $\dashv$), so that $A \dashv B$ is $\overrightarrow{AB}$. It satisfies the Chasles relation,

**Lemma** $\mathrm{vect\_chasles}^2: \forall (A B C : E), (A \dashv B) + (B \dashv C) = A \dashv C$.

and is such that $\mathrm{vect}^2$ is bijective for all point $O$ (where $O$ stands for a possible “origin” of the affine space). Equivalently, the type $E$ can also be equipped with an external law $\mathrm{transl}^2: E \to V \to E$, denoted with the infix symbol $\dashv\vdash$), that is “associative”,

**Lemma** $\mathrm{transl\_assoc}^2: \forall (A : E) (u v : V), (A \dashv\vdash u) \dashv\vdash v = A \dashv\vdash (u + v)$.

and is such that $\mathrm{transl}^2 O$ is bijective for all point $O$. Either $\mathrm{vect}$ or $\mathrm{transl}$ can be built from the other using the bijectivity property.

The *AffineSpace* structure is another entry point in the hierarchy of algebraic structures. Indeed, the type $E$ of points does not need any prior structure. But it is parameterized by a *ModuleSpace*, which is itself parameterized by a *Ring* (see Figure 3). Any module space $V$ can be equipped with an *AffineSpace* structure on itself using subtraction, $\mathrm{ms\_vect}^2 (A B : V) := B - A$, and then the translation is simply the addition, $\mathrm{ms\_transl\_eq}^2 (A u : V) : A \dashv\vdash u = A + u$ (the $\mathrm{ms\_}$ prefix stands for “module space”).

Symmetrically, any affine space $E$ on some vector space $V$ can be equipped with a *ModuleSpace* structure by choosing any origin $O : E$, defined as the zero. Then, addition of points would be $\lambda (A B : E) \Rightarrow O \dashv\vdash ((O \dashv\vdash A) + (O \dashv\vdash B)) : E$, and the scalar multiplication of a point would be $\lambda (x : K) (A : E) \Rightarrow O \dashv\vdash (\mathrm{scal} x (O \dashv\vdash A)) : E$. Such laws yield a valid canonical structure, and the resulting module space is obviously isomorphic to the initial $V$. But, this generates a cycle in the instantiation rules of the
canonical structures that further confuses the typing system in the presence of module spaces used as affine spaces, and thus it is not provided in the library.

We now describe new operators on finite families specific to affine spaces. Let $V$ be a module space on some ring $K$ and $E$ be an affine space on $V$. Given points $0: E$ and $A: E^n$, $\mathrm{vectF}^0 A$ is the $n$-family of vectors $0 \rightarrow A i$. The law $\mathrm{transl}$ can be lifted similarly into $\mathrm{translF}^0$. The operators $\mathrm{vectF}^0$ and $\mathrm{translF}^0$ are each other inverses, and they both commute with most operations on finite families. Given $B: E^n+1$ and some index $i0$, $\mathrm{frameF}^0 B$ is the $n$-family of vectors $\mathrm{vectF}(B i0)$ ($\mathrm{skipF}^0 i0 B$), i.e. $(\overrightarrow{B_{i0}B_j})_{j \in [0..n]\setminus\{i_0\}}$.

Given scalars $L: K^n$, points $A: E^n$ and $G: E$, the predicate $\mathrm{aff\_comb}^L A G$ stipulates that $\mathrm{lin\_comb}^L G$ ($\mathrm{vectF} G A$) vanishes. When the scalar $\Sigma L$ is invertible (i.e. nonzero in a field, see Section 4.3.2), lemma $\mathrm{ac\_uniq\_r}^L$ (where "ac" stands for "affine combination") states the uniqueness of such point $G$, and $\mathrm{baryc\_EX}^L$ states its "constructive" existence (this uses axiom $\mathrm{ex\_EX}^L$, see Section 3.4). Mathematically, this states

$$ \sum_{i=0}^{n-1} L_i \neq 0 \implies \exists! G \in E, \sum_{i=0}^{n-1} L_i \overrightarrow{GA_i} = 0, $$

or equivalently, that for any $O \in E$, $\sum_{i=0}^{n-1} L_i \overrightarrow{OG} = \sum_{i=0}^{n-1} L_i \overrightarrow{OA_i}$.

Then, $\mathrm{barycenter}^L A$ is a total function extracting that point $G$ when invertible ($\Sigma L$) holds, and outputing some point in the inhabited type $E$ otherwise. Note that when $\Sigma L = 1$ and the real scalars $L_i$'s are nonnegative (and thus lower than 1), the combination is called *convex combination*. The coefficients can be called *weights*, with the analogy between the barycenter and the center of mass of the masses $L_i$'s located at the points $A_i$'s. When $L$ is a constant scalar family, the barycenter is called *isobarycenter*, moreover when there are only two points, it is called the *middle*.

For instance in $\mathbb{R}^d$, the middle $M$ of two points $A$ and $B$ can be computed through $M = A + \frac{B-A}{2} = \frac{A+B}{2}$. The first formula uses the coordinate frame $(A, \overrightarrow{AB})$ with the single coordinate $\frac{1}{2}$. The second formula is the barycentric view, with two coordinates $(\frac{1}{2}, \frac{1}{2})$ summing up to 1, in which $A$ and $B$ play a symmetric role.

Interesting results include the case of Kronecker weights, the possibility to skip zero weights, insensitivity to permutations, and commutativity of barycenters (for a table/matrix of points),

* Lemma $\mathrm{baryc\_l\_kron\_r}^L: \forall \{n\} (A: E^n) (i: [0..n])$, $\mathrm{barycenter}$ (kron $i$) $A = A i$.

* Lemma $\mathrm{baryc\_skip\_zero}^L: \forall \{n\} i0 L (A: E^n+1)$, invertible $(\Sigma L) \rightarrow L i0 = 0 \rightarrow \mathrm{barycenter} L A = \mathrm{barycenter} (\mathrm{skipF} i0 L) (\mathrm{skipF} i0 A)$.

* Lemma $\mathrm{baryc\_permutF}^L$: $\forall \{n\} \{p: [0..n] \rightarrow [0..n]\} \{L\} \{A: E^n\}$, injective $p \rightarrow \mathrm{invertible} (\Sigma L) \rightarrow \mathrm{barycenter} (\mathrm{permutF}^L p L) (\mathrm{permutF} p A) = \mathrm{barycenter} L A$.

* Lemma $\mathrm{baryc2\_rr\_R}^L: \forall \{n1 n2\} L1 L2 (A12: E^n \{n1,n2\})$, $\Sigma L1 \neq 0 \rightarrow \Sigma L2 \neq 0 \rightarrow \mathrm{barycenter} L2 (\mathrm{barycenter} L1 A12) = \mathrm{barycenter} L1 (\mathrm{mapF} (\mathrm{barycenter} L2) A12)$.

The latter needs a commutative ring of scalars (stated here for $\mathbb{R}$); its proof is quite long and cumbersome, it is based on a similar result that is much
easier to prove for linear combinations (lc2 RR R). In module spaces equipped with the affine space structure over themselves, when $\Sigma L = 1$, baryc_ms_eq[✓] states that barycenter L A = lin comb L A. We also provide specialized helpers such as barycenter_ms[✓] that expects points in a module space and not in an affine space. This is due to an incorrect canonical structures instantiation for functions to affine spaces. This is explained in more details in Figure 14 and Section 10.2.

Affine mappings are functions between affine spaces that preserve barycenters ("am" in the names of statements stands for "affine map"). Equivalently, they are functions $f: E_1 \to E_2$ such that $\vec{u}_1 \in V_1 \mapsto f(O_1)f(O_1 + \vec{u}_1) \in V_2$ is a linear map for all $O_1$. In Rocq (see Section 3.1 for more details on the Context keyword),

Context {K : Ring}.

Context {V1 V2 : ModuleSpace K}.

Context {E1 : AffineSpace V1}.

Context {E2 : AffineSpace V2}.

Variable $f : E_1 \to E_2$.

Variable $O1 : E1$.

Definition aff_map[✓] : Prop := $\forall n L (A1 : 'E1^n)$,
invertible ($\Sigma L \to f (\text{barycenter } L A1) = \text{barycenter } L (\text{mapF } f A1)$.

Definition fct_lms[✓] : V1 → V2 := $\lambda u1 \Rightarrow f O1 \dashrightarrow f O1$

Lemma am_lms Equiv : aff_map f ↔ lin_map (fct_lms f O1).

Note that fct_lms is independent of O1 (see am_flm_orig_indep[✓]). Interesting results include comp_am_ms[✓], stating that the composition of affine mappings is an affine mapping, and am_bij_compat[✓], stating that the inverse of a bijective affine mapping is an affine mapping. Moreover, in the case of modules space equipped with the affine space structure, the gathering of a family of affine mappings is an affine mapping,

Context {K : Ring}.

Context {E1 E2 : ModuleSpace K}.

Lemma am_gather_equiv_ms[✓] :
$\forall \{n\} (f : ' (E1 \to E2) ^n), aff_map_ms (gather f) \leftrightarrow (\forall i, aff_map_ms (f i)).$

We also provide specialized helpers such as `aff_map_ms[✓]` that expects a function between module spaces and not affine spaces. Again, this is due to an incorrect canonical structures instantiation for functions to affine spaces (see Figure 14 and Section 10.2 for more details).

Affine subspaces are defined as the subsets that are closed under the barycenter operation (similar to what is done in Section 4.3.4).

Context {K : Ring}.

Context {V : ModuleSpace K}.

Context {E : AffineSpace V}.

Definition compatible_as[✓] (PE : E → Prop) : Prop :=
$\forall n L (A : 'E^n), \text{invertible } (\Sigma L) \to \text{inclF A PE} \to PE \text{ (barycenter } L A).$
In modules spaces equipped with the affine space structure, the kernel of affine mappings is shown to be an affine subspace (Ker\_cas$^2$).

## 5.2 Finite-dimensional Affine Algebra and Barycentric Coordinates

Definitions for finite-dimensional affine spaces are similar to those for vector spaces in Section 4.4. To start with, the affine span of a family of points is the subset made of all barycenters of these points (instead of the linear combinations of a family of vectors for the linear span), and affine independence of a nonempty family of points states that the family of vectors from the first point to all the others ($\overrightarrow{A_0A_i}$) $i \in [1..n]$ is linearly independent. Then, affine generator, affine basis and affine dimension are the same as in linear algebra, except that the cardinality of the (affine) basis is one more than the corresponding module space dimension.

* Context $\{K : \text{Ring}\}$.

* Context $\{V : \text{ModuleSpace K}\}$.

* Context $\{E : \text{AffineSpace V}\}$.

* Variant `aff_span$^2$` $\{n\}$ ($A : 'E^n$) : $E \to \text{Prop} :=$

Aff\_span : $\forall L$, invertible ($\Sigma L \to \text{aff\_span } A$ (barycenter $L A$)).

* Definition `aff_indep$^2$` $\{n\}$ : '$E^n+1 \to \text{Prop} := \lambda A \Rightarrow \text{lin\_indep} (\text{frameF ord0 } A)$.

* Variable $PE : E \to \text{Prop}$.

* Definition `aff_gen$^2$` $\{n\}$ ($A : 'E^n$) : $\text{Prop} := PE = \text{aff\_span } A$.

* Definition `aff_basis$^2$` $\{n\}$ ($A : 'E^n+1$) : $\text{Prop} := \text{aff\_gen } A \wedge \text{aff\_indep } A$.

* Variant `has_aff_dim$^2$` $n : \text{Prop} :=$

Aff\_dim : $\forall (A : 'E^n+1), \text{aff\_basis } A \to \text{has\_aff\_dim } n$.

As expected, `aff_span_cas$^2$` states that affine spans are affine subspaces, and `aff_indep_all Equiv$^2$` states that affine independence can be equivalently expressed with any point, and not just the first. Moreover, `aff_indep_skipF$^2$` states that removing any point keeps the affine independence property, and `aff_indep_permutF$^2$` states that this is the same when permuting the points. Then, `aff_dim_uniq$^2$` states that the affine dimension is unique, and `aff_indep_aff_basis$^2$` states that an affinely independent family is a basis if it is of the correct length. Finally, `has_aff_dim_Rn$^2$` states that $\mathbb{R}^n$ is of affine dimension $n$.

A *simplex* of dimension $d$ is defined by $d+1$ affinely independent points of $\mathbb{R}^d$ (i.e. an affine basis of dimension $d$), that are its vertices. The simplices are tetrahedra in three dimensions, triangles in the plane, line segments on the line, and the point itself in zero dimension.

Consider now a family of points $A : 'E^n$. Then, for any point in their affine span, `aff` states the (constructive) existence of weights of sum 1 such that the point is the barycenter with these weights (with the help of axiom `ex_EX$^2$`, see Section 3.4),

* Lemma `aff`EX$^2$:

$\forall \{n\} \{A : 'E^n\} \{B\}, \text{aff\_span } A B \to \{L \mid \Sigma L = 1 \wedge B = \text{barycenter } L A\}.$
Then, the barycentric coordinates relative to $n$ points $\mathbf{A}$, denoted $\ell^A$, are the $n$ total functions extracting those weights when the point is in the affine span, and zero otherwise. In Rocq,

**Definition** barycCoord $\{n\}$ ($A : 'E^n$) ($B : E$) : $[0..n) \to K :=$
match in_dec $\{aff\_span A\}$ B with
| left HB $\Rightarrow$ proj1_sig ($aff\_span\_EX$ HB)
| right _ $\Rightarrow$ 0
end.

When the family of points is nonempty and affinely independent, the barycentric coordinates are shown to be the unique weights of sum 1 ($bc$ stands for “barycentric coordinates”),

**Lemma** bc_uniq $\{\forall \{n\} \{A : 'E^n.+1\}, aff\_indep A \to \forall \{B L\}, aff\_span A B \to \sum L = 1 \to barycenter L A = B \to L = baryc\_coord A B.$

Now, let $E$ be a module space over a ring $K$ equipped with the affine structure over itself. Let $\mathbf{A}$ be $n+1$ affinely independent points and let $B$ be in the affine span of $\mathbf{A}$. From $bc\_uniq$, we have $B = \sum_{i \in [0..n]} \ell_i^A(B)A_i$, and thus

$$ \overrightarrow{A_0 B} = B - A_0 = \sum_{i \in [1..n]} \ell_i^A(B)A_i + (\ell_0^A(B) - 1)A_0 = \sum_{i \in [1..n]} \ell_i^A(B) \overrightarrow{A_0 A_i}, $$

since $\sum_{i \in [0..n]} \ell_i^A(B) = 1$. This means that $(\ell_i^A(B))_{i \in [1..n]}$ are the $n$ coordinates of $B$ in the frame $(A_0, (\overrightarrow{A_0 A_i})_{i \in [1..n]})$, where $A_0$ plays a specific role. To sum up, with an additional component (such that the sum becomes 1), the barycentric coordinates allow to give all points in $\mathbf{A}$ a symmetric role.

Again, we also provide specialized helpers such as `aff_indep_ms` for module spaces of functions equipped with the affine space structure over themselves. We also give results about the interaction with affine mappings such as `bc_am_ms`, stating that the barycentric coordinates in an affine basis are affine mappings (i.e. $(B \mapsto \ell^A(B))$ is affine), and `am_ext_gen_full_ms`, stating that two affine mappings are equal as soon as they are equal on an affine generator. Moreover, `am_affIRM`, states that the image of an affine span by an injective affine mapping is the affine span of the images, and `am_aff_indep_ms` states that the images of affinely independent points by an injective affine mapping are affinely independent. Finally, `am_has_aff_dim_ms` states that injective affine mappings preserve the affine dimension.

# 6 Formalization of Finite Elements

In this Section we define the Rocq type of a finite element. In Section 2.4, a finite element is defined as a triple $(K, P, \Sigma)$ together with some properties, that are included in a record type. Furthermore, this record is designed to ask minimal effort to build instances without weakening the deductible properties. Note also that an instance of this record is built in Section 9.
In the mathematical literature, the degenerate cases $n_{dof} = 0$ and $d = 0$ are usually prohibited. Indeed, the absence of degree of freedom is of no help in resolving actual PDE problems, and PDE never live in the space $\mathbb{R}^0$ reduced to a single point (which is different from considering a specific equation at some point of a space of nonzero dimension). However, they can be included in the mathematical theory and perfectly fit the formalization. Specifically avoiding them would lead to dispensable hypotheses in lemmas. Therefore, we have chosen to address them. Section 6.1 is first devoted to the general record type, and then, Section 6.2 presents the two degenerate cases ($n_{dof} = 0$ and $d = 0$).

## 6.1 Rocq Definition as a Record

Let us start by presenting some definitions useful for the geometric element $K$. As specified in Section 2.4, $K$ is a (nondegenerate) polytope. In this article, we focus on the simplex shape only, but anticipating future works, we also provide the possibility of using cuboids (there is a field called **shape** in the record type below),

**Inductive shape_type**<sup>□</sup> := Simplex | Cuboid.

The number of vertices (field `nvtx`) is given by

**Definition nvtx_of_shape**<sup>□</sup> (d : $\mathbb{N}$) (shp : shape_type) : $\mathbb{N}$ :=
match shp with Simplex ⇒ d.+1 | Cuboid ⇒ 2^d end.

Indeed, a simplicial polytope in dimension $d$ has $d+1$ vertices, such as a triangle for $d=2$ or a tetrahedra for $d=3$.

The vertices (field `KVertices`) are a family of points in $\mathbb{R}^d$ of the appropriate size `nvtx_ofShape`. Then, $K$ is specified as their convex hull (field `KGeom`), i.e. the subset of all convex combinations (see Section 5.1),

**Variant convex_hull**<sup>□</sup> {E : ModuleSpace R_Ring} {$n$} (vtx : 'E^n) : E → Prop :=
| Cvx : ∀ L, (∀ i, 0 ≤ L i) → Σ L = 1 → convex_hull vtx (barycenter L vtx).

Note that the set $K$ is constructed from the vertices, and thus the geometry is entirely defined by the vertices. The nondegeneration of the geometrical element is not explicitly specified. Instead, we rely on the affine independence of the vertices which is assumed in most results, and is sufficient in the simplicial case, see Section 10.3.

Considering the case of unknown functions from $\mathbb{R}^d$ to $\mathbb{R}$ only ($q = 1$, see Section 2.4), the approximation space $P$ is a finite-dimensional vector subspace of $\mathcal{F}(\mathbb{R}^d, \mathbb{R})$, the space of functions from $\mathbb{R}^d$ to $\mathbb{R}$ represented in Rocq as $FRd^d$ (fields `PApprox` for the subset and `PApprox_has_dim` for its dimension). The degrees of freedom $\Sigma$ are a collection of $n_{dof}$ linear forms in $\mathcal{F}(\mathbb{R}^d, \mathbb{R}) \to \mathbb{R}$ (fields `ndof`, `S_dof`, and `S_dof_lm`).

Then, the record type `FE` is designed to ask minimal effort to build instances.

**Record FE**<sup>□</sup> (d : $\mathbb{N}$) : Type := mk FE {
  shape : shape_type;
  nvtx : $\mathbb{N}$ := nvtx_ofShape d shape;
  KVertices : [0..nvtx) → 'R^d;
  K_geometry : 'R^d → Prop := convex_hull KVertices; (* geometric element K *)
  ndof : $\mathbb{N}$;
Pricc : FRd d → Prop; (* approximation space P *)
Pricc has dim : has dim Pricc ndof;
Sdof : [0..ndof) → FRd d → R; (* degrees of freedom Σ *)
Sdof lm : ∀ i, lin map (Sdof i);
unisolvence inj : KerS0 Papprox (gather Sdof);
}.

where mk FE is the type constructor of FE.

The fields ndof, PApprox, and Sdof specify the data parts of P and Σ. The
proof parts reduce to three fields: first, PApprox has dim where the dimension of
PApprox is equal to the size of Σ, then Sdof lm for the linearity of degrees
of freedom, and last unisolvence_inj for the sole injectivity of ΦΣ, expressed as
KerS0 PApprox (gather Sdof).

However, the unisolvence property specified in Section 2.4 is the bijectivity of ΦΣ
from P onto R<sup>n<sub>dof</sub></sup>. This bijectivity property implies dim P = n<sub>dof</sub>. Conversely, from
lemma lmS_bijS_KerS0 (see Section 4.3.5), if dim P = n<sub>dof</sub>, then ΦΣ is bijective iff
it is injective, i.e. iff its kernel restricted to P is reduced to {0}, which is exactly the
field unisolvence_inj. Thus, to satisfy the unisolvence property, we provide the lemma
giving bijectivity back, whose proof is merely an application of lmS_bijS_KerS0 and of
the property fields of the record,

Lemma unisolvence :

`∀ \{d\} (fe : FE d), bijS (Pricc fe) fullset (gather (Sdof fe)).`

The following extensionality result states that it is sufficient to check equality of
data fields to ensure the equality of two instances of the record type FE.

Lemma FE_ext :

$$
\begin{align*}
\forall \{d\} \{fe1 fe2 : FE d\} (Hs : \text{shape fe1} = \text{shape fe2}) (Hn : \text{ndof fe1} = \text{ndof fe2}), \\
K\_vertices fe1 = \text{castF } (K\_vertices fe2) \rightarrow P\_approx fe1 = P\_approx fe2 \rightarrow \\
S\_dof fe1 = \text{castF } (S\_dof fe2) \rightarrow fe1 = fe2.
\end{align*}
$$

6.2 Degenerate Cases

We now focus on the two degenerate cases, and build the corresponding instances.

When $n_{\text{dof}} = 0$, the approximation space $P$ is of dimension $0$, i.e. it reduces to the
sole zero function (see `has_dim_0_rev`), and the empty family of degrees of freedom $\Sigma$
is zero (see `hatOF_eq` in Section 4.3.1). Of course, all degrees of freedom are linear,
since there are none. Unisolvence is also straightforward since the only approximation
function is the zero function.

Moreover, lemma FE_ndof_0_uniq states that given any geometric element K
.shape and vertices), there is only one finite element on K with no degree of freedom.

When $d = 0$, there is only one spatial point ($\mathbb{R}^0 = \{0\}$), and only one geometric element $K$, equal to the whole space ($\text{hatm0T_eq}$).$^{18}$ The functional space $\mathcal{F}(\mathbb{R}^0, \mathbb{R})$ is

<footnote n="18"><sup>18</sup>This really is “the realm of Pointland”, from Edwin A. Abbott’s *Flatland* (1884).</footnote>
the one-dimensional vector space of constant functions (FR0_has_dim[1]), since there is only one input value; it is generated by the constant function 1 (FR0_bin_gen[1]). And if not of dimension 0 (i.e. not the previous case above), the approximation space P is also of dimension 1 (FR0_sub_dim[1]), i.e. equal to the whole space F(R⁰, R) (FR0_sub_eq[1]). The dual space F(R⁰, R) → R is generated by the linear form (f ↦→ f(0)) evaluating constant functions at the single spatial point (FR0_dual_bin_gen[1]), and two linear forms are equal if they coincide on the constant function 1 (FR0_dual_ext[1]). Then, the unisolvence property is equivalent to the noncancellation of the degree of freedom on 1 (FE_d_0_S_dof_n0[1] and unisolvence_inj_d_0[1]).

We can now construct a record for $d = 0$ and any shape shp,

**Definition FE_d_0[1]: FE 0 :=**
$$mk FE shp 0 1 P \_approx \_has \_dim \_d \_0 \_S \_dof \_lm \_d \_0 \_unisolvence \_inj \_d \_0.$$

Again, lemma FE_d_0_uniq[1] states that given any shape (which is meaningless here), there is only one finite element on the single zero-dimensional geometric element, up to a scaling factor, that is f(0), on the nonzero single degree of freedom.

# 7 Geometric Transformations

This section addresses geometrical issues such as how to map a geometrical element onto another, and how to control the mapping of the hyperplanes carrying the hyperfaces of a geometrical element.

One of the main ingredients of the finite element method is to bring computations on mesh cells of a wide range of aspect ratios back to the sole case of the unit rectangular reference element (see Section 2.3.3). The function mapping the reference element to any given cell of the mesh is called the *geometric transformation*; it is supposed to be a smooth invertible mapping, such as a family of polynomial functions with invertible Jacobian. The action of such geometric transformation can be seen as a change in the local representation of the spatial dimensions. A salient aspect of simplicial geometries is that geometric transformations are actually affine mappings, and that they can be considered globally, on the whole space $R^d$.

Section 7.1 describes the reference element. The geometric transformation mapping the reference element onto a generic one is built in Section 7.2, it is based on Lagrange polynomials of degree at most 1. Then, Section 7.3 is about face hyperplanes and how they are mapped by the geometric transformation. Finally, Section 7.4 focuses on the geometric transformation of a full finite element.

## 7.1 Reference Geometry

The reference element associated with the simplicial Lagrange finite elements is the unit rectangular simplex K. Its extremities are called the *reference vertices*, denoted $\hat{v}$. Quantities related to the reference element are all referred to with a “hat” notation.

Let $d$ be the spatial dimension. The $d+1$ reference vertices $\hat{v}$ in $R^d$ are defined by

$$\hat{v}_0 \stackrel{\text{def.}}{=} 0, \quad \text{and} \quad \forall i \in [1..d], \quad \hat{v}_i \stackrel{\text{def.}}{=} \delta_{i-1},$$
where $\mathbf{0}$ is the zero family, and $(\delta_i)_{i \in [0..d-1]}$ is the Kronecker family, such that $\delta_{ij} = 0$ if $i \neq j$ and $\delta_{ii} = 1$, here of length $d$. It is straightforwardly formalized as

Context $\{d : \mathbf{nat}\}$.

**Definition** $vtx\_ref^\square (i : [0..d]) : '\mathbb{R}^d :=$
match ord\_eq\_dec^\square i ord0 with
| left _ $\Rightarrow$ 0
| right Hi $\Rightarrow$ kron (lower\_S Hi)
end.

where $\mathbf{kron}^\square$ is the Kronecker delta function, and $\mathbf{lower\_S}^\square$ is the decrement of an ordinal from the proof it is nonzero.

The reference vertices are shown to be affinely independent (see Section 5.2),

**Lemma** $vtx\_ref\_aff\_indep^\square : aff\_indep\_ms vtx\_ref$.

They also span the affine space $\mathbb{R}^d$, and thus constitute an affine basis of $\mathbb{R}^d$. As a consequence (see Section 5.2), any point $\hat{x} \in \mathbb{R}^d$ is a barycenter of the reference vertices $\hat{v}$, and the weights are its reference barycentric coordinates $\hat{\ell}(\hat{x}) \stackrel{\text{def.}}{=} \ell^{\hat{v}}(\hat{x})$,

**Lemma** $bc\_vtx\_ref\_uniq^\square : \forall (x\_ref : '\mathbb{R}^d) (L : '\mathbb{R}^d.+1), \Sigma L = 1 \rightarrow$
$barycenter\_ms L vtx\_ref = x\_ref \rightarrow L = baryc\_coord\_ms vtx\_ref x\_ref$.

The equality of affine mappings can be checked on an affine basis only, thus we have the following extensionality result

**Lemma** $am\_ext\_vtx\_ref^\square : \forall \{f g : '\mathbb{R}^d \rightarrow '\mathbb{R}^d\}, aff\_map\_ms f \rightarrow aff\_map\_ms g \rightarrow$
$(\forall i, f (vtx\_ref i) = g (vtx\_ref i)) \rightarrow f = g$.

## 7.2 Lagrange Polynomials and Geometric Transformation

Lagrange polynomials are interpolation polynomials. The expression of Lagrange polynomials is simple when $d = 1$, or when $k = 1$. The former case is treated in Section 8.2.2. In the affine case (degree at most 1), Lagrange polynomials have a very simple expression when the interpolation points are the reference vertices (Section 7.2.1). In the general case, the affine transformation that maps the reference vertices to some $d+1$ (generic) vertices (Section 7.2.2) is the convenient way to define the (generic) Lagrange polynomials (Section 7.2.3).

The affine Lagrange polynomials are actually shown to be the barycentric coordinates with respect to the chosen generic vertices, and the geometric transformation is at the core of the construction and use of the simplicial Lagrange finite elements.

### 7.2.1 $\hat{\mathcal{L}}$ Reference Lagrange Polynomials

Given some dimension $d$, the reference Lagrange polynomials $\hat{\mathcal{L}}: \mathbb{R}^d \rightarrow \mathbb{R}^{d+1}$ of degree at most 1 are defined for all $\hat{x} \in \mathbb{R}^d$ as

$$ \hat{\mathcal{L}}_0(\hat{x}) \stackrel{\text{def.}}{=} 1 - \sum_{j \in [0..d)} \hat{x}_j, \quad \text{and} \quad \forall i \in [1..d], \hat{\mathcal{L}}_i(\hat{x}) \stackrel{\text{def.}}{=} \hat{x}_{i-1}. $$
This family of functions is formalized as

Context {d : nat}.

Definition LagPd1_ref[?](i : [0..d]) : FRd d :=
match ord_eq_dec i ord0 with
| left _ ⇒ λ x_ref ⇒ 1 − Σ x_ref
| right Hi ⇒ λ x_ref ⇒ x_ref (lower_S Hi)
end.

The reference Lagrange polynomials are shown to be of sum 1 (LagPd1_ref_sum[?]),
and to satisfy

$$
\forall \hat{\mathbf{x}} \in \mathbb{R}^d, \quad \sum_{i \in [0..d]} \hat{\mathcal{L}}_i(\hat{\mathbf{x}}) \hat{v}_i = \hat{\mathbf{x}} \qquad (14)
$$

(LagPd1_ref_barycF[?], scalars adding up to 1, this linear combination is actually a barycenter). Thus, from bc_vtx_ref_uniq[?](see Section 7.1), the reference Lagrange polynomials are indeed the barycentric coordinates in the system of reference vertices, $\hat{\mathcal{L}} = \hat{\ell}$ (see LagPd1_ref_bcf[?]). As such, the reference Lagrange polynomials are affine mappings (see Section 5.2, LagPd1_ref_am[?]), they do interpolate at the reference vertices, $\hat{\mathcal{L}}_j(\hat{v}_i) = \delta_{ij}$ (see LagPd1_ref_kron_vtx[?]), and $\hat{\mathcal{L}}$ is injective, and surjective onto the hyperplane of $\mathbb{R}^{d+1}$ of equation $\sum_{i \in [0..d]} x_i = 1$,

Lemma LagPd1_ref_inj[?]: injective (gather LagPd1_ref).

Lemma LagPd1_ref_sujS[?]: surjS fullset (λ x ⇒ Σ x = 1) (gather LagPd1_ref).

In Section 8.2.2, they are also shown to form a basis of the space $\mathbb{P}_1^d$ of $d$-variate affine polynomials.

7.2.2 Affine Geometric Transformation

<figure><img src="image_5.png" /><figcaption>Fig. 4 Affine geometric transformation $T^v$ in the case $d = 2$.
Given $\mathbf{v} = (\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2)$, $T^v$ maps the unit rectangular triangle $\hat{K}$ to the triangle $K$ of vertices $\mathbf{v}$.</figcaption></figure>

Given $d+1$ (generic) vertices $\mathbf{v}$ of a (generic) geometric element $K \subset \mathbb{R}^d$, an affine mapping sending reference vertices $\hat{\mathbf{v}}$ onto $\mathbf{v}$ preserves barycenters, that is barycentric coordinates (see the definition of aff_map in Section 5.1). Thus, the affine geometric transformation $\hat{T}^v$ outputs the point with the same barycentric coordinates as the
input point in the system $\hat{\mathbf{v}}$, but in the system of vertices $\mathbf{v}$ (see Figure 4 for an example in the case $d = 2$). It is defined by (see (7) in Section 2.3.3)

$$T^{\mathbf{v}} : \hat{\mathbf{x}} \in \mathbb{R}^d \longmapsto \sum_{i \in [0..d]} \hat{\mathcal{L}}_i(\hat{\mathbf{x}}) \mathbf{v}_i \in \mathbb{R}^d,$$

and in Rocq (emphasizing the barycenter, since the weights add up to 1),

Context $\{d : \mathbf{nat}\}$.

Definition $\mathbf{T}_{\text{geom}}^{\mathcal{Z}}$ (vtx : $[0..d] \to \mathbb{R}^{\mathcal{d}}$) : $\mathbb{R}^{\mathcal{d}} \to \mathbb{R}^{\mathcal{d}} := \lambda \ x_{\text{ref}} \Rightarrow \text{barycenter}_{\text{ms}}$ ($\lambda$ i $\Rightarrow$ LagPd1\_ref i $x_{\text{ref}}$) vtx.

The geometric transformation $T^{\mathbf{v}}$ is an affine mapping ($\mathbf{T}_{\text{geom}}^{\mathcal{Z}}$) since

$$T^{\mathbf{v}} - T^{\mathbf{v}} \mathbf{0} = \left( \hat{\mathbf{x}} \longmapsto \sum_{i \in [0..d]} (\hat{\mathcal{L}}_i(\hat{\mathbf{x}}) - \hat{\mathcal{L}}_i(\mathbf{0})) \mathbf{v}_i \right)$$

is a linear map (because the $\hat{\mathcal{L}}_i$'s are affine mappings). We have $T^{\mathbf{v}} \hat{\mathbf{v}} = \mathbf{v}$ ($\mathbf{T}_{\text{geom}}_{\text{vtxF}}^{\mathcal{Z}}$), which actually means for all $i \in [0..d]$, $T^{\mathbf{v}} \hat{\mathbf{v}}_i = \mathbf{v}_i$, and this provides the formula for the composition of geometric transformations,

$$\forall \mathbf{v}, \mathbf{w} \in (\mathbb{R}^d)^{d+1}, \quad T^{\mathbf{v}} \circ T^{\mathbf{w}} = T^{\mathbf{v}} \mathbf{w}, \qquad (15)$$

that is formalized as

Lemma $\mathbf{T}_{\text{geom\_comp}}^{\mathcal{Z}}$ : $\forall$ (vtx1 vtx2 : $[0..d] \to \mathbb{R}^{\mathcal{d}}$), $\mathbf{T}_{\text{geom}}$ vtx1 $\circ$ $\mathbf{T}_{\text{geom}}$ vtx2 = $\mathbf{T}_{\text{geom}}$ (mapF ($\mathbf{T}_{\text{geom}}$ vtx1) vtx2).

Indeed, from $\mathbf{am\_ext\_vtx\_ref}^{\mathcal{Z}}$, it suffices to check that both functions agree on $\hat{\mathbf{v}}$, which is the case.

Moreover, $T^{\mathbf{v}}$ is the identity map ($\mathbf{T}_{\text{geom\_ref}}^{\mathcal{Z}}$). Introducing $(s \cdot) : \mathbb{R}^d \to \mathbb{R}^d$ the scaling by $s$, defined for all $i \in [0..d)$ by $((s \cdot) \mathbf{x})_i \stackrel{\text{def.}}{=} s \cdot x_i$, we have ($\mathbf{T}_{\text{geom\_scal\_ref}}^{\mathcal{Z}}$)

$$\forall s \in \mathbb{R}^*, \quad T^{(s \cdot)} \hat{\mathbf{v}} = (s \cdot). \qquad (16)$$

And the geometric elements $\hat{K}$ and $K$ being the convex hull of their vertices, we also have $T^{\mathbf{v}}(\hat{K}) = K$ ($\mathbf{T}_{\text{geom\_K\_geom\_eq}}^{\mathcal{Z}}$), see Figure 4 for an example when $d = 2$.

Assume now that the vertices $\mathbf{v}$ are affinely independent. Then, the barycentric coordinates in the system $\mathbf{v}$ are unique, and injectivity and surjectivity of $T^{\mathbf{v}}$ follow that of the reference Lagrange polynomials $\hat{\mathcal{L}}$ (see Section 7.2.1). Thus, $T^{\mathbf{v}}$ is bijective onto the whole $\mathbb{R}^d$, and we can build its inverse $(T^{\mathbf{v}})^{-1}$ (see Section 4.1). In Rocq,

Context $\{vtx : [0..d] \to \mathbb{R}^{\mathcal{d}}\}$.
Hypothesis $Hvtx : \text{aff\_indep\_ms}$ vtx.

Lemma $\mathbf{T}_{\text{geom\_bij}}^{\mathcal{Z}}$ : bijective ($\mathbf{T}_{\text{geom}}$ vtx).
Definition $\mathbf{T}_{\text{geom\_inv}}^{\mathcal{Z}}$ : $\mathbb{R}^{\mathcal{d}} \to \mathbb{R}^{\mathcal{d}} := f_{\text{inv}} \mathbf{T}_{\text{geom\_bij}}$.
Note that $T^v$ is also proved to be a bijection from $\hat{K}$ onto $K$ ($T_{\text{geom\_bijS\_K\_geom}}$). The inverse geometric transformation $(T^v)^{-1}$ is also an affine map $(T_{\text{geom\_inv\_am}})$, and we have $(T^v)^{-1} \mathbf{v} = \hat{\mathbf{v}}$ (see $T_{\text{geom\_inv\_vtx}}$) and $(T^v)^{-1}(K) = \hat{K}$ (see $T_{\text{geom\_inv\_K\_geom\_eq}}$).

As explained in Section 3, we have made implicit some variables. Here, the function $T_{\text{geom}}$ requires only the vertices (and not the dimension), while $T_{\text{geom\_inv}}$ requires only the proof of affine independence of the vertices (and neither the vertices, nor the dimension). It reduces the number of parameters to be given, hence it improves the readability.

### 7.2.3 $\mathcal{L}^v$ Lagrange Polynomials

The geometric transformation, or rather its inverse, is a convenient way to bring computations back to the simpler case of the reference simplex. Given $d+1$ affinely independent vertices $\mathbf{v}$, the associated affine Lagrange polynomials $\mathcal{L}^v: \mathbb{R}^d \to \mathbb{R}^{d+1}$ are defined as

$$ \forall i \in [0..d], \forall \mathbf{x} \in \mathbb{R}^d, \quad \mathcal{L}_i^v(\mathbf{x}) \stackrel{\text{def.}}{=} \hat{\mathcal{L}}_i \circ (T^v)^{-1}(\mathbf{x}). $$

This is formalized as

Context $\{d: \text{nat}\}$.

Context $\{\text{vtx}: [0..d] \to \mathbb{R}^d\}$.

Hypothesis $\text{Hvtx}: \text{aff\_indep\_ms vtx}$.

Definition $\text{LagPd1}^v$ ($i: [0..d]$): $\text{FRd} : = \lambda \mathbf{x} \Rightarrow \text{LagPd1\_ref i} (T_{\text{geom\_inv}} \text{Hvtx} \mathbf{x})$.

Unsurprisingly, all affine properties of the reference affine Lagrange polynomials $\hat{\mathcal{L}}$ are transported to $\mathcal{L}^v$ by the affine geometric transformation $(T^v)^{-1}$. Their sum is 1 ($\text{LagPd1\_sum}^v$). Since they satisfy the property ($\text{LagPd1\_barycF}^v$)

$$ \forall \mathbf{x} \in \mathbb{R}^d, \quad \sum_{i \in [0..d]} \mathcal{L}_i^v(\mathbf{x})\mathbf{v}_i = \mathbf{x}, $$

they are also the barycentric coordinates in the system of vertices $\mathbf{v}$, i.e. $\mathcal{L}^v = \mathcal{L}^v$ (see $\text{LagPd1\_bcF}^v$), they are affine mappings ($\text{LagPd1\_am}^v$), and they interpolate at the vertices, $\mathcal{L}_j^v(\mathbf{v}_i) = \delta_{ij}$ (see $\text{LagPd1\_kron\_vtx}^v$).

Another property is that $\mathcal{L}^v$ forms a basis of the space $\mathbb{P}_1^d$ of $d$-variate affine polynomials (see Section 8.2.2).

## 7.3 Geometric Transformations and Face Hyperplanes

Faces play an important role in the numerical resolution of PDE as they carry the boundary conditions of the problem to solve, and also the continuity conditions between elements. In the present work, we deal with affine mappings that can operate on the whole space, they are total functions. This is why we rather consider the
(hyper)planes carrying the faces of the geometric element (see Figure 6). In Section 9.3, both steps of the proof of unisolvence involve face hyperplanes.

In geometry, faces in any dimension are rather called hyperfaces. A hyperplane is an affine subspace with one dimension less than its ambient space (e.g. a plane in 3D, a line in 2D, or a point in 1D). Given a hyperface (boundary) of a simplex, the hyperplane containing that hyperface is called face hyperplane. This is the generalization to any dimension of what is a plane containing one of the triangular face of a tetrahedron. Affine mappings transport affine subspaces, and when they are injective, the affine dimension is preserved (see Section 5.2), thus the image of a (face) hyperplane is a (face) hyperplane.

Reference and generic face hyperplanes are first discussed in Section 7.3.1, then we study two variants of the geometric transformation based on the way they deal with face hyperplanes in Sections 7.3.2 and 7.3.3.

### 7.3.1 Face Hyperplanes

Face hyperplanes are referred to by the vertex opposite them. The reference face hyperplane $\hat{\mathcal{H}}_i$ opposite the reference vertex $\hat{v}_i$ is the hyperplane generated by the other reference vertices $(\hat{v}_j)_{j \neq i}$ (see Figure 6 on the left). It is formalized as

**Definition Hface_ref** (i : [0..d]) : $\mathbb{R}^d \to \text{Prop} := \text{aff\_span\_ms} (\text{skipF i vtx\_ref}).$

An important property is that the reference face hyperplanes are the kernels of the reference affine Lagrange polynomials, $\hat{\mathcal{H}}_i = \ker \mathcal{L}_i$, in $\text{Rocq}$,

**Lemma Hface_ref_is_Ker** : $\forall i$, $\text{Hface\_ref i} = \text{Ker} (\text{LagPd1\_ref i})$.

Since the kernel of an affine mapping is an affine subspace and $\mathcal{L}_i(\hat{v}_j) = 0$ for $j \neq i$, it is sufficient to show the second inclusion $\ker \mathcal{L}_i \subset \hat{\mathcal{H}}_i$. Let $\hat{x}$ such that $\mathcal{L}_i(\hat{x}) = 0$. Then, from (14), we can skip the zero weight in the sum/barycenter, and $\hat{x}$ is a barycenter of the sole other reference vertices $(\hat{v}_j)_{j \neq i}$. It belongs to their span $\hat{\mathcal{H}}_i$. $\square$

A noticeable property is that the hyperfaces of a simplex are in affine bijection with simplices in one dimension less. Extended to hyperplanes, all face hyperplanes are in bijection with the whole space with one dimension less. Given $d > 0$, let $[\hat{i}]$ be the injection from $\mathbb{R}^{d-1}$ to $\mathbb{R}^d$, whose range is $\hat{\mathcal{H}}_i$, and that is defined for all $\hat{x}' \in \mathbb{R}^{d-1}$ by

$$ \hat{I}_0(\hat{x}') \stackrel{\text{def.}}{=} \left( 1 - \sum_{j \in [0..d-1]} \hat{x}', \hat{x}' \right), \qquad (17) $$

$$ \forall i \in [1..d], \quad \hat{I}_i(\hat{x}') \stackrel{\text{def.}}{=} (\hat{x}'_0, \dots, \hat{x}'_{i-2}, 0, \hat{x}'_{i-1}, \dots, \hat{x}'_{d-2}). \qquad (18) $$

Then, $\hat{I}_i$ is a bijective affine mapping from $\mathbb{R}^{d-1}$ onto $\hat{\mathcal{H}}_i$. This is formalized as

**Context** $\{d1 : \text{nat}\}$.

Let $d := d1+1$.

**Definition inj_Hface_ref** (i : [0..d]) : $\mathbb{R}^d \to \mathbb{R}^d := \text{match ord\_eq\_dec i ord0 with}$
| left _ ⇒ part1F0

| right Hi ⇒ insertOF (lower_S Hi)

end.

Lemma `inj_Hface_ref_am_2`: $\forall i$, $aff_map_ms$ (`inj_Hface_ref` $i$).

Lemma `inj_Hface_ref_bijS_2`: $\forall i$, `bijS` `fullset` (`Hface_ref` $i$) (`inj_Hface_ref` $i$).

where the `Let` command is for local definitions (see Section 3.1). Note that the positivity of $d$ is granted structurally to avoid building a function with a type '$\mathbb{R}^d$. $-1 \to '$\mathbb{R}^d$' depending on the proof that $d > 0$. Note also that `$lower_S Hi$` is indeed of type `[0..d]`. The operators `$part1F0_2$` and `$insertOF_2$` ensure that the output is in the correct kernel, or reference face hyperplane. Indeed, the former is the insertion in front of the complement to 1 (see Section 4.3.2), thus its output is of sum 1, and the latter is the insertion of 0 in the specified location (see Section 4.3.1), thus its output is in $\ker \hat{\mathcal{L}}_i$.

The face hyperplanes of generic simplices are obtained through the affine geometric transformation. Given $d + 1$ vertices $\mathbf{v}$, the face hyperplane opposite vertex $\mathbf{v}_i$ for $i \in [0..d]$ is defined as $\mathcal{H}_i^v \stackrel{\text{def.}}{=} T^v(\hat{\mathcal{H}}_i)$. And when the vertices are affinely independent, $T^v$ is bijective, and we also have

$$ \ker \mathcal{L}_i^v = \ker (\hat{\mathcal{L}}_i \circ (T^v)^{-1}) = ((T^v)^{-1})^{-1} (\ker \hat{\mathcal{L}}_i) = T^v(\hat{\mathcal{H}}_i) = \mathcal{H}_i^v. $$

Moreover, the restriction of $T^v$ to $\hat{\mathcal{H}}_i$ is a bijection onto $\mathcal{H}_i^v$. This is formalized as

**Context** $\{d : \text{nat}\}$.

**Definition** `Hface_2` (`vtx`: `[0..d]` $\to$ `$\mathbb{R}^d$`)(`i`: `[0..d]`) : `$\mathbb{R}^d$` $\to$ `Prop` := image (`T_geometry` `vtx`) (`Hface_ref` $i$).

**Context** `$vtx$: `[0..d]` $\to$ `$\mathbb{R}^d$`.

**Hypothesis** `Hvtx`: `aff_indep_ms` `vtx`.

Lemma `Hface_is_Ker_2`: $\forall i$, `Hface` `vtx` $i$ = $\text{Ker}$ (`LagPd1` `Hvtx` $i$).

Lemma `$T_geometry_bijS_Hface_2$`: $\forall i$, `$bijS$` (`Hface_ref` $i$) (`Hface` `vtx` $i$) (`T_geometry` `vtx`).

The hyperplane $\mathcal{H}_i^v$ is also the affine span of the vertices ($\mathbf{v}_j$)$_{j \neq i}$ (`Hface_aff_span_2`).

### 7.3.2 Injection onto a Face Hyperplane

The first step of the inductive proof of unisolvence in Section 9.3 deals with the increment of the dimension for a generic simplex: a result on the reference rectangle unit simplex in one dimension less is transported to the first face hyperplane $\mathcal{H}_0^v$ opposite the first vertex $\mathbf{v}_0$. It is thus necessary to transport information on some simplex in dimension $d-1$ to a face hyperplane $\mathcal{H}_i^v$ in dimension $d$ (with $d > 0$ and $i \in [0..d]$). This is simply obtained by using the composition $[I_i^v]^{\text{def.}} = T^v \circ \hat{\mathcal{I}}_i$, see an illustration for $d = 3$ in Figure 5. In Rocq,

**Context** `$\{d1 : \text{nat}\}$.

**Let** $d := d1+1$.
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Injection onto a face hyperplane $I_0^V$ in the case $d = k = 3$.</figcaption></figure>

The reference triangle $\hat{K}^{d-1}$ is first mapped onto the 0-th face of $\hat{K}$, opposite vertex $\hat{v}_0$, and then to the 0-th face of $K$, that is contained in the face hyperplane $\mathcal{H}_0^V$ supported by $(v_1, v_2, v_3)$. The hyperplanes are depicted in light blue and the faces in a darker blue.

After reading Section 9.1 about nodes, note that the correspondence between the reference nodes in the triangle ($d-1=2$) and the face nodes of the tetrahedron ($d=3$) is illustrated by the colors. We have $I_0^V(\hat{a}_{(i,j)}^{d-1}) = a_{(3-(i+j),i,j)}^V$, for all $(i,j) \in \mathcal{A}_2^3$.

$$ \textbf{Variable } vtx : [0..d] \to \mathbb{R}^d. $$

$$ \textbf{Definition } \text{inj\_Hface}^\square (i : [0..d]) : \mathbb{R}^d1 \to \mathbb{R}^d := T\_geom vtx \circ \text{inj\_Hface\_ref } i. $$

where the positivity of $d$ is granted structurally because of $\text{inj\_Hface\_ref}^\square$.

As a composition of affine mappings, $I_i^V$ is an affine mapping ($\text{inj\_Hface\_am}^\square$). Moreover, as expected, when the vertices $v$ are affinely independent, $I_i^V$ is bijective from $\mathbb{R}^{d-1}$ onto $\mathcal{H}_i^V$, as the composition of bijective functions. In Rocq,

$$ \textbf{Lemma } \text{inj\_Hface\_bijS}^\square : \forall i, \text{bijS fullset } (\text{Hface vtx } i) (\text{inj\_Hface vtx } i). $$

### 7.3.3 Geometric Transformation Swapping Vertices

The first step of the proof of unisolvence in Section 9.3 is also based on a factorization result established in Section 8.2.3. This result involves the first face hyperplane $\mathcal{H}_0^V$ opposite the first vertex $v_0$, but it is more convenient to first prove it for the last reference face hyperplane $\hat{\mathcal{H}}_d$ opposite last reference vertex $\hat{v}_d$ (for $d > 0$).

More generally, it can be useful to transport information on some face of some simplex to any face of any other simplex (see an illustration for $d=3$ in Figure 6). To that purpose, given a permutation $\pi$ of $[0..d]$ (i.e. a bijection from $[0..d]$ onto itself), we consider the permuted geometric transformation $T^{\pi(v)}$, where the notation $\pi(v)$ actually means $(v_{\pi(i)})_{i \in [0..d]}$, and the specialization to the transposition $\tau_d^i$ swapping the last index $d$ with any $i \in [0..d]$. This is formalized as

$$ \textbf{Context } \{d : \text{nat}\}. $$
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Geometric transformation $T_{\tau_d^0(v)}$ in the case $d=k=3$.</figcaption></figure>

The transposition $\tau_d^0$ maps $(0, 1, 2, 3)$ onto $(3, 1, 2, 0)$. The reference simplex $\hat{K}^d$ is mapped onto the generic simplex $K$. The reference vertices are mapped onto the generic vertices with a change of indices: for all $i \in [0..d]$, we have $T_{\tau_d^0(v)}(\hat{v}_i) = v_{\tau_d^0(i)}$. Note that, as $\tau_d^0(d) = 0$, the reference face hyperplane $\hat{H}_d$ (opposite reference vertex $\hat{v}_d$) is mapped onto $H_0^v$ (opposite vertex $v_0$). In this figure, only the visible nodes are depicted.

**Definition T_geometry_permutF** (vtx : [0..d] $\rightarrow$ $\mathbb{R}^d$) (pi_d : [0..d] $\rightarrow$ [0..d]) : $\mathbb{R}^d \rightarrow \mathbb{R}^d := T_geometry (permutF pi_d vtx).

**Definition T_geometry_transpF** vtx (i0 : [0..d]) : $\mathbb{R}^d \rightarrow \mathbb{R}^d := T_geometry_permutF vtx (transp_ord ord_max i0).

where transp_ord implements any transposition and ord_max is defined in Section 4.2.

As expected, when the vertices are affinely independent, the permuted geometric transformation $T^{\pi(v)}$ is a bijective affine mapping ($T_geometry_permutF_bij$), and $T_geometry_permutF_am$), and it maps $\hat{H}_i$ onto $H_{\pi(i)}^v$. We also have, for all $i \in [0..d]$, $\mathcal{L}_i^v = \hat{\mathcal{L}}_{\pi^{-1}(i)} \circ (T^{\pi(v)})^{-1}$. In Rocq,

**Context** $vtx : [0..d] \rightarrow \mathbb{R}^d$.

**Hypothesis** $Hvtx : aff_indep_ms vtx$.

**Context** $pi_d : [0..d] \rightarrow [0..d]$.

**Hypothesis** $Hpi_d : injective pi_d$.

**Let** $pi_d_bij := injF_bij Hpi_d$.

**Let** $pi_d_inv := f_inv pi_d_bij$.

**Lemma T_geometry_permutF_Hface** :

$$ \forall i, \text{ image } (T_geometry_permutF vtx pi_d) (Hface_ref i) = Hface vtx (pi_d i). $$

**Lemma LagPd1_eq_permutF** : $\forall i, \text{ LagPd1 Hvtx } i = (\text{permutF pi_d_inv LagPd1_ref i}) \circ T_geometry_permutF_inv Hvtx Hpi_d$.
where `injF_bij` is a result from the `MathComp` library extrapolating bijectivity from injectivity for functions from a finite type to itself. These results are also available for the specialization to transpositions swapping the last index, where the inverse of the permutation becomes the transposition itself as it is involutive (`T_geometry_transpF_Hface`), and `LagPd1_eq_transpF`.

## 7.4 Geometric Transformation of a Finite Element

The geometric transformation defined in Section 7.2.2 allows to change the geometry of an element. We study now how are affected the other two items in the finite element triple, the approximation space and the degrees of freedom (see Sections 2.4 and 6).

In a more general setting, given types `T1`, `T2`, `A`, consider any function $\phi : T1 \to T2$. Then, `pullback` $\phi := \lambda (f : T2 \to A) \Rightarrow f \circ \phi$ is a function that maps (i.e. pulls back) the functions $T2 \to A$ to $T1 \to A$. When the type `A` is endowed with a module space structure, then `pullback` $\phi$ is linear (`pullback_lms`). Moreover, given a proof $H\phi : \text{bijective} \phi$, then `pullback` $\phi$ is bijective too (`pullback_bij`). Given an additional type `B`, the dual transformation maps (i.e. pushes forward) functions $\sigma$ of type ($T1 \to A \to B$ to ($T2 \to A \to B$), it is defined by `pushforward` $\phi \sigma := \lambda f \Rightarrow \sigma$ (`pullback` $\phi f`). Again, when the types `A` and `B` represent module spaces, `pushforward` $\phi \sigma$ is linear as soon as $\sigma$ is linear (`pushforward_lms`).

Given the spatial dimension `d`, assuming that types `T1`, `T2` are $\mathbb{R}$ (`d` and `A`, `B` are the reals $\mathbb{R}$ (`d` and $\mathbb{R}$, `d`), we want to express the transformation of any input FE into a new FE of the same type, with the transformed vertices.

**Context** `{d : nat}`.

**Context** `{\phi : \mathbb{R}^d \to \mathbb{R}^d}`.

**Hypothesis** $H\phi : \text{bijective} \phi$.

**Variable** `fe_in`: `FE d`.

For that, we need to provide all the items of the new instance `TFE` (for "transformed finite element") of the record `FE` defined in Section 6, proofs included. The fields `shape` and `ndof`, and thus `nvtx`, are kept unchanged. As wanted, the transformed vertices are the images of the input vertices ($v_i^{\text{out}} \stackrel{\text{def.}}{=} \phi(v_i^{\text{in}})$, for all $i \in [0..\text{ndof})$). The transformed approximation space is defined as the preimage of the input approximation space by `pullback`, or equivalently the input approximation space is the pullback of the transformed approximation space ($P^{\text{out}} \stackrel{\text{def.}}{=} \{p^{\text{in}} \circ \phi^{-1} | p^{\text{in}} \in P^{\text{in}}\}$, e.g. see [25, Section 3.4] in the affine case, as below). And the transformed degrees of freedom are the images of the input degrees of freedom by `pushforward`: $\sigma_i^{\text{out}}(f) \stackrel{\text{def.}}{=} \sigma_i^{\text{in}}(f \circ \phi)$, for all $i \in [0..\text{ndof})$, and for all $f: \mathbb{R}^d \to \mathbb{R}$). In `Rocq`,

**Definition** `TFE_KVertices` $\mathbb{R}$: $[0..\text{nvtx fe_in}) \to \mathbb{R}^d := \text{mapF} \phi (\text{KVertices fe_in})$.

**Definition** `TFE_PApprox` $\mathbb{R}$: $\text{FRd d} \to \text{Prop} := \text{preimage} (\text{pullback } \phi) (\text{P}\_\text{approx fe_in})$.

**Definition** `TFE_SDof` $\mathbb{R}$: $[0..\text{ndof fe_in}) \to \text{FRd d} \to \mathbb{R} := \text{mapF} (\text{pushforward } \phi) (\text{S}\_\text{dof fe_in})$.
The proof TFE_S_dof_lms of linearity of the transformed degrees of freedom is basically pushforward_lms applied to S_dof_lms fe_in. The proof TFE_P ciplink has dim that the transformed approximation space is a vector space of the same dimension ndof fe_in follows from the fact that pullback $\phi$ is a bijective linear mapping, and that such a transformation preserves the dimension of vector subspaces (lm has dim rev). And finally, the proof TFE_unisolvence_inj of unisolvence directly follows from the injectivity of pullback $\phi$ and the unisolvence of the input finite element. Then, the transformed FE is

**Definition TFE**: FE d := mk FE (shape fe_in) TFE_K vertices (ndof fe_in) TFE_P approx has dim TFE_S_dof_lms TFE_unisolvence_inj.

In the simplicial case, the transformation function $\phi$ is an affine mapping. It can be fully specified by the vertices of the input finite element and those desired for the transformed FE. Indeed, given an input simplicial FE with affinely independent vertices $v^{in}$, and output vertices $v^{out}$ that are also affinely independent, then the geometric transformations $T^{v^{in}}$ and $T^{v^{out}}$ defined in Section 7.2.2 are both bijective, and the composition $\phi \stackrel{\text{def.}}{=} T^{v^{out}} \circ (T^{v^{in}})^{-1}$ is a bijective affine transformation that maps $v^{in}$ onto $v^{out}$ as illustrated in Figure 7. The corresponding FE transformation TFE_vtx allows to transport any simplicial FE of type FE d to any simplicial geometry.

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Affine transformation of a finite element, from vertices to vertices. This figure explains the chosen $\phi$ in order to construct TFE_vtx.</figcaption></figure>

As a consequence, when the input vertices are the reference vertices $\hat{v}$ of Section 7.1, for any given affinely independent vertices $v$, the transformation function simplifies into $\phi \stackrel{\text{def.}}{=} T^v$, and the corresponding FE transformation TFE_vtx_from_ref allows to transport any reference simplicial FE of type FE d to any simplicial geometry.

The various transformations with their inputs/hypotheses are given in Table 8.

# 8 Multi-indices and Polynomial Space $\mathbb{F}_k^d$

Following the definition of the Rocq record for a finite element of Section 6, we have detailed in the previous section results about the geometry and the vertices. We now
<table><caption>**Table 8** Inputs and hypotheses for the finite element transformation TFE<sup>def</sup>, and its successive specializations TFE_vtx<sup>def</sup> and TFE_vtx_from_ref<sup>def</sup>.</caption><thead><tr><th colspan="2">Input / Hypothesis</th><th>Output</th></tr></thead><tbody><tr><td>fe_in : FE d</td><td>φ : ℝ<sup>d</sup> → ℝ<sup>d</sup> bijective</td><td>TFE : FE d</td></tr><tr><td>fe_in : FE d</td><td>fe_in is a simplex</td><td>TFE_vtx : FE d</td></tr><tr><td>v<sup>in</sup> affine independent</td><td>v<sup>out</sup> affine independent</td><td>(with φ = T<sup>v<sup>out</sup></sup> ◦ (T<sup>v<sup>in</sup></sup>)<sup>-1</sup>)</td></tr><tr><td>fe_in : FE d</td><td>fe_in is a simplex</td><td>TFE_vtx_from_ref : FE d</td></tr><tr><td>v<sup>in</sup> = v̂</td><td>v affine independent</td><td>(with φ = T<sup>v</sup>)</td></tr></tbody></table>

focus on the approximation space. Indeed, we need to define a finite dimensional function space, in which the approximate solution is sought in each element. The most common one is the space of polynomials of bounded total degree, that we define in Section 8.2. The subtlety here is that we have multi-dimensional polynomials, that is to say polynomials from $ℝ^d$ to $ℝ$. We rely on multi-indices defined in Section 8.1 for defining one basis of multi-dimensional polynomials of bounded degree. The definition and properties of these multi-indices is key to the following, as they are also used for numbering the nodes in Section 9.1.

## 8.1 Multi-indices

Given a dimension $d \in \mathbb{N}$, a multi-index $\alpha$ is a family of size $d$ of natural numbers, that is an element of $\mathbb{N}^d$. Its sum is denoted by $|\alpha| \stackrel{\text{def.}}{=} \sum_{i=0}^{d-1} \alpha_i \in \mathbb{N}$. Given a maximal sum value $k \in \mathbb{N}$, $A_k^d$ is defined in Section 2.3.4 as the set of multi-indices of size $d$ and of sum at most $k$, that is (8):

$$A_k^d \stackrel{\text{def.}}{=} \{\alpha \in \mathbb{N}^d \mid |\alpha| \le k\}.$$

This is a multifaceted object. This is basically a combinatorics object and defined as is. It is also related to a monomial (see (9)), so that all the defined monomials form a basis of the $d$-dimensional polynomials of degree bounded by $k$ (see Section 8.2).

It has also a geometric interpretation as the indices of the nodes in the simplicial Lagrange finite elements. More precisely, see Figures 2 and 9. The Lagrange nodes are the dots of any color; they are defined such that they are evenly distributed in the simplex. Their formal definition is given in Section 9.1.1, but let us give here an intuition (see (25) below). Each reference node $\hat{a}_{\alpha}$ has coordinates for $i \in [0..d-1]$ $(\hat{a}_{\alpha})_i = \frac{\alpha_i}{k}$. It is thus natural to represent a Lagrange node in a simplex by its multi-index $\alpha$ in $A_k^d$. The correspondences between a multi-index, a node and a monomial are used throughout the document.

Section 8.1.1 formally defines $A_k^d$ (all the families of size $d$ of sum smaller or equal to $k$) and $C_k^d$ (all the families of size $d$ of sum equal to $k$, see (19)) as functions. Section 8.1.2 focuses on the chosen order for $A_k^d$. Section 8.1.3 defines an inverse for the function Adk. Section 8.1.4 defines various multi-indices bijections related to faces.
<figure><img src="image_9.png" /><figcaption>Fig. 9 Lagrange nodes $(\hat{\mathbf{a}}_{\alpha})_{\alpha \in \mathcal{A}_k^d}$ of the reference simplex $\hat{\mathbf{K}}_d$ when $d \in \{2, 3\}$ and $k = 3$.</figcaption></figure>

Each node is depicted as a colored ball, and corresponds to a unique element of $\mathcal{A}_3^d$. The colors correspond to degrees $l \le 3$ of polynomials, or equivalently to the sizes of multi-indices (i.e. in $C_l^d$ for $l \le 3$). In magenta, the node $\hat{\mathbf{a}}_0$ corresponds to constant polynomials (with degree 0) in $\mathbb{P}_0^d$, and to the multi-index $0$ in the singleton $C_0^d$. In green, the nodes correspond to nonconstant affine polynomials (with degree 1), and to the multi-indices $\delta_1, \dots, \delta_d$ in $C_1^d$. In red, the nodes correspond to non-affine quadratic polynomials (with degree 2), and to multi-indices in $C_2^d$. In blue, the nodes correspond to nonquadratic cubic polynomials (with degree 3), and to multi-indices in $C_3^d$. We observe in this picture that $\mathcal{A}_3^d = C_0^d \uplus C_1^d \uplus C_2^d \uplus C_3^d$.

### 8.1.1 Formal Definition of Multi-indices

Now let us formalize of this definition. Even if it is a mathematical set, we formalize it as a family of multi-indices of type $[0..\text{pbinom d k}] \to \mathbb{N}^d$, recall the definition of $\text{pbinom}$ in Section 4.5. It is indeed known that the set $\mathcal{A}_k^d$ is of size $\binom{d}{d+k}$, so we choose to define $\mathcal{A}_k^d$ as a family of size $\binom{d}{d+k}$ of multi-indices of size $d$. The family $\mathcal{A}_k^d$ is defined by layers: the multi-indices of sum 0, then of sum 1, and so on until the multi-indices of sum $k$. More precisely,

$$ \mathcal{C}_k^d \stackrel{\text{def.}}{=} \{\alpha \in \mathbb{N}^d \mid |\alpha| = k\}, \quad \text{and} \quad \mathcal{A}_k^d = \bigcup_{l=0}^{k} C_l^d. \qquad (19) $$

In words, $\mathcal{C}_k^d$ is the family of multi-indices of size $d$ and of sum equal to $k$. Then, $\mathcal{A}_k^d$ is the (disjoint) union of $(\mathcal{C}_l^d)_{l \in [0..k]}$, built in Rocq by concatenation of the $\mathcal{C}_l^d$'s.

Now let us formalize the definition of $\mathcal{C}_k^d$. A difficulty is the handling of $d = 0$. In fact, $\mathcal{C}_k^0$ does not exist (except for $k = 0$), as there is no vector of size 0 and of nonzero sum. We therefore only consider $\mathcal{C}_k^{d_1+1}$, that is the set of vectors of size $d_1 + 1$ and of sum $k$. Its type is therefore $[0..\text{pbinom d1 k}] \to \mathbb{N}^d \cdot 1$.

It is constructed as a fixpoint on the dimension $d$. For $d = 0$, there is a single vector of size 1 and of sum $k$ that is $\{k\}$. Let us assume we have computed all the $\mathcal{C}_i^d$'s. Then,
the family $C_k^{d+1}$ is made of multi-indices of size $d+1$. Let $\alpha \in C_k^{d+1}$; then we split it as $\alpha = (\alpha_0, \alpha')$ and we have that $0 \le \alpha_0 \le k$ and $\alpha_0 + |\alpha'| = k$. So we deduce that the set $C_k^{d+1}$ is the (disjoint) union, for $0 \le i \le k$, of $(k-i, \alpha')$ for $\alpha' \in C_i^d$. By using the `concatnF` function of Section 4.3.1, we define `CSdk` that is $C_k^{d+1}$,

Fixpoint `CSdk` (d1 k : nat) : [0..pbinom d1 k] $\rightarrow$ 'nat^d1.+1 :=
match d1 with
  | 0 $\Rightarrow$ $\lambda_{\_\_}$ $\Rightarrow$ k
  | S d2 $\Rightarrow$ castF_(concatnF ($\lambda$ i $\Rightarrow$ mapF (insertF0 (k-i)) (CSdk d2 i)))
end.

The `insertF0` function only concatenates a value (here $k-i$) as a first element in a family of multi-indices of size $s$, making a family (of the same size) of multi-indices of size $s+1$. This construction is also illustrated in Figure 10.

<figure><img src="image_10.png" /><figcaption>**Fig. 10** Partition $C_k^d = \cup_{i=0}^k \{(k-i, \beta) \mid \beta \in C_i^{d-1}\}$, with $d=k=3$.</figcaption></figure>

Geometrically, $\{(k-i, \beta) \mid \beta \in C_i^{d-1}\}$ is the intersection of the hyperplane of equation $\alpha_0 = k-i$ (depicted in light pink) and of $C_k^d$ (of equation $\alpha_0 + \alpha_1 + \alpha_2 = k$).

The reference Lagrange node $\hat{a}_{(\alpha_1, \alpha_2, \alpha_3)}$ belonging to the face $\mathcal{H}_0$ is in blue and corresponds to the element $(\alpha_1, \alpha_2, \alpha_3) \in C_k^d$. For instance, the set $\{(1, 2, 0), (1, 1, 1), (1, 0, 2)\}$ is depicted by the blue nodes linked by a dashed arrow.

Note that the set $C_k^d$ defined in (19) may have been defined otherwise. For instance, we could have concatenated the $(i, \beta)$ for $\beta \in C_{k-i}^{d-1}$, or the $(\beta, k-i)$ for $\beta \in C_i^{d-1}$, and so on. All those are mathematically correct and the sets are equivalent but do not respect the same numbering, hence the same ordering (see Section 8.1.2 for more on ordering).

The `Rocq` definition of $A_k^d$ is then straightforward, it is the concatenation of the $C_l^d$ as in (19),

Definition Adk (d k : $\mathbb{N}$) : '$\mathbb{N}^{\{(pbinom d k).+1,d\}}$ :=
match d with
  | 0 $\Rightarrow$ 0
  | S d1 $\Rightarrow$ castF_(concatnF (CSdk d1))
end.

The underscore is here a rather complex proof (compared to most proofs inside casts) as it is `binom_rising_sum_1` (see (13) in Section 4.5) to ensure that the sum of the size of the $C_i^{d'}$'s makes the expected size of $A_k^d$.

Interesting results include the values for $k = 1$ (`Ad1_eq_2`) and $d = 1$ (`A1k_eq_2`), and the fact that the last layer of $A_k^d$ corresponds to $C_k^d$, that is to say the maximal sum is at the end of the family,<fnref n="19" />

**Lemma** Adk_last layer: $\forall d k idk, d \neq 0 \rightarrow k \neq 0 \rightarrow \Sigma (\text{Adk d k idk}) = k \leftrightarrow (\text{pbinom d k.} - 1) + 1 \leq idk$.

## 8.1.2 Ordering Multi-indices

There are several ways to construct $A_k^d$ when seen as a set. When seen as a family, we need to solve the ordering of these multi-indices: in which order do we enumerate them? And this gives us the way we build them. Even if our design choices have been explained in Section 8.1.1, the explanation lies with the multi-index ordering.

Another point of view on a multi-index is the exponent of a monomial in dimension $d$, recall (9). Ordering multi-indices can therefore similarly be seen as ordering monomials. This is a known problem with many solutions, that is to say many possible orderings, e.g. see [48, Chap. 2, Sec. 2]. A common ordering is the lexicographic order and another is the graded reverse lexicographic order (a.k.a. grevlex order), often used when Gröbner bases are concerned. But neither suited us for reasons explained later in Section 8.1.4. When the implementation of the FEM will be tackled, we may need other orderings, but program verification is out of the scope of this paper.

**Ordering Definition: Grsymlex.**

Let us now define the chosen order, graded symmetric lexicographic, or simply “grsymlex” order,

$$ \alpha <^{\text{grsymlex}} \beta \quad \overset{\text{def.}}{\iff} \quad \begin{cases} |\alpha| < |\beta|, & \text{or} \\ |\alpha| = |\beta| & \wedge \beta <^{\text{lex}} \alpha. \end{cases} \qquad (20) $$

This amounts to first compare the sum of multi-indices, and in case of equality, use the symmetric of the lexicographic order. In other words, when $|\alpha| = |\beta|$, we have $\alpha <^{\text{grsymlex}} \beta$ iff $\beta_i < \alpha_i$ for the first index $i$ where $\alpha_i$ and $\beta_i$ differ. An equivalent recursive point of view is

$$ \alpha <^{\text{grsymlex}} \beta \iff \begin{cases} |\alpha| < |\beta|, & \text{or} \\ |\alpha| = |\beta| & \wedge d \geq 2 \wedge \alpha' <^{\text{grsymlex}} \beta' \\ \text{where we denote } \alpha = (\alpha_0, \alpha') \text{ and } \beta = (\beta_0, \beta'). \end{cases} $$

This is defined for a generic order $<$, and we use it with `lt`, the strict order on $\mathbb{N}$.

The Rocq definition is exactly (20), but we have also the previous equivalence on $\mathbb{N}$,

<footnote n="19">¹⁹Note that when positivity of natural numbers is not needed structurally, we choose to use non-cancellation instead of positiveness.</footnote>
**Lemma grsymlex_lt_S[‡] :** $\forall \{n: nat\} (x y : ^{nat}n.+1)$,
$grsymlex_lt x y \leftrightarrow \Sigma x < \Sigma y \lor (\Sigma x = \Sigma y \land grsymlex_lt (skipF0 x) (skipF0 y))$.

We prove that this is a strict and total order (`grsymlex_lt_strict_total_order[‡]`), and a monomial order (`grsymlex_br_PLUSCompAT1[‡]`). This ordering is shown in Figure 11 on the $A_k^d$ where examples of values ordered by grsymlex are provided.

<figure><img src="image_11.png" /><figcaption>**Fig. 11** Grsymlex ordering of $A_k^d$ when $d \in \{2, 3\}$ and $k = 3$. The increase in the order is represented by dashed arrows (when possible). For $A_3^2$ ($d=2$), we have $(0,0) < \text{grsymlex (1,0)} < \text{grsymlex (0,1)} < \text{grsymlex (2,0)} < \text{grsymlex (1,1)} < \text{grsymlex (0,2)} < \text{grsymlex (3,0)} < \text{grsymlex (2,1)} < \text{grsymlex (1,2)} < \text{grsymlex (0,3)}$. When $d \ge 3$, we have for instance for $C_3^3$ $(3,0,0) < \text{grsymlex (2,1,0)} < \text{grsymlex (2,0,1)} < \text{grsymlex (1,2,0)} < \text{grsymlex (1,1,1)} < \text{grsymlex (1,0,2)} < \text{grsymlex (0,3,0)} < \text{grsymlex (0,2,1)} < \text{grsymlex (0,1,2)} < \text{grsymlex (0,0,3)}$.</figcaption></figure>

Ordering $A_k^d$.

There is left to prove that the construction of the $A_k^d$ and the grsymlex ordering match.

**Lemma Adk_sortF[‡] :** $\forall d k$, sortedF `grsymlex_lt` (`Adk d k`).

where `sortedF[‡] leT A` is the predicate stating that the family `A: T^n` is ordered according to the binary relation `leT: T → T → Prop`.

As $A_k^d$ is built by concatenation of the $C_l^{d'}$'s, we rely on the `concatN_F_order[‡]` lemma of Section 4.3.1. Given the properties of the order, we have two things left to prove. First, the last element of $C_l^d$ is smaller than the first element of $C_{l+1}^d$ for $l < k$. As the first one has sum $l$ and the second one has sum $l+1$, they are ordered by grsymlex. This is due to the fact that $A_k^d$ is built by layer of increasing sum. The second and more complex proof is that each $C_l^d$ is ordered. As the definition of $C_k^d$ is by induction on $d$, so is this proof. By using again the `concatN_F_order` lemma, the induction and an easy result on `insertF0` provides the fact that the families to be concatenated are sorted. The hard part is that the families are enumerated in the right order, and this is due to our choice to concatenate the $(k-i, \beta)$ for $\beta \in C_i^{d+1}$ for $i$ from 0 to $k$.
### 8.1.3 Bijection and Multi-indices

Now let us focus on the correctness of $\mathcal{A}_k^d$, that returns multi-indices of size $d$ and of sum smaller or equal to $k$. There is left to prove that all such multi-indices are provided and that $\mathcal{A}_k^d$ makes a bijection between the integers from 0 to $\binom{d}{d+k}-1$ and the multi-indices of size $d$ and of sum smaller or equal to $k$.

**Lemma Adk_inj**: $\forall d k$, injective (Adk d k).

**Lemma Adk_surj**: $\forall d k$ ($b: \mathbb{N}^\dagger d$), $\sum b \le k \to \{\text{idk} | b = \text{Adk d k idk}\}$.

The injectivity is deduced from the previous lemma: as $\mathcal{A}_k^d$ is sorted with a strict order, it is injective. The surjectivity (all multi-indices of sum smaller or equal to $k$) belong to $\mathcal{A}_k^d$ relies on the surjectivity of $\mathcal{C}_k^d$, that is done by induction. There is no Rocq specific difficulty except a tedious handling of casts.

From this constructive definition, we then define an inverse of Adk that takes a multi-index of size $d$ and returns an integer between 0 and $\binom{d}{d+k}-1$,

**Definition Adk_inv**: $(d k : \mathbb{N}) : \mathbb{N}^\dagger d \to [0..\text{pbinom d k}] :=$
$\lambda b \Rightarrow \text{match le_dec}(\Sigma b) k \text{ with}$
$\quad \left| \begin{array}{l} \text{left H} \Rightarrow \text{proj1_sig}(\text{Adk_surj d k b H}) \\ \text{right } \Rightarrow \text{ord0} \end{array} \right.$
$\text{end.}$

We choose to make it a total function, that takes any multi-index. When the sum of the multi-index is too large, a default value (zero) is returned. We proved the correctness of Adk_inv, with a hypothesis on the multi-index when needed,

**Lemma Adk_inv_correct_r**: $\forall d k b$, $\sum b \le k \to \text{Adk d k (Adk_inv d k b)} = b$.

**Lemma Adk_inv_correct_l**: $\forall d k \text{idk}$, $\text{Adk_inv d k (Adk d k idk)} = \text{idk}$.

### 8.1.4 Face Multi-indices

As multi-indices are related to Lagrange nodes (see Figure 9 and Section 9.1), it is useful to characterize the multi-indices associated with the nodes on the faces and to construct bijections similar to the geometric mappings of Section 7.2.3. When $d > 0$, for $i \in [0..d]$, there exist bijections between the family of multi-indices $\mathcal{A}_k^{d-1}$ and either the layer $\mathcal{C}_k^d$ (case $i=0$), and the multi-indices of $\mathcal{A}_k^d$ with a zero component.

More precisely, we define

$$ \mathcal{A}_{k,0}^d \stackrel{\text{def.}}{=} \mathcal{C}_k^d, \quad \text{and} \quad \forall i \in [1..d], \quad \mathcal{A}_{k,i}^d \stackrel{\text{def.}}{=} \{\alpha \in \mathcal{A}_k^d | \alpha_{i-1} = 0\}. \quad (21) $$

The Rocq definition is the property `in_A_Sdki` on the index `idk` such that `Adk d k idk` belongs to $\mathcal{A}_{k,i}^d$. Then, the aforementioned mappings are defined for all $\alpha' \in \mathcal{A}_k^{d-1}$ as

$$ \mathbb{1}_{k,0}^{d}(\alpha') \stackrel{\text{def.}}{=} (k - |\alpha'|, \alpha') \in \mathcal{A}_{k,0}^{d}, \quad (22) $$

$$ \forall i \in [1..d], \quad \mathbb{1}_{k,i}^{d}(\alpha') \stackrel{\text{def.}}{=} (\alpha'_0, \dots, \alpha'_{i-2}, 0, \alpha'_{i-1}, \dots, \alpha'_{d-2}) \in \mathcal{A}_{k,i}^{d}. \quad (23) $$

These functions can be straightforwardly formalized in Rocq. Nevertheless, it is more convenient for us to use indices instead of multi-indices. More precisely, we
define `inj_ASdki` of type `[0..pbinom d1 k] → [0..pbinom d k]` that takes the index of $\alpha'$ and provides the index of the result of either $\|_{k,0}^d$ or $\|_{k,i}^d$ (22) and (23) are merged as a single definition). Note that this reminds (17) and (18) related to faces. This is exploited later in Section 9.1.2 that details in which faces are the Lagrange nodes.

Now, we have all the definitions to give and ensure the ordering requirements, that are the reasons we choose the `grsymlex_lt` order. We want the order to be:

(i) a monomial order, i.e. a total order that is compatible with the monoid structure of the monomials, that is for all $\alpha, \beta, \gamma \in \mathcal{A}_k^d$, $X^\alpha < X^\beta$ implies $X^\alpha X^\gamma < X^\beta X^\gamma$, which means that $\alpha < \beta$ implies $\alpha + \gamma < \beta + \gamma$. It is also generally required that for all $\alpha \neq 0$, $1 = X^0 < X^\alpha$, that is $0 < \alpha$. It is proved in Section 8.1.2.

(ii) increasing with respect to the degree (from $k$ to $k+1$), i.e. multi-indices of sum at most $k$ should be numbered before those of sum $k+1$, that is for all $\alpha \in \mathcal{A}_k^d$ and for all $\beta \in \mathcal{C}_{k+1}^d$, we should have $\alpha < \beta$. It allows us to easily sort monomials with respect to their total degree and is proved by Lemma **Adk_last_layers** in Section 8.1.1.

(iii) increasing with respect to the dimension (from $d-1$ to $d$), i.e. “natural” bijections between ($d-1$)-multi-indices of sum at most $k$ and $d$-multi-indices of sum $k$ or with a zero entry ($\|_{k,i}^d$ for $i \in [0..d]$ defined in Section 8.1.4)), should be increasing. It helps to relate the face nodes and the volume nodes during the computations and is proved by Lemma **inj_ASdki_monot**.

(iv) the multi-indices corresponding to the reference vertices $\hat{v}$ (see Section 7.1) are increasing, i.e. $0 < k\delta_0 < k\delta_1 < \dots < k\delta_{d-1}$. It allows to have positive simplices (i.e. positively oriented) when computing the integrals. This is formally proved by two lemmas. First, `grsymlex_lt_z1` states that $\alpha \neq 0$ implies $0 < \alpha$. Second, we prove that $k\delta$ is sorted for the grsymlex order (for a positive $k$).

**Lemma grsymlex_sortedF_itemF:**
$\forall \{n\} k, k \neq 0 \rightarrow \text{sortedF grsymlex_lt } (\lambda i \Rightarrow \text{itemF n i k})$.

## 8.2 Polynomial Space $\mathbb{F}_k^d$

In this section, we formalize the polynomial space $\mathbb{F}_k^d$ defined in (10) and its properties. It includes all the polynomials on $\mathbb{F}^d$ with total degree at most $k$. In Section 8.2.1, we define this subset and its first properties, including it is a finite-dimensional vector space and its dimension. Section 8.2.2 is about specific values of $d$ and $k$. Section 8.2.3 is about the factorization of polynomials.

### 8.2.1 $\mathbb{F}_k^d$ Definition and Basis

We have defined in Section 8.1 the multi-indices $\mathcal{A}_k^d$, that are all the families of size $d$ with sum at most $k$. We therefore consider each of these multi-index as the exponent of a monomial in order to construct each monomial. The polynomials are then all the linear combinations of these monomials.

Following (9), monomials are first formalized as

**Definition Monom_dk:** $d k : [0..pbinom d k] \rightarrow \text{FRd d} := \lambda \text{idk} x \Rightarrow \text{powF_P}(\text{Adk d k idk}) x$.
where `powF_P` is the product of univariate monomials. We remind that `FRd` is the vector space of functions mapping from the real space $\mathbb{R}^d$ to $\mathbb{R}$. Therefore, `Monom_dk` is a family of such functions. The size of this family is the cardinal of $\mathcal{A}_k^d$, that is $(\binom{d+k}{k})$, i.e. $(\text{pbinom}_d k).+1$. Then, following (10), $\mathbb{P}_k^d$ is formalized as their linear span,

**Definition** Pdk d k : $FRd \to \text{Prop} := \text{lin\_span} (\text{Monom\_dk d k})$.

As a linear span, it is trivially a vector space. We proved many lemmas, including

**Lemma** Pdk\_pullback\_am : $\forall \{d n\} k (p : FRd \to \mathbb{R}^n \to \mathbb{R}^d)$, $aff\_map\_ms f \to Pdk \ d \ p \to Pdk \ n \ k$ (pullback $f \ p$).

The composition of a polynomial of degree at most $k$ with an affine mapping is also a polynomial of degree at most $k$. Note that the affine mapping $f$ is from $\mathbb{R}^n$ to $\mathbb{R}^d$, therefore $p$ is a polynomial on $\mathbb{R}^d$ while $p \circ f$ is on $\mathbb{R}^n$. The proof relies on two lemmas. First, the product of polynomials has degree at most the sum of the degrees of the polynomials. And second, an affine mapping to the power $n$ has degree at most $n$. In the previous statement, one should understand that “degree $k$” means “belongs to $\mathbb{P}_k^d$”. We have not defined the degree *per se*, as the nicest definition of degree gives the value $-\infty$ to the zero polynomial and we did not really need it.

The most complex proof of this section is that the monomials `Monom_dk d k` form a linearly independent family. This is required to prove that it is a basis of $\mathbb{P}_k^d$, and then its cardinal provides the value of the dimension.

We could have used a formal derivative (on the sequences of coefficients), but it was simpler to rely on “true” derivatives. More precisely, we rely on the derivative of real functions of a real variable provided by the `Coquelicot` library [38]. We first define the derivative with respect to a given variable index $i \in [0..d-1]$, then an iterated derivative and finally the derivative with respect to a multi-index $\alpha$. The next challenge is then to switch the derivatives with respect to several variables. This holds when the function is regular enough, but a generic theory of multi-variable differentiation was out of scope of this article, so we decided to only derive polynomials (of any degree). On this subset, we proved we could switch iterative derivatives at will and compute any derivative with respect to a multi-index. Actually, we do not need the real value of the derivative, only whether it is zero or not.

We proved that, for all $\alpha, \beta \in \mathbb{N}^d$, there exists $C_{\alpha,\beta} \in \mathbb{R}$ such that

$$ \partial^{\beta} \mathbf{X}^{\alpha} \stackrel{\text{def.}}{=} \frac{\partial^{|\beta|} \mathbf{X}^{\alpha}}{\partial x_0^{\beta_0} \dots \partial x_{d-1}^{\beta_{d-1}}} = C_{\alpha, \beta} \mathbf{X}^{\alpha-\beta} \quad \text{and} \quad (C_{\alpha, \beta} = 0 \iff \exists i, \alpha_i < \beta_i), \quad (24) $$

where $\alpha - \beta$ denotes the family $(\alpha_i - \beta_i)_{i \in [0..d)}$, with the subtraction being the `Rocq` total function `sub` on $\mathbb{N}$, that returns 0 when $\alpha_i < \beta_i$ (but in this case $\mathbf{X}^{\alpha-\beta}$ is useless because $C_{\alpha,\beta} = 0$).

Now let us go back to prove that $(\mathbf{X}^{\alpha})_{\alpha \in \mathcal{A}_k^d}$ is linear independent. Let us assume that we have a family of scalars $\mathbf{L}$ such that $\sum_{\alpha \in \mathcal{A}_k^d} L_{\alpha} \mathbf{X}^{\alpha} = 0$ (a function equality). Given any $\beta \in \mathcal{A}_k^d$, we have to prove that $L_{\beta} = 0$. By linearity of differentiation with
respect to $\beta$, we obtain $\sum_{\alpha \in A_k^d} L_\alpha \partial^\beta \mathbf{X}^\alpha = 0$. From (24), we deduce that

$$ \sum_{\substack{\alpha \in A_k^d \\ \forall i, \beta_i \le \alpha_i}} L_\alpha C_{\alpha, \beta} \mathbf{X}^{\alpha-\beta} = L_\beta C_{\beta, \alpha} + \sum_{\substack{\alpha \in A_k^d \\ \forall i, \beta_i \le \alpha_i \\ \exists j, \beta_j < \alpha_j}} L_\alpha C_{\alpha, \beta} \mathbf{X}^{\alpha-\beta} = 0. $$

Now let us evaluate these functions in zero, we get

$$ L_\beta C_{\beta, \beta} + \sum_{\substack{\alpha \in A_k^d \\ \forall i, \beta_i \le \alpha_i \\ \exists j, \beta_j < \alpha_j}} L_\alpha C_{\alpha, \beta} \mathbf{0}^{\alpha-\beta} = L_\beta C_{\beta, \beta} = 0. $$

As $C_{\beta, \beta} \ne 0$, we deduce that $L_\beta = 0$ and that the family $(\mathbf{X}^\alpha)_{\alpha \in A_k^d}$ is linearly independent. Since $\mathbb{P}_k^d$ is the span of the monomials, it is then straightforward that they form a basis, and the dimension is the cardinal of $A_k^d$,

Lemma Monom_dk_basis: basis (Pdk d k) (Monom_dk d k).

Lemma Pdk_has_dim: has_dim (Pdk d k) (pbinom d k). + 1.

### 8.2.2 $\mathbb{P}_k^d$ for Small Values of $d$ and $k$

Now that we have defined a generic $\mathbb{P}_k^d$, it is useful to focus on the small values of $d$ and $k$ for specific properties and links with the previously-defined polynomials.

#### Small $d$.

First, for $d = 0$, $\mathbb{R}^0$ being the unit type $\{0\}$ (see Section 4.3.1), all functions of type $\mathrm{FRd}$ 0 are constant functions, and thus they are also polynomials.

Lemma P0k_eq: $\forall k$ ($p : \mathrm{FRd} 0$), $\mathrm{Pdk} 0 k p \leftrightarrow p = \lambda \_ \Rightarrow p 0$.

Lemma P0k_full: $\forall k$, $\mathrm{Pdk} 0 k = \mathrm{fullset}$.

Then, for $d = 1$, we have the standard polynomials of a single variable. In addition to the canonical basis Monom_dk, we define another basis [21, p.8]. Let $\mathbf{a}$ be a family of $k+1$ (pairwise distinct) reals called nodes. Then, the Lagrange polynomial basis, called LagP1k in Rocq, is mathematically defined as

$$ \forall i \in [0..k], \forall x \in \mathbb{R}, \quad \mathbb{L}_i(x) \stackrel{\text{def.}}{=} \prod_{j \in [0..k] \setminus \{i\}} \frac{x - a_j}{a_i - a_j}. $$

The Rocq definition is a total function that does not require the injectivity of the nodes $\mathbf{a}$. However, the interesting properties hold when they are pairwise distinct: each $\mathbb{L}_i$ evaluates to 1 at its corresponding node $a_i$, and to 0 at all others (LagP1k_kron); these polynomials are linearly independent and span the polynomial space $\mathbb{P}_k^1$, forming a basis (LagP1k_basis). Moreover, the components of any polynomial on this basis are its values taken at the nodes,
**Hypothesis** $Hnode : injective node$.

**Lemma** LagP1k_decomp: $\forall \{p\}, Pdk \ 1 \ k \ p \to p = lin\_comb (mapF \ p \ node) (LagP1k \ node)$.

### Small $k$.

First, for $k = 0$, zero degree polynomials are constant functions, but contrarily to the case $d = 0$, they do not span the entire set of functions.

**Lemma** Pd0_eq: $\forall \{d\} (p : FRd \ d), Pdk \ 0 \ p \leftrightarrow p = \lambda \Rightarrow p \ 0$.

The most interesting case is $k = 1$. It corresponds to the affine mappings,

**Lemma** Pd1_am_1 Equiv: $\forall \{d\} (f : FRd \ d), Pdk \ 1 \ f \leftrightarrow aff\_map\_ms \ f$.

We also provide two other bases for $Pdk \ d \ 1$. More precisely, the reference Lagrange polynomials $\hat{\mathcal{L}}$ of Section 7.2.1 (in Rocq, LagPd1_ref) form a basis of $Pdk \ d \ 1$. It is also the case for the generic Lagrange polynomials $\mathcal{L}^v$ of Section 7.2.3 (in Rocq, LagPd1), provided the vertices are affinely independent. Moreover, just as for $d = 1$, the components of any polynomial on this latter basis are its values taken at the vertices, that are also the nodes in this case (see Section 9.1.1),

**Hypothesis** Hvtx : aff\_indep\_ms \ vtx.

**Lemma** LagPd1_decomp_node: $\forall \{p\}, Pdk \ 1 \ p \to p = lin\_comb (mapF \ p \ (node \ vtx)) (castF \ _ (LagPd1 \ Hvtx))$.

## 8.2.3 Factorization of Polynomials

Polynomials vanishing on a face hyperplane are shown to be factorizable by the corresponding Lagrange polynomial. These results are based on the Euclidean division of any polynomial by a monomial.

For $d, k \in \mathbb{N}$, given a polynomial $p$ in $\mathbb{P}_{k+1}^{d+1}$, its Euclidean division by the monomial $X_d$ writes $p = p_0 + X_d * p_1$ where $p_0$ belongs to $\mathbb{P}_{k+1}^d$, i.e. does not depend on the last variable $x_d$, and $p_1$ belongs to $\mathbb{P}_k^{d+1}$, i.e. has a total degree one less than that of $p$. More precisely, $PSdSk\_split$ states that,

$$ \forall p \in \mathbb{P}_{k+1}^{d+1}, \exists p_0 \in \mathbb{P}_{k+1}^d, \exists p_1 \in \mathbb{P}_k^{d+1}, \forall \mathbf{x} \in \mathbb{R}^{d+1}, p(\mathbf{x}) = p_0(\tilde{\mathbf{x}}) + x_d * p_1(\mathbf{x}), $$

where $\mathbf{x} = (\tilde{\mathbf{x}}, x_d)$ with $\tilde{\mathbf{x}} \in \mathbb{R}^d$. The proof path we followed was not found in textbooks, and is original up to our knowledge. The idea is first that the existence of these two polynomials is transported by linear combination. As $\mathbb{P}_{k+1}^{d+1}$ is the span of monomials, there is left to prove that it applies to monomials. For $\alpha \in \mathcal{A}_{k+1}^{d+1}$, we have two cases. If $\alpha_d = 0$, then we pose $p_0 = \mathbf{X}^\alpha$ and $p_1 = 0$. If $\alpha_d > 0$, then we pose $p_0 = 0$ and $p_1 = \mathbf{X}^{\alpha-\delta_d}$.

This lemma is interesting as it could have been the definition of $\mathbb{P}_k^d$. Instead of defining $\mathbb{P}_k^d$ as the span of the monomials (see Section 8.2.1) defined from the multi-indices, we could have defined $\mathbb{P}_k^d$ by induction as $\mathbb{P}_{k+1}^{d+1} = \mathbb{P}_{k+1}^d + X_d * \mathbb{P}_k^{d+1}$, with initializations from Section 8.2.2.
When $d > 0$, the following lemma states that any polynomial of degree at most $k+1$ (i.e. in $\mathbb{P}_{k+1}^d$) cancels on the last reference face hyperplane $\hat{\mathcal{H}}_d$ (see Section 7.3.1) iff it factors into $\hat{\mathcal{L}}_d * q$ for some polynomial $q$ of degree at most $k$ (i.e. in $\mathbb{P}_k^d$). More formally,

$$ \forall p \in \mathbb{P}_{k+1}^d, \quad (\forall \mathbf{x} \in \hat{\mathcal{H}}_d, p(\mathbf{x}) = 0) \quad \Longrightarrow \quad \exists q \in \mathbb{P}_{k+1}^d, p = \hat{\mathcal{L}}_d * q. $$

In Rocq,

**Lemma factor_zero_on_last_Hface_ref** $\square$: $\forall \{p\}$, Pdk d k.+1 p $\rightarrow$ $(\forall x, Hface\_ref\_ord\_max x \rightarrow p x = 0) \rightarrow \exists q$, Pdk d k q $\wedge$ p = LagPd1\_ref\_ord\_max * q.

The proof relies on the Euclidean division PdSk_split above. Given $p \in \mathbb{P}_{k+1}^d$, there exists $\tilde{p}_0 \in \mathbb{P}_k^{d-1}$ and $q \in \mathbb{P}_k^d$ such that $p = \tilde{p}_0 + X_{d-1} * q = \tilde{p}_0 + \hat{\mathcal{L}}_d * q$. Thus, for any $\tilde{x} \in \mathbb{R}^{d-1}$, using $Hface\_ref\_is\_Ker$ $\square$ from Section 7.3.1, we have

$$ \mathbf{x} \stackrel{\text{def.}}{=} (\tilde{x}, 0) \in \hat{\mathcal{H}}_d = \ker \hat{\mathcal{L}}_d \subset \mathbb{R}^d \quad \text{and} \quad 0 = p(\mathbf{x}) = \tilde{p}_0(\tilde{x}) + 0 * q(\mathbf{x}). $$

Hence, $\tilde{p}_0 = 0$ and $p = \hat{\mathcal{L}}_d * q$. $\square$

The slight technical point is the handling of the injection and projection between $\mathbb{R}^{d-1}$ and $\mathbb{R}^d$, respectively using `insertF_0 ord_max` and `widenF_S` (see Section 4.2).

The previous result is actually valid for any reference face hyperplane, but such generality is more profitable on a generic element. Given vertices `vtx` and a proof `Hvtx`: `aff_indep_ms` `vtx` of their affine independence, a similar result holds for cancellation on any face hyperplane of the generic element,

**Lemma factor_zero_on_Hface** $\square$: $\forall \{p\}$ i, Pdk d k.+1 p $\rightarrow$ $(\forall x, Hface\_vtx i x \rightarrow p x = 0) \rightarrow \exists q$, Pdk d k q $\wedge$ p = LagPd1\_Hvtx i * q.

The proof (both in paper and in Rocq) uses the affine geometric mapping $T^{\tau_d^i(v)}$ that maps the references vertices $\hat{v}$ onto the transposition of vertices $v$ where the vertex $v_i$ is swapped with the last one $v_d$ (see Section 7.3.3). Given $p \in \mathbb{P}_{k+1}^d$ vanishing on $\mathcal{H}_i^v$, let $\hat{p} \stackrel{\text{def.}}{=} p \circ T^{\tau_d^i(v)}$. Then, $\hat{p} \in \mathbb{P}_{k+1}^d$ because $T^{\tau_d^i(v)}$ is affine, and it vanishes on $\hat{\mathcal{H}}_d$ because $T^{\tau_d^i(v)}(\hat{\mathcal{H}}_d) = \mathcal{H}_i^v$ (see Section 7.3.3). Hence, the previous result states the existence of $\hat{q} \in \mathbb{P}_k^d$ such that $\hat{p} = \hat{\mathcal{L}}_d * \hat{q}$. Finally, the bijectivity of $T^{\tau_d^i(v)}$ provides

$$ p = \hat{p} \circ (T^{\tau_d^i(v)})^{-1} = \hat{\mathcal{L}}_d \circ (T^{\tau_d^i(v)})^{-1} * \hat{q} \circ (T^{\tau_d^i(v)})^{-1} = \mathcal{L}_i^v * \hat{q} \circ (T^{\tau_d^i(v)})^{-1}, $$

and $q \stackrel{\text{def.}}{=} \hat{q} \circ (T^{\tau_d^i(v)})^{-1} \in \mathbb{P}_k^d$ answers the question. $\square$

# 9 Lag$\mathbb{P}_k^d$ Simplicial Lagrange Finite Element

The Lagrange finite elements are nodal finite elements. This means that the degrees of freedom $\Sigma$ are evaluations at specific points in the geometric element, called *nodes*. These nodes are the interpolation points where the approximate solution of the PDE
is computed. The Lagrange nodes are characterized by the fact that they are evenly distributed over the element, thus, given a spatial dimension $d > 0$, they are fully defined by choosing the vertices **v** and the approximation degree $k$, as illustrated in Figure 12. The simplicial Lagrange finite elements correspond to simplicial geometries, and to the approximation space $\mathbb{P}_k^d$ of polynomials in $d$ space variables and of (total) degree at most $k$ (see Section 8.2).

For nonzero $k$, Section 9.1 is dedicated to the Lagrange nodes and their properties, Section 9.2 to the Lagrange degrees of freedom, and Section 9.3 to the proof of unisolvence. The degenerate case $k = 0$ is treated separately in Section 9.4. And finally, the $\text{Log}\mathbb{P}_k^d$ simplicial Lagrange finite element is built in Section 9.5 for any $d$ and $k$, and Section 9.6 is devoted to the face unisolvence property.

<figure><img src="image_12.png" /><figcaption>Fig. 12 Lagrange nodes of the reference simplex for $d \in \{1, 2, 3\}$ (from top to bottom), and $k$ in $\{0, 1, 2, 3, 4\}$ (from left to right). The nodes corresponding to the highest degree are depicted in red, and the others in blue. For readability, the multi-indices are indicated in all cases when $d = 1, 2$, and only for the highest degree when $d = 3$.</figcaption></figure>

## 9.1 Simplicial Lagrange Nodes

One of the strengths of FEM is to bring all computations in mesh cells of various size and aspect ratio back to the sole case of the reference unit cell. The same operating mode can be used in the process of building finite elements such as the simplicial Lagrange FE: definitions and properties are first established in the case of the reference simplex (see Section 7.1), and then transported to the generic case by means of the bijective affine geometric transformation $T^v$ (see Section 7.2.2).

Following Section 7, reference and generic vertices are respectively denoted $\hat{v}$ and $v$. The reference vertices are always affinely independent, while this is assumed for **v** when necessary. Section 9.1.1 is dedicated to the definition of the Lagrange nodes.
The reference nodes only depend on $d$ and $k$. They are denoted $\hat{\mathbf{a}}$, and sometimes $\hat{\mathbf{a}}^d$, or $\hat{\mathbf{a}}^k$, to emphasize the dimension, or the degree of approximation. Given (generic) vertices $\mathbf{v}$, the corresponding nodes are denoted $\mathbf{a}^v$. The interaction of nodes with face hyperplanes is discussed in Section 9.1.2. Section 9.1.3 focuses on the increment in the spatial dimension, and Section 9.1.4 on the increment in the degree of approximation, as shown in Figure 12. In the latter section, reference and generic subvertices are respectively denoted $\hat{\mathbf{v}}$ and $\mathbf{v}$, and reference and generic subnodes $\hat{\mathbf{a}}$ and $\mathbf{a}^v$. Recall from Section 8.1 that multi-indices used to number nodes are denoted $\alpha$. Usually, $\alpha$ belongs to $\mathcal{A}_k^d$, and $\alpha'$ with a prime is used to denote a “smaller” multi-index, either in $\mathcal{A}_k^{d-1}$ or in $\mathcal{A}_{k-1}^d$.

In this section, the spatial dimension $d$ and the approximation degree $k$ are assumed positive, and sometimes, we need $k > 1$. However, in the Rocq formalization, all definitions are total functions, also set for $d = 0$ or $k = 0$. For instance, the inverse of noninvertible elements in rings (such as $0$ in $\mathbb{R}$) is set to $0$, but this special value is not exploited, while some properties may need additional hypotheses.

### 9.1.1 Lagrange Reference and Generic Nodes

The Lagrange reference nodes (see Figure 12) and the generic nodes associated with a family $\mathbf{v}$ of $d+1$ vertices in $\mathbb{R}^d$ are defined as

$$ \forall \alpha \in \mathcal{A}_k^d, \quad \hat{\mathbf{a}}_{\alpha} \stackrel{\text{def.}}{=} \left(\frac{1}{k}\right) \alpha, \quad \text{and} \quad \mathbf{a}_{\alpha}^{\mathbf{v}} \stackrel{\text{def.}}{=} T^{\mathbf{v}}(\hat{\mathbf{a}}_{\alpha}), \qquad (25) $$

where $(\frac{1}{k} \cdot)$ is the scaling by $\frac{1}{k}$, see Section 7.2.2. In the sequel, the latter is simplified under the compact form $\mathbf{a}^{\mathbf{v}} = T^{\mathbf{v}}(\hat{\mathbf{a}})$. In Rocq,

**Definition node_ref**[2]: $[0..\text{pbinom d k}] \to '\mathbb{R}^d := \lambda \text{idk} \Rightarrow \text{scal} (\text{/ INR k}) (\text{mapF INR (Adk d k idk)})$.

**Definition node**[2] ($\text{vtx} : [0..d] \to '\mathbb{R}^d$) : $[0..\text{pbinom d k}] \to '\mathbb{R}^d := \text{mapF (T_geometry \text{vtx}) node\_ref}$.

where INR : $\mathbb{N} \to \mathbb{R}$ is the injection of $\mathbb{N}$ into $\mathbb{R}$. The finite family operator $\text{mapF}$[2] (see Section 4.2), and scal (from the Coquelicot library, see also Section 4.3.3), allow higher-level proofs, directly on the whole family of points, just as in pen-and-paper proofs. Since the geometric transformation for the reference vertices $\hat{\mathbf{v}}$ is the identity map (see Section 7.2.2), we have $\mathbf{a}^{\mathbf{v}} = \hat{\mathbf{a}}$.

From the definition of the reference Lagrange polynomials (see Section 7.2.1) and of the geometric transformation $T^{\mathbf{v}}$ (see Section 7.2.2), we have for all $\alpha \in \mathcal{A}_k^d$,

$$ \hat{\mathcal{L}}_0(\hat{\mathbf{a}}_{\alpha}) = 1 - \frac{|\alpha|}{k}, \quad \forall i \in [1..d], \quad \hat{\mathcal{L}}_i(\hat{\mathbf{a}}_{\alpha}) = \frac{\alpha_{i-1}}{k}, \qquad (26) $$

$$ \text{and} \quad \mathbf{a}_{\alpha}^{\mathbf{v}} = \left(1 - \frac{|\alpha|}{k}\right) \mathbf{v}_0 + \sum_{i=0}^{d-1} \frac{\alpha_i}{k} \mathbf{v}_{i+1}. $$

Thus, vertices are obviously specific nodes: $\mathbf{v}_0 = \mathbf{a}_0^{\mathbf{v}}$ and $\mathbf{v}_{i+1} = \mathbf{a}_{i+1}^{\mathbf{v}}$ (see Figure 12). Moreover, when $k = 1$, the nodes are only and exactly the vertices (node\_vtx\_d1[2]).
The reference nodes $\hat{a}$ are distinct (node\_ref\_inj$^*$). Thus, there are $\text{card}(\mathcal{A}_k^d)$, i.e. $\binom{d+k}{d}$, distinct reference nodes. When the vertices $v$ are affinely independent, the weights $1-\frac{|\alpha|}{k}$ and $(\frac{\alpha_i}{k})_{i \in [0..d]}$ are the barycentric coordinates of the nodes with respect to the vertices (their sum is 1, see Section 7.2.1). Moreover, in this case, the geometric transformation is bijective, the Lagrange nodes are also distinct, and their number is also given by

**Definition** nnode\_LagPdk$^*$ (d k : $\mathbb{N}$) : $\mathbb{N}$ := (pbinom d k). + 1.

It stands as the number of Lagrange degrees of freedom in Section 9.2.

## 9.1.2 Lagrange Nodes and Faces

Some Lagrange nodes live on the (hyper)faces (i.e. the boundaries) of the simplex, and thus on the face hyperplanes that carry the faces (see Section 7.3.1 and Figure 5). These nodes can simply be characterized by their multi-index: either of sum $k$, or with a zero component.

The face hyperplanes being the kernels of the Lagrange polynomials (see Section 7.3.1), this is straightforward from (26), and the definition of $\mathcal{A}_{k,i}^d$ (see Section 8.1.4) for the reference nodes $\hat{a}$. And then, when the vertices $v$ are affinely independent, the bijective $T^v$ transports reference vertices, face hyperplanes and nodes to their generic counterparts. For all $i \in [0..d]$ and $\alpha \in \mathcal{A}_k^d$, we have

$$ \mathbf{a}_{\alpha}^{v} \in \mathcal{H}_{i}^{v} \iff T^{v}(\hat{\mathbf{a}}_{\alpha}) \in T^{v}(\hat{\mathcal{H}}_{i}) \iff \hat{\mathbf{a}}_{\alpha} \in \hat{\mathcal{H}}_{i} \iff \alpha \in \mathcal{A}_{k,i}^{d}. \quad (27) $$

This is formalized in Rocq as

**Context** {d1 k : nat}.

**Let** d := d1 + 1.

**Hypothesis** Hk : k ≠ 0.

**Context** {vtx : [0..d] → $\mathbb{R}^d$}.

**Hypothesis** Hvtx : aff\_indep\_ms vtx.

**Lemma** node\_Hface$^*$ : ∀ i (idk : [0..pbinom d k]), Hface vtx i (node vtx idk) ↔ in\_ASdki i idk.

Note that $d > 0$ is granted here structurally by d := d1+1, which is needed for the predicate in\_ASdki$^*$.

## 9.1.3 Incrementing the Spatial Dimension

We present now results used in Section 9.3 for the proof of unisolvence in the induction step incrementing the spatial dimension (and to get face unisolvence in Section 9.6).

In Section 7.3.2, the injection $\hat{I}_i : \mathbb{R}^{d-1} \to \mathbb{R}^d$ for $i \in [0..d]$ maps $\mathbb{R}^{d-1}$ onto $\hat{\mathcal{H}}_i$, either by inserting in front one minus the sum ($i = 0$), or by inserting 0 at position $i-1$ ($i \neq 0$). It also maps the reference nodes in $\mathbb{R}^{d-1}$ onto the reference nodes in $\mathbb{R}^d$ belonging to $\hat{\mathcal{H}}_i$ (see Figure 5 for the case $i = 0$),

$$ \forall \alpha' \in \mathcal{A}_{k}^{d-1}, \quad \hat{I}_{i}(\hat{\mathbf{a}}_{\alpha'}^{d-1}) = \hat{\mathbf{a}}_{\frac{d}{2}, i}^{d}(\alpha'), $$
where the $d-1$ and $d$ superscripts indicate the dimension in which the reference nodes live. In Section 8.1.4, the injection $\mathbb{I}_{k,i}^d : \mathbb{N}^{d-1} \to \mathbb{N}^d$ maps $\mathcal{A}_k^{d-1}$ onto $\mathcal{A}_{k,i}^d$, either by inserting in front $k$ minus the sum ($i=0$), or by inserting 0 at position $i-1$ ($i \neq 0$). The proof is merely a matter of commutativity with the scaling by $\frac{1}{k}$. From (25), (17), and (22), when $i=0$, we have

$$ \hat{I}_0(\hat{\mathbf{a}}_{\alpha'}^{d-1}) = \left( 1 - \sum_{i=0}^{d-1} \frac{\alpha'_i}{k}, \left( \frac{1}{k} \cdot \right)_{\alpha'} \right) = \left( \frac{1}{k} \cdot \right)_{\alpha'} (k - |\alpha'|, \alpha') = \hat{\mathbf{a}}_{\mathbb{I}_{k,0}^d(\alpha')}^d. $$

And this is even simpler when $i \neq 0$ as the insertion of 0 is linear.

Then, by applying the geometric transformation on both sides, we obtain

$$ \forall \alpha' \in \mathcal{A}_k^{d-1}, \quad I_i^{\mathbf{v}}(\hat{\mathbf{a}}_{\alpha'}^{d-1}) = \mathbf{a}_{\mathbb{I}_{k,i}^d(\alpha')}^{\mathbf{v}}, \qquad (28) $$

and in Rocq,

Context $\{d1 k : nat\}$.

Let $d := d1.+1$.

Hypothesis $Hk : k \neq 0$.

Variable $vtx : [0..d] \to \mathbb{R}^d$.

Lemma $inj\_Hface\_node^\square : \forall \{i\} (idk : [0..pbinom d1 k]), inj\_Hface vtx i (node\_ref idk) = node vtx (inj\_ASdki i idk)$.

Note that the vertices vtx do not need here to be affinely independent.

### 9.1.4 Incrementing the Degree of Approximation: Lagrange Subnodes

Subnodes are the key ingredient in Section 9.3 for the proof of unisolvence in the induction step incrementing the degree of approximation. And this works well because Lagrange nodes are evenly distributed.

What we call reference and generic “subvertices” $\hat{\mathbf{v}}$ and $\mathbf{v}$, and reference and generic “subnodes” $\hat{\mathbf{a}}$ and $\mathbf{a}$ are defined by (see Figure 13)

$$ \hat{\mathbf{v}} \stackrel{\text{def.}}{=} \left( \frac{k-1}{k} \cdot \right) \hat{\mathbf{v}}, \qquad \mathbf{v} \stackrel{\text{def.}}{=} T^{\mathbf{v}} \hat{\mathbf{v}}, \qquad (29) $$

$$ \forall \alpha' \in \mathcal{A}_{k-1}^d, \quad \hat{\mathbf{a}}_{\alpha'} \stackrel{\text{def.}}{=} \left( \frac{k-1}{k} \cdot \right) \hat{\mathbf{a}}_{\alpha'}^{k-1}, \qquad \mathbf{a}_{\alpha'} \stackrel{\text{def.}}{=} T^{\mathbf{v}} \hat{\mathbf{a}}_{\alpha'}, \qquad (30) $$

where the $k-1$ superscript indicates that these reference nodes correspond to the degree of approximation $k-1$. And in Rocq, with vertices vtx of type $[0..d] \to \mathbb{R}^d$,

Definition $sub\_vtx\_ref^\square k : [0..d] \to \mathbb{R}^d := \text{mapF}(\text{scal}(\text{INR} k.-1 / \text{INR} k)) \text{ vtx\_ref}$.

Definition $sub\_vtx^\square k vtx : [0..d] \to \mathbb{R}^d := \text{mapF}(\text{T\_geom vtx}) (\text{sub\_vtx\_ref} k)$.

Definition $sub\_node\_ref^\square : [0..pbinom d k.-1] \to \mathbb{R}^d := \text{mapF}(\text{scal}(\text{INR} k.-1 / \text{INR} k)) \text{ node\_ref}$.
**Definition sub_node** k vtx : [0..pbinom d k.−1] → 'R^d := mapF (T_geometry vtx) sub_node_ref.

<figure><img src="image_13.png" /><figcaption>Fig. 13 Passage from the reference subnodes of $LogP_3^d$ (degree $\le 3$) to the reference nodes of $LogP_2^d$ (degree $\le 2$) in the case $d = k = 3$. Plot of subvertices $(\hat{v}_i)_{i \in [0..d]}$ and subnodes $(\hat{a}_{\alpha'})_{\alpha' \in A_2^3}$.</figcaption></figure>

The nodes $(\hat{a}_{\alpha})_{\alpha \in A_3^3}$ are set in the tetrahedron defined by the vertices $(\hat{v}_0, \hat{v}_1, \hat{v}_2, \hat{v}_3)$.

The subvertices are $(\hat{v}_0, \hat{a}_{2\delta_0}, \hat{a}_{2\delta_1}, \hat{a}_{2\delta_2})$. The subnodes $(\hat{a}_{\alpha'})_{\alpha' \in A_2^3}$ are defined with respect to the tetrahedron whose vertices are $(\hat{v}_i)_{i \in [0..d]}$. These are the original nodes $(\hat{a}_{\alpha})_{\alpha \in A_3^3}$, except the (small blue) nodes corresponding to $k = 3$ (with indices in $C_3^3$).

Thus, geometrically, passing from $P_k^d$ to $P_{k-1}^d$ amounts to remove the face hyperplane $\hat{H}_0$ from the tetrahedron defined by $(\hat{v}_0, \hat{v}_1, \hat{v}_2, \hat{v}_3)$.

Consider now $k > 1$. The reference subvertices $\hat{v}$ are not reference vertices, but they are affinely independent. In the same way, when the vertices $v$ are affinely independent (and $T^v$ is bijective), the subvertices $v$ are also affinely independent. On the contrary, the reference subnodes $\hat{a}$ are some of the reference nodes, and the subnodes $\hat{a}^v$ are some of the nodes $\hat{a}^v$. Indeed, from (25) and (30), and since $k > 1$, we have for all $\alpha' \in A_{k-1}^d$,

$$ \hat{a}_{\alpha'} = \left( \frac{k-1}{k} \right) \left( \frac{1}{k-1} \right) \alpha' = \left( \frac{1}{k} \right) \alpha' = \hat{a}_{\alpha'}^k, \quad \text{and} \\ \underline{a}_{\alpha'}^v = T^v \underline{a}_{\alpha'} = T^v \hat{a}_{\alpha'}^k = \underline{a}_{\alpha'}^{v,k}. \qquad (31) $$

This is formalized in Rocq as

**Lemma sub_node_eq** : sub_node k vtx = widenF _ (node vtx).

where `widenF _` only keeps the first multi-indices in $A_{k-1}^d$ (see Section 4.2). Moreover, the subnodes $\underline{a}^v$ are also the nodes $\underline{a}^{v,k-1}$, of degree $k-1$, associated with the sub-vertices,
**Lemma sub_node_node**$^2$ : sub_node k vtx = node (sub_vtx k vtx).

Indeed, from (25), (29), (30), (15) and (16), (the latter with $\frac{k-1}{k} \neq 0$), we have

$$ 
\begin{aligned}
\forall \alpha' \in \mathcal{A}_{k-1}^d, \quad \mathbf{a}_{\alpha'}^{\mathbf{v},k-1} &= T^{T^{\mathbf{v}} \hat{\mathbf{v}} \hat{\mathbf{a}}_{\alpha'}^{k-1}} = T^{\mathbf{v}} \circ T^{\mathbf{v} \hat{\mathbf{v}} \hat{\mathbf{a}}_{\alpha'}^{k-1}} = T^{\mathbf{v}} \circ T^{\left(\frac{k-1}{k}\right) \mathbf{v} \hat{\mathbf{a}}_{\alpha'}^{k-1}} \\
&= T^{\mathbf{v}} \left(\frac{k-1}{k}\right) \cdot \hat{\mathbf{a}}_{\alpha'}^{k-1} = \underline{\mathbf{a}}_{\alpha'}^{\mathbf{v}}.
\end{aligned}
\quad (32) $$

And finally, when the vertices $\mathbf{v}$ are affinely independent, the subnodes $\underline{\mathbf{a}}^{\mathbf{v}}$ are not in the face hyperplane $\mathcal{H}_0^{\mathbf{v}}$ opposite the vertex $\mathbf{v}_0$,

**Context** $\{vtx : [0..d] \to \mathbb{R}^d\}$.

**Hypothesis** $Hvtx : aff_indep_ms vtx$.

**Lemma sub_node_out_Hface_0**$^2$ : $\forall$ idk, $\neg$ Hface $Hvtx$ ord0 (sub_node k vtx idk).

Indeed, from (27), the reference nodes $\hat{\mathbf{a}} \in \hat{\mathcal{H}}_0$ correspond to multi-indices in $C_k^d$, which is disjoint from $\mathcal{A}_{k-1}^d$, see (19). Thus, for all $\alpha' \in \mathcal{A}_{k-1}^d$, $\hat{\mathbf{a}}_{\alpha'} \notin \hat{\mathcal{H}}_0$, and from (27), we have $\underline{\mathbf{a}}_{\alpha'}^{\mathbf{v}} \notin \mathcal{H}_0^{\mathbf{v}}$.

## 9.2 Lagrange Degrees of Freedom

The nodal Lagrange degrees of freedom are the point evaluation at the Lagrange nodes defined above. They are expressed in Rocq as

**Context** $\{d : nat\}$.

**Variable k** : nat.

**Variable vtx** : $[0..d] \to \mathbb{R}^d$.

**Definition** $\text{FE\_LagPSdSk\_S\_dof}^2$ : $[0..\text{nnode\_LagPdk d k}) \to \text{FRd d} \to \mathbb{R}:= \lambda$ idk p $\Rightarrow$ p (node vtx idk).

The linearity of the Lagrange degrees of freedom ($\text{FE\_LagPSdSk\_S\_dof\_lm}^2$) simply comes from the generic result $\text{pt\_eval\_lm}^2$ for any point evaluation function.

From Section 6, the unisolvence property to establish for the construction of an instance of type $\text{FE}^2$ is that the restriction to $\mathbb{P}_k^d$ of the linear mapping $\Phi_{\Sigma^{\text{Log}_{\mathbb{P}_k^d}}} \text{ has a zero kernel.}$ Mathematically, this property is expressed as the following predicate $\Pi$, for all $d, k \in \mathbb{N}^*$ (zero values are handled otherwise),

$$ \Pi(d, k) \stackrel{\text{def.}}{=} \left[ \forall \mathbf{v} \in (\mathbb{R}^d)^{d+1} \text{ affinely independent, } \ker \Phi_{\Sigma^{\text{Log}_{\mathbb{P}_k^d}}} \cap \mathbb{P}_k^d = \{0\} \right]. \quad (33) $$

We recall that for a function $f : E \to F$ and $P \subset E$, the property $\ker f \cap P = \{0\}$ is formalized in Rocq as $\text{KerS0}^2$ P f. Moreover, it is convenient to start unisolvence proofs by applying $\text{KerS0\_gather\_equiv}^2$ to expose the internals of $\text{KerS0\_gather\_equiv}$ (see Sections 4.3.5 and 6). This means that to prove $\Pi(d, k)$, one has simply to prove that for any $p \in \mathbb{P}_k^d$ such that $\forall \alpha \in \mathcal{A}_k^d, p(\mathbf{a}_\alpha) = 0$, we have $p = 0$.
## 9.3 Unisolvence: Double Induction for Nonzero $d$ and $k$

When $d$ and $k > 0$, the unisolvence property stated in (33) is formalized in Rocq as

Let $\text{UP d k} := \forall \text{ vtx}, \text{aff\_indep\_ms vtx} \to \\ \text{KerS0 (Pdk d k) (gather (FE\_LagPSdSk\_S\_dof k vtx))}$. \\ 
\textbf{Theorem } FE\_LagPSdSk\_unisolvence\_inj[\square] : $\forall \{d k\}, d \neq 0 \to k \neq 0 \to \text{UP d k}$.

The proof is conducted as a double induction on $d$ and $k$, both starting at 1 as explained in Section 4.5. Thus, it is sufficient to prove three lemmas, first the two base cases

Lemma FE_LagP1Sk_unisolvence_inj[\square] : $\forall \{k\}, k \neq 0 \to \text{UP 1 k}$.

Lemma FE_LagPSd1_unisolvence_inj[\square] : $\forall \{d\}, d \neq 0 \to \text{UP d 1}$.

and then the (double) induction step

Lemma FE_LagPSdSk_unisolvence_inj\_ind[\square] : \\
$\forall \{d k\}, d \neq 0 \to k \neq 0 \to \text{UP d.+1 k} \to \text{UP d.k.+1} \to \text{UP d.+1 k.+1}$.

When $d = 1$ and $k \neq 0$, or when $k = 1$ and $d \neq 0$, both proofs follow directly from the fact that the components of any polynomial in $\mathbb{P}_k^1$, or in $\mathbb{P}_1^d$, on the corresponding Lagrange basis are simply given by its values at the Lagrange nodes, see lemmas *LagP1k_decomp[\square]* and *LagPd1_decomp_node[\square]*, both in Section 8.2.2. Thus, canceling the values at the nodes is sufficient to cancel the polynomial everywhere. Note that when $d = 1$, the unisolvence proof is actually established for any $k$ in the more general case of any pairwise distinct nodes (*FE_LagP1k_unisolvence_inj[\square]*), which is the case for the evenly distributed Lagrange nodes.

When $d$ and $k > 0$, let $\mathbf{v}$ be $d+2$ affinely independent vertices in $\mathbb{R}^{d+1}$, and let $p \in \mathbb{P}_{k+1}^{d+1}$ such that $p(\mathbf{a}_{\alpha}^{\mathbf{v}}) = 0$ for all $\alpha \in \mathcal{A}_{k+1}^{d+1}$ ($\mathbf{a}^{\mathbf{v}}$ are the Lagrange nodes associated with the vertices $\mathbf{v}$, see Section 9.1.1).

### Step 1: Factorization.

Let $p_0 \stackrel{\text{def.}}{=} p \circ I_0^\nu \in \mathcal{F}(\mathbb{R}^d, \mathbb{R})$, where $I_0^\nu : \mathbb{R}^d \to \mathcal{H}_0^\nu \subset \mathbb{R}^{d+1}$ is the bijective affine transformation mapping the reference simplex $\hat{K}^d$ onto the face included in $\mathcal{H}_0^\nu$, opposite vertex $\mathbf{v}_0$ (see its definition and main properties in Section 7.3.2).

Then, $p_0 \in \mathbb{P}_{k+1}^d$ (from *inj_Hface_comp[\square]*), and *inj_Hface_node[\square]* states that $I_0^\nu$ transports the reference nodes of $\hat{K}^d$ onto the nodes of $\mathcal{H}_0^\nu$, see (28). Hence, $p_0$ vanishes on those reference nodes, and applying the injectivity induction hypothesis $\Pi(d, k+1)$ provides its cancellation on the whole $\mathbb{R}^d$.

Then, the bijectivity of $I_0^\nu$ provides $p|_{\mathcal{H}_0^\nu} = p_0 \circ (I_0^\nu)^{-1} = 0$, and there exists $q$ in $\mathbb{P}_k^{d+1}$ such that $p = \mathcal{L}_0^\nu * q$ (from *factor_zero_on_Hface[\square]*, see Section 8.2.3).

Note that the vanishing of $p$ on the face hyperplane $\mathcal{H}_0^\nu$ is actually established at no extra cost for any face hyperplane $\mathcal{H}_i^\nu$ for $i \in [0..d]$ (*unisolvence_inj_Hface[\square]*), since it is needed in Section 9.6 for face unisolvence.
### Step 2: Cancellation.

By hypothesis, $p = \mathcal{L}_0^v * q$ vanishes on all Lagrange nodes $a_\alpha^v$, i.e. for all $\alpha \in \mathcal{A}_{k+1}^{d+1}$. Hence, $q$ vanishes on all nodes not in $\ker \mathcal{L}_0^v = \mathcal{H}_0^v$, thus corresponding to the multi-indices in $\mathcal{A}_k^{d+1} (= \mathcal{A}_{k+1}^{d+1} \setminus \mathcal{C}_{k+1}^{d+1})$, see (27), (21), and (19). This means that $q$ vanishes on the nodes $a_{\alpha'}^{v,k+1}$ for all $\alpha' \in \mathcal{A}_k^{d+1}$.

The subvertices $\underline{v}$, and the associated subnodes $(\underline{a}_{\alpha'}^v)_{\alpha' \in \mathcal{A}_k^{d+1}}$ defined in Section 9.1.4 are precisely designed to allow the application of the injectivity induction hypothesis $\Pi(d+1, k)$ to $q$. Indeed, these subnodes $\underline{a}_{\alpha'}^v$ are in the same time the above-mentioned nodes $a_{\alpha'}^{v,k+1}$ (where $q \in \mathbb{P}_k^{d+1}$ vanishes), see (31), and the nodes $a_{\alpha'}^v$ of degree $k$ for the subvertices $\underline{v}$ (for which $\Pi(d+1, k)$ can be applied), see (32). And thus, we obtain the cancellation of $q$, and of $p$. $\square$

As the induction step is done, we have proved that $\Pi$ (that is $\text{UP}$ in Rocq) holds for nonzero $d$ and $k$. This hard proof is required to build the Lagrange FE in Section 9.5.

## 9.4 Degenerate Case $k = 0$

The expression and proof path of unisolvence being different in the case $k = 0$, it is treated separately. When $k = 0$, the approximation polynomials are constant, thus any point in the geometric element could serve as the single node. The isobarycenter of the vertices is chosen for symmetry reasons, and because it is not on the boundary, thus not interacting with neighbor elements when considering a mesh. Thus, vertices are no longer nodes, and the affine independence hypothesis on the vertices is no longer needed for the unisolvence property.

The nodes and degrees of freedom are formalized as

**Definition node_d0** (vtx : $[0..d] \to \mathbb{R}^d$) : $\mathbb{R}^d := \text{isobarycenter_ms}$ vtx.

**Definition FE_LagPd0_S_dof** (vtx : $[0..d] \to \mathbb{R}^d$) : $[0..\text{nnode_LagPdk d 0}) \to \text{FRd d} \to \mathbb{R} := \lambda \_p \Rightarrow p(\text{node_d0 vtx})$.

And the predicate $\Pi$ to grant unisolvence simplifies into

$$ \forall d > 0, \quad \Pi(d, 0) \stackrel{\text{def.}}{=} [\ker \Phi_{\Sigma}^{\text{Log}_{\mathbb{P}_0^d}} \cap \mathbb{P}_0^d = \{0\}] . $$

In Rocq,

Let $\text{UP0 d} := \forall \text{vtx}, \text{KerS0 (Pdk d 0)} (\text{gather (FE_LagPd0_S_dof vtx)})$.

**Lemma FE_LagPd0_unisolvence_inj** : $\forall \{d\}, \text{UP0 d}$.

In this case, the approximation space $\mathbb{P}_0^d$ is the one-dimensional space of constant functions. The proof follows directly from the fact that a constant function vanishes iff it is zero at some point.

## 9.5 Construction of the $\text{Log}_{\mathbb{P}_k^d}$ Simplicial Lagrange Finite Element

Once the unisolvence properties are established, we can build the records corresponding to Lagrange finite elements $\text{Log}_{\mathbb{P}_k^d}$. First, for $d, k > 0$,
Context {$d : k : nat$}.

Hypothesis $Hd : d \neq 0$.

Hypothesis $Hk : k \neq 0$.

Context {$vtx : [0..d] \to 'R^d$}.

Hypothesis $Hvtx : aff_indep_ms vtx$.

Definition $FE_LagPSdSk^2 : FE d :=$
$$ mk FE Simplex vtx (nnode_LagPdk d k) (Pdk_has_dim d k) \\ (FE_LagPSdSk_S_dof_lmk vtx) (FE_LagPSdSk_unisolvence_inj Hd Hk \_ Hvtx). $$

This is the main construction of $\text{Lag}_{P(k)}^d$, based on the Lagrange degrees of freedom defined in Section 9.2 and the unisolvence proved in Section 9.3.

Then, for any $d$ and for $k=0$ from Section 9.4, where the vertices no longer need to be affinely independent,

Context {$d : nat$}.

Variable $vtx : [0..d] \to 'R^d$.

Definition $FE_LagPd0^2 : FE d :=$
$$ mk FE Simplex vtx (nnode_LagPdk d 0) (Pdk_has_dim d 0) \\ (FE_LagPd0_S_dof_lmk vtx) (FE_LagPd0_unisolvence_inj vtx). $$

However, we have to keep in mind that nondegenerate simplices will need affinely independent vertices.

And for $d=0$ from Section 6.2, where the only existing vertex does not need to be specified and $k$ is no longer relevant,

Definition $FE_LagP0k^2 : FE 0 := FE_d 0 Simplex one_not_zero_R$.

Finally, for the case $d=1$ for any $k$, the nodes need only be pairwise distinct, and need not be evenly distributed,

Context {$k : nat$}.

Context {$node_1k : [0..nnode_LagPdk 1 k) \to 'R^1$}.

Hypothesis $Hnode : injective node_1k$.

Variable $vtx : 'R^2 1$.

Definition $FE_LagP1k^2 : FE 1 :=$
$$ mk FE Simplex vtx (nnode_LagPdk 1 k) (Pdk_has_dim 1 k) \\ (FE_LagP1k_S_dof_lmk node_1k) (FE_LagP1k_unisolvence_inj Hnode). $$

This construction holds for any couple of vertices (in $\mathbb{R}$). However, we note that non-degenerate geometric elements will need affinely independent vertices (i.e. distinct). Moreover, in practice, it will be wise to assume that the nodes are located between vertices and include them. This will be useful for the construction of Lagrange FE on cuboids, see Section 10.4.

Note also that the arguments to build the Lagrange FE are not the same in all cases: either none, vertices, or nodes, with some property or not. Consequently, the unification of all definitions into a single function proves to be an attractive misconception.
Another seemingly nice idea would have been to construct only the reference Lagrange FE $\mathbb{P}_k^d$ and build $\mathbb{P}_k^d$ using a geometric transformation of Section 7.4. This proof path would have failed as the second step of the unisolvence proof applies the induction hypothesis to subvertices, which are not reference vertices. Said otherwise, we need the induction hypothesis on generic vertices to be applicable.

However, we can check for any $d, k > 0$ and independent vertices $v$ that the transformation induced by $T^v$ of the Lagrange FE specialized to the reference vertices gives the corresponding Lagrange FE back, see Figure 7, i.e. in Rocq (where the annotation allows access to all arguments, including implicit ones, here $d$)

**Context** $\{d, k : \text{nat}\}$.

**Hypothesis** $Hd : d \neq 0$.

**Hypothesis** $Hk : k \neq 0$.

**Definition** $\text{FE}_{\text{LagPSdSk\_ref}}^\square := \text{FE}_{\text{LagPSdSk}} Hd Hk (\text{@vtx\_ref\_aff\_indep d})$.

**Lemma** $\text{FE}_{\text{LagPSdSk\_eq\_from\_ref}}^\square : \forall \{vtx : [0..d] \to \mathbb{R}^d\} (Hvtx : \text{aff\_indep\_ms vtx})$, $\text{FE}_{\text{LagPSdSk}} Hd Hk Hvtx = \text{TFE}_{vtx\_from\_ref} \text{FE}_{\text{LagPSdSk\_ref}} Hvtx$.

where the transformed FE $\text{TFE}_{vtx\_from\_ref}$ is presented in Table 8. The proof directly follows from the extensionality result $\text{FE}_{\text{ext}}^\square$ for $\text{FE}^\square$ records, which stipulates that it is sufficient to check that the nonproof fields match, namely vertices, approximation spaces and degrees of freedom. This is quite straightforward since the geometric transformation $T^v$ is affine and bijective, thus it preserves barycenters (such as vertices and nodes), and the polynomial spaces $\mathbb{P}_k^d$. The slight difficulty is the handling of the casts due to different expressions for the same dimension when they are not structurally equal. The proof is actually the same when $k = 0$ ($\text{FE}_{\text{LagPd0\_eq\_from\_ref}}^\square$), and from the uniqueness result $\text{FE}_{d\_0\_uniq}^\square$ of Section 6.2 (with the scaling factor 1), it is trivial when $d = 0$.

## 9.6 Face Unisolvence

Assume that $d > 1$ and $k > 0$. The face unisolvence property states that

$$ \forall i \in [0..d], \forall p \in \mathbb{P}_k^d, (\forall \alpha \in \mathcal{A}_k^d, \mathbf{a}_\alpha \in \mathcal{H}_i^v \implies p(\mathbf{a}_\alpha) = 0) \iff p|_{\mathcal{H}_i} = 0. $$

This property states that a polynomial in $\mathbb{P}_k^d$ vanishes on a face hyperplane if and only if it is zero at the Lagrange nodes that lie in this face hyperplane. This is useful to enforce inter-element continuity using discrete conditions. Indeed, let $f$ be a function that is in $\mathbb{P}_k^d$ in two adjacent mesh cells $K_1$ and $K_2$, sharing the (nonempty) face $\mathcal{F} = K_1 \cap K_2$. Assume that face unisolvence holds on the hyperplane containing $\mathcal{F}$. Thus, by linearity of the evaluation of functions at some point, in order for $f$ to be continuous, it suffices to enforce that the two polynomials $f|_{K_1}$ and $f|_{K_2}$ take the same values at the Lagrange nodes in the face $\mathcal{F}$.

Face unisolvence holds for $\mathbb{P}_k^d$, thanks to the fact that the Lagrange nodes are evenly distributed on each face. As we do not need here to structurally grant the bounds on $d$ and $k$, this is formalized as