Merge pull request #1104 from Antoinemarteau/moment-based-reffes

Moment based reffes - Implementation of scalar Bernstein basis in barycentric coordinates.

Author: JordiManyer

NEWS.md

@@ -7,16 +7,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
-- Documentation and refactoring of Gridap.Polynomials. Since PR[#TODO](https://github.com/gridap/Gridap.jl/pull/#TODO).
+### Added
+
+- Added AMR-related methods `mark` and `estimate` to `Adaptivity` module. Implemented Dorfler marking strategy. Since PR[#1063](https://github.com/gridap/Gridap.jl/pull/1063).
+
+- Documentation and refactoring of Gridap.Polynomials. Since PR[#1072](https://github.com/gridap/Gridap.jl/pull/#1072).
 - Two new families of polynomial bases in addition to `Monomial`, `Legendre` (former `Jacobi`) and `ModalC0`: `Chebyshev` and `Bernstein`
-- `MonomialBasis` and `Q[Curl]GradMonomialBasis` have been generalized to `Legendre`, `Chebyshev` and `Bernstein` using the new `UniformPolyBasis` and `CompWiseTensorPolyBasis` respectively.
+- `MonomialBasis` and `Q[Curl]GradMonomialBasis` have been generalized to `Legendre`, `Chebyshev` and `Bernstein` using the new `CartProdPolyBasis` and `CompWiseTensorPolyBasis` respectively.
 - `PCurlGradMonomialBasis` has been generalized to `Legendre` and `Chebyshev` using the new `RaviartThomasPolyBasis`.
-- New aliases and high level constructor for `UniformPolyBasis` (former MonomialBasis): `MonomialBasis`, `LegendreBasis`, `ChebyshevBasis` and `BernsteinBasis`.
+- New aliases and high level constructor for `CartProdPolyBasis` (former MonomialBasis): `MonomialBasis`, `LegendreBasis`, `ChebyshevBasis` and `BernsteinBasis`.
 - New high level constructors for Nedelec and Raviart-Thomas polynomial bases:
     - Nedelec on simplex `PGradBasis(PT<:Polynomial, Val(D), order)`
     - Nedelec on n-cubes `QGradBasis(PT<:Polynomial, Val(D), order)`
     - Raviart on simplex `PCurlGradBasis(PT<:Polynomial, Val(D), order)`
     - Raviart on n-cubes `QCurlGradBasis(PT<:Polynomial, Val(D), order)`
+- Added `BernsteinBasisOnSimplex` that implements Bernstein polynomials in barycentric coordinates, since PR[#1104](https://github.com/gridap/Gridap.jl/pull/#1104).
 
 ## [0.18.10] - 2025-03-04
 
@@ -50,7 +55,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Existing Jacobi polynomial bases/spaces were renamed to Legendre (which they were).
 - `Monomial` is now subtype of the new abstract type`Polynomial <: Field`
-- `MonomialBasis` is now an alias for `UniformPolyBasis{...,Monomial}`
+- `MonomialBasis` is now an alias for `CartProdPolyBasis{...,Monomial}`
 - All polynomial bases are now subtypes of the new abstract type `PolynomialBasis <: AbstractVector{<:Polynomial}`
 - `get_order(b::(Q/P)[Curl]Grad...)`, now returns the order of the basis, +1 than the order parameter passed to the constructor.
 - `NedelecPreBasisOnSimplex` is renamed `NedelecPolyBasisOnSimplex`

docs/generate_diagrams.jl

@@ -53,29 +53,31 @@ end
   a1 <|-left- PolynomialBasis
 
   together {
-    struct UniformPolyBasis {
+    struct CartProdPolyBasis {
       +get_exponents
       +get_orders
     }
     struct CompWiseTensorPolyBasis
     struct RaviartThomasPolyBasis
     struct NedelecPolyBasisOnSimplex
+    struct BernsteinBasisOnSimplex
     struct ModalC0Basis {
       +get_orders
     }
   }
 
-  PolynomialBasis <|-- UniformPolyBasis
+  PolynomialBasis <|-- CartProdPolyBasis
   PolynomialBasis <|-- CompWiseTensorPolyBasis
   PolynomialBasis <|-- RaviartThomasPolyBasis
   PolynomialBasis <|-- NedelecPolyBasisOnSimplex
+  PolynomialBasis <|-- BernsteinBasisOnSimplex
   PolynomialBasis <|-- ModalC0Basis
 
   object "(<:Polynomial)Basis" as m1
   object "QGrad[<:Polynomial]Basis\nQCurlGrad[<:Polynomial]Basis" as m2
   object "PCurlGrad[<:Polynomial]Basis" as m4
   object "PGradMonomialBasis" as m5
-  UniformPolyBasis <-down- m1
+  CartProdPolyBasis <-down- m1
   CompWiseTensorPolyBasis <-down- m2
   RaviartThomasPolyBasis <-down- m4
   NedelecPolyBasisOnSimplex <-down- m5

docs/make.jl

@@ -24,6 +24,7 @@ pages = [
   "Developper notes" => Any[
     "dev-notes/block-assemblers.md",
     "dev-notes/pullbacks.md",
+    "dev-notes/bernstein.md",
     "dev-notes/autodiff.md",
   ],
 ]

docs/src/Polynomials.md

@@ -31,10 +31,10 @@ polynomial spaces. The order ``K`` 1D monomial is
 ```math
 x \rightarrow x^K,
 ```
-and the order ``\boldsymbol{K}=(K_1, K_2, \dots, K_D)`` D-dimensional monomial is defined by
+and the order ``\bm{K}=(K_1, K_2, \dots, K_D)`` D-dimensional monomial is defined by
 ```math
-\boldsymbol{x} = (x_1, x_2, \dots, x_D) \longrightarrow
-\boldsymbol{x}^{\boldsymbol{K}} = x_1^{K_1}x_2^{K_2}...x_D^{K_D} = \Pi_{i=1}^D
+\bm{x} = (x_1, x_2, \dots, x_D) \longrightarrow
+\bm{x}^{\bm{K}} = x_1^{K_1}x_2^{K_2}...x_D^{K_D} = \Pi_{i=1}^D
 x_i^{K_i}.
 ```
 
@@ -52,7 +52,7 @@ This module implements the normalized shifted [`Legendre`](@ref) polynomials,
 shifted to be orthogonal on ``[0,1]`` using the change of variable ``x
 \rightarrow 2x-1``, leading to
 ```math
-P^*_{n}(x)=\frac{1}{\sqrt{2n+1}}P_n(2x-1)=\frac{1}{\sqrt{2n+1}}(-1)^{n}\sum _{i=0}^{n}{\binom{n}{i}}{\binom{n+i}{i}}(-x)^{i}.
+P^*_{n}(x)=\frac{1}{\sqrt{2n+1}}P_n(2x-1)=\frac{1}{\sqrt{2n+1}}(-1)^{n}∑ _{i=0}^{n}{\binom{n}{i}}{\binom{n+i}{i}}(-x)^{i}.
 ```
 
 #### Chebyshev polynomials
@@ -65,9 +65,9 @@ polynomials ``U_n`` can be recursively defined by
 ```
 or explicitly defined by
 ```math
-T_{n}(x)=\sum _{i=0}^{\left\lfloor {n}/{2}\right\rfloor }{\binom
+T_{n}(x)=∑ _{i=0}^{\left\lfloor {n}/{2}\right\rfloor }{\binom
     {n}{2i}}\left(x^{2}-1\right)^{i}x^{n-2i},\qquad
-U_{n}(x)=\sum _{i=0}^{\left\lfloor {n}/{2}\right\rfloor }{\binom
+U_{n}(x)=∑ _{i=0}^{\left\lfloor {n}/{2}\right\rfloor }{\binom
     {n+1}{2i+1}}\left(x^{2}-1\right)^{i}x^{n-2i},
 ```
 where ``\left\lfloor {n}/2\right\rfloor`` is `floor(n/2)`.
@@ -82,12 +82,28 @@ The analog second kind shifted Chebyshev polynomials can be implemented by
 
 #### Bernstein polynomials
 
-The [`Bernstein`](@ref) polynomials forming a basis of ``\mathbb{P}_K`` are
-defined by
+The univariate [`Bernstein`](@ref) polynomials forming a basis of ``ℙ_K``
+are defined by
+```math
+B^K_{n}(x) = \binom{K}{n} x^n (1-x)^{K-n}\qquad\text{ for } 0 ≤ n ≤ K.
+```
+
+The ``D``-multivariate Bernstein polynomials of degree ``K`` are defined by
 ```math
-B^K_{n}(x) = \binom{K}{n} x^n (1-x)^{K-n}\qquad\text{ for } 0\leq n\leq K.
+B^{D,K}_α(\bm{x}) = \binom{K}{α} λ(\bm{x})^α\qquad\text{for all }α ∈\mathcal{I}_K^D
 ```
-They are positive on ``[0,1]`` and sum to ``1``.
+where
+- ``\mathcal{I}^D_K=\{ α ∈ ℤ_+^{N} \quad\big|\quad |α| = K \}`` where ``N = D+1``,
+- ``\binom{|α|}{α} = \frac{|α|!}{α!} = \frac{K!}{α_1 !α_2 !\dotsα_N!}`` and ``K=|α|=∑_{1 ≤ i ≤ N} α_i``,
+and ``λ(\bm x)`` is the set of barycentric coordinates of ``\bm x`` relative to a given simplex,
+see the developer notes on the [Bernstein bases algorithms](@ref).
+
+The superscript ``D`` and ``K`` in ``B^{D,K}_α(x)`` can be omitted because they
+are always determined by ``α`` using ``{D=\#(α)-1}`` and ``K=|α|``. The set
+``\{B_α\}_{α∈\mathcal{I}_K^D}`` is a basis of ``ℙ^D_K``, implemented by
+[`BernsteinBasisOnSimplex`](@ref).
+
+The Bernstein polynomials sum to ``1``, and are positive in the simplex.
 
 #### ModalC0 polynomials
 
@@ -97,7 +113,7 @@ polynomials ``\phi_K`` for which ``\phi_K(0) = \delta_{0K}`` and ``\phi_K(1) =
 17 of [Badia-Neiva-Verdugo 2022](https://doi.org/10.1016/j.camwa.2022.09.027).
 
 When `ModalC0` is used as a `<:Polynomial` parameter in
-[`UniformPolyBasis`](@ref) or other bases (except `ModalC0Basis`), the trivial
+[`CartProdPolyBasis`](@ref) or other bases (except `ModalC0Basis`), the trivial
 bounding box `(a=Point{D}(0...), b=Point{D}(1...))` is assumed, which
 coincides with the usual definition of the ModalC0 bases.
 
@@ -106,90 +122,90 @@ coincides with the usual definition of the ModalC0 bases.
 
 #### P and Q spaces
 
-Let us denote ``\mathbb{P}_K(x)`` the space of univariate polynomials of order up to ``K`` in the varible ``x``
+Let us denote ``ℙ_K(x)`` the space of univariate polynomials of order up to ``K`` in the varible ``x``
 ```math
-\mathbb{P}_K(x) = \text{Span}\big\{\quad x\rightarrow x^i \quad\big|\quad 0\leq i\leq K \quad\big\}.
+ℙ_K(x) = \text{Span}\big\{\quad x\rightarrow x^i \quad\big|\quad 0 ≤ i ≤ K \quad\big\}.
 ```
 
-Then, ``\mathbb{Q}^D`` and ``\mathbb{P}^D`` are the spaces for Lagrange elements
+Then, ``ℚ^D`` and ``ℙ^D`` are the spaces for Lagrange elements
 on D-cubes and D-simplices respectively, defined by
 ```math
-\mathbb{Q}^D_K = \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^\alpha \quad\big|\quad 0\leq
-    \alpha_1, \alpha_2, \dots, \alpha_D \leq K \quad\big\},
+ℚ^D_K = \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^α \quad\big|\quad 0 ≤
+    α_1, α_2, \dots, α_D  ≤ K \quad\big\},
 ```
 and
 ```math
-\mathbb{P}^D_K = \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^\alpha \quad\big|\quad 0\leq
-    \alpha_1, \alpha_2, \dots, \alpha_D \leq K;\quad \sum_{d=1}^D \alpha_d \leq
+ℙ^D_K = \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^α \quad\big|\quad 0 ≤
+    α_1, α_2, \dots, α_D  ≤ K;\quad ∑_{d=1}^D α_d  ≤
     K \quad\big\}.
 ```
 
-To note, there is ``\mathbb{P}_K = \mathbb{P}^1_K = \mathbb{Q}^1_K``.
+To note, there is ``ℙ_K = ℙ^1_K = ℚ^1_K``.
 
 #### Serendipity space Sr
 
 The serendipity space, commonly used for serendipity finite elements on n-cubes,
 are defined by
 ```math
-\mathbb{S}r^D_K = \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^\alpha \quad\big|\quad 0\leq
-    \alpha_1, \alpha_2, \dots, \alpha_D \leq K;\quad
-    \sum_{d=1}^D \alpha_d\;\mathbb{1}_{[2,K]}(\alpha_d) \leq K \quad\big\}
+\mathbb{S}r^D_K = \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^α \quad\big|\quad 0 ≤
+    α_1, α_2, \dots, α_D  ≤ K;\quad
+    ∑_{d=1}^D α_d\;\mathbb{1}_{[2,K]}(α_d)  ≤ K \quad\big\}
 ```
-where ``\mathbb{1}_{[2,K]}(\alpha_d)`` is ``1`` if ``\alpha_d\geq 2`` or else
+where ``\mathbb{1}_{[2,K]}(α_d)`` is ``1`` if ``α_d ≥ 2`` or else
 ``0``.
 
 #### Homogeneous P and Q spaces
 
 It will later be useful to define the homogeneous Q spaces
 ```math
-\tilde{\mathbb{Q}}^D_K = \mathbb{Q}^D_K\backslash\mathbb{Q}^D_{K-1} =
-    \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^\alpha \quad\big|\quad 0\leq \alpha_1,
-    \alpha_2, \dots, \alpha_D \leq K; \quad \text{max}(\alpha) = K \quad\big\},
+\tilde{ℚ}^D_K = ℚ^D_K\backslashℚ^D_{K-1} =
+    \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^α \quad\big|\quad 0 ≤ α_1,
+    α_2, \dots, α_D  ≤ K; \quad \text{max}(α) = K \quad\big\},
 ```
 and homogeneous P spaces
 ```math
-\tilde{\mathbb{P}}^D_K = \mathbb{P}^D_K\backslash \mathbb{P}^D_{K-1} =
-    \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^\alpha \quad\big|\quad 0\leq \alpha_1,
-    \alpha_2, \dots, \alpha_D \leq K;\quad \sum_{d=1}^D \alpha_d = K \quad\big\}.
+\tilde{ℙ}^D_K = ℙ^D_K\backslash ℙ^D_{K-1} =
+    \text{Span}\big\{\quad \bm{x}\rightarrow\bm{x}^α \quad\big|\quad 0 ≤ α_1,
+    α_2, \dots, α_D  ≤ K;\quad ∑_{d=1}^D α_d = K \quad\big\}.
 ```
 
 
-#### Nédélec spaces
+#### Nédélec ``curl``-conforming spaces
 
 The Kᵗʰ Nédélec polynomial spaces on respectively rectangles and
 triangles are defined by
 ```math
-\mathbb{ND}^2_K(\square) = \left(\mathbb{Q}^2_K\right)^2 \oplus
-    \left(\begin{array}{c} y^{K+1}\,\mathbb{P}_K(x)\\ x^{K+1}\,\mathbb{P}_K(y) \end{array}\right)
+ℕ𝔻^2_K(\square) = \left(ℚ^2_K\right)^2 ⊕
+    \left(\begin{array}{c} y^{K+1}\,ℙ_K(x)\\ x^{K+1}\,ℙ_K(y) \end{array}\right)
 ,\qquad
-\mathbb{ND}^2_K(\bigtriangleup) =\left(\mathbb{P}^2_K\right)^2 \oplus\bm{x}\times(\tilde{\mathbb{P}}^2_K)^2,
+ℕ𝔻^2_K(\bigtriangleup) =\left(ℙ^2_K\right)^2 ⊕\bm{x}\times(\tilde{ℙ}^2_K)^2,
 ```
 where ``\times`` here means ``\left(\begin{array}{c} x\\ y
 \end{array}\right)\times\left(\begin{array}{c} p(\bm{x})\\ q(\bm{x})
 \end{array}\right) = \left(\begin{array}{c} y p(\bm{x})\\ -x q(\bm{x})
-\end{array}\right)`` and ``\oplus`` is the direct sum of vector spaces.
+\end{array}\right)`` and ``⊕`` is the direct sum of vector spaces.
 
 Then, the Kᵗʰ Nédélec polynomial spaces on respectively hexahedra and
 tetrahedra are defined by
 ```math
-\mathbb{ND}^3_K(\square) = \left(\mathbb{Q}^3_K\right)^3 \oplus \bm{x}\times(\tilde{\mathbb{Q}}^3_K)^3,\qquad
-\mathbb{ND}^3_K(\bigtriangleup) =\left(\mathbb{P}^3_K\right)^3 \oplus \bm{x}\times(\tilde{\mathbb{P}}^3_K)^3.
+ℕ𝔻^3_K(\square) = \left(ℚ^3_K\right)^3 ⊕ \bm{x}\times(\tilde{ℚ}^3_K)^3,\qquad
+ℕ𝔻^3_K(\bigtriangleup) =\left(ℙ^3_K\right)^3 ⊕ \bm{x}\times(\tilde{ℙ}^3_K)^3.
 ```
 
-``\mathbb{ND}^D_K(\square)`` and ``\mathbb{ND}^D_K(\bigtriangleup)`` are of
-order K+1 and the curl of their elements are in ``(\mathbb{Q}^D_K)^D``
-and ``(\mathbb{P}^D_K)^D`` respectively.
+``ℕ𝔻^D_K(\square)`` and ``ℕ𝔻^D_K(\bigtriangleup)`` are of
+order K+1 and the curl of their elements are in ``(ℚ^D_K)^D``
+and ``(ℙ^D_K)^D`` respectively.
 
-#### Raviart-Thomas spaces
+#### Raviart-Thomas and Nédélec ``div``-conforming Spaces
 
 The Kᵗʰ Raviart-Thomas polynomial spaces on respectively D-cubes and
 D-simplices are defined by
 ```math
-\mathbb{ND}^D_K(\square) = \left(\mathbb{Q}^D_K\right)^D \oplus \bm{x}\;\tilde{\mathbb{Q}}^D_K, \qquad
-\mathbb{ND}^D_K(\bigtriangleup) = \left(\mathbb{P}^D_K\right)^D \oplus \bm{x}\;\tilde{\mathbb{P}}^D_K,
+ℕ𝔻^D_K(\square) = \left(ℚ^D_K\right)^D ⊕ \bm{x}\;\tilde{ℚ}^D_K, \qquad
+ℕ𝔻^D_K(\bigtriangleup) = \left(ℙ^D_K\right)^D ⊕ \bm{x}\;\tilde{ℙ}^D_K,
 ```
 these bases are of dimension K+1 and the divergence of their elements are in
-``\mathbb{Q}^D_K`` and ``\mathbb{P}^D_K`` respectively.
+``ℚ^D_K`` and ``ℙ^D_K`` respectively.
 
 
 #### Filter functions
@@ -205,20 +221,18 @@ The signature of the filter functions should be
     (term,order) -> Bool
 
 where `term` is a tuple of `D` integers containing the exponents of a
-multivariate monomial, that correspond to the multi-index ``\alpha`` previously
+multivariate monomial, that correspond to the multi-index ``α`` previously
 used in the P/Q spaces definitions.
 
-When using [`hierarchical`](@ref isHierarchical) 1D bases, the following
-filters can be used to define associated polynomial spaces:
-
-| space       | filter                                                           |
-| :-----------| :--------------------------------------------------------------- |
-| ℚᴰ          | `_q_filter(e,order) = maximum(e) <= order`                       |
-| ℚᴰₙ\\ℚᴰₙ₋₁  | `_qs_filter(e,order) = maximum(e) == order`                      |
-| ℙᴰ          | `_p_filter(e,order) = sum(e) <= order`                           |
-| ℙᴰₙ\\ℙᴰₙ₋₁  | `_ps_filter(e,order) = sum(e) == order`                          |
-| 𝕊rᴰₙ        | `_ser_filter(e,order) = sum( [ i for i in e if i>1 ] ) <= order` |
+The following example filters can be used to define associated polynomial spaces:
 
+| space       | filter                                                       | possible family                       |
+| :-----------| :------------------------------------------------------------| :------------------------------------ |
+| ℚᴰ          | `_q_filter(e,order) = maximum(e) <= order`                   | All                                   |
+| ℚᴰₙ\\ℚᴰₙ₋₁  | `_qh_filter(e,order) = maximum(e) == order`                  | [`Monomial`](@ref)                    |
+| ℙᴰ          | `_p_filter(e,order) = sum(e) <= order`                       | All                                   |
+| ℙᴰₙ\\ℙᴰₙ₋₁  | `_ph_filter(e,order) = sum(e) == order`                      | [`Monomial`](@ref)                    |
+| 𝕊rᴰₙ        | `_ser_filter(e,order) = sum( i for i in e if i>1 ) <= order` | [`hierarchical`](@ref isHierarchical) |
 
 ## Types for polynomial families
 
@@ -252,7 +266,7 @@ ModalC0
 
 ```@docs
 PolynomialBasis
-get_order
+get_order(::PolynomialBasis)
 MonomialBasis(args...)
 MonomialBasis
 LegendreBasis(args...)
@@ -263,12 +277,16 @@ BernsteinBasis(args...)
 ```
 !!! warning
     Calling `BernsteinBasis` with the filters (e.g. a `_p_filter`) rarely
-    yields a basis for the associated space (e.g. ``\mathbb{P}``).  Indeed, the
+    yields a basis for the associated space (e.g. ``ℙ``).  Indeed, the
     term numbers do not correspond to the degree of the polynomial, because the
     basis is not [`hierarchical`](@ref isHierarchical).
 
 ```@docs
 BernsteinBasis
+BernsteinBasisOnSimplex
+BernsteinBasisOnSimplex(::Val,::Type,::Int,vertices=nothing)
+bernstein_terms
+bernstein_term_id
 PGradBasis
 QGradBasis
 PCurlGradBasis
@@ -279,10 +297,10 @@ QCurlGradBasis
 ![](./assets/poly_2.svg)
 
 ```@docs
-UniformPolyBasis
-UniformPolyBasis(::Type, ::Val{D}, ::Type, ::Int, ::Function) where D
-UniformPolyBasis(::Type, ::Val{D}, ::Type{V}, ::NTuple{D,Int}, ::Function) where {D,V}
-get_orders
+CartProdPolyBasis
+CartProdPolyBasis(::Type, ::Val{D}, ::Type, ::Int, ::Function) where D
+CartProdPolyBasis(::Type, ::Val{D}, ::Type{V}, ::NTuple{D,Int}, ::Function) where {D,V}
+get_orders(::CartProdPolyBasis)
 get_exponents
 CompWiseTensorPolyBasis
 NedelecPolyBasisOnSimplex