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

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

docs/src/dev-notes/bernstein.md

@@ -0,0 +1,144 @@
+# Bernstein bases algorithms
+
+### Barycentric coordinates
+
+A ``D``-dimensional simplex ``T`` is defined by ``N=D+1`` vertices ``\{v_1,
+v_2, …, v_N\}=\{v_i\}_{i∈1:N}``. The barycentric coordinates
+``λ(\bm{x})=\{λ_j(\bm{x})\}_{1 ≤ j ≤ N}`` are uniquely
+defined by:
+```math
+\bm{x} = ∑_{1 ≤ j ≤ N} λ_j(\bm{x})v_j \quad\text{and}\quad
+∑_{1≤ j≤ N} λ_j(\bm{x}) = 1,
+```
+as long as the simplex is non-degenerate (vertices are not all in one
+hyperplane).
+
+Assuming the simplex polytopal (has flat faces), this change of coordinates is
+affine, and is implemented using:
+```math
+λ(\bm{x}) = M\left(\begin{array}{c} 1\\ x_1\\ ⋮\\ x_D \end{array}\right)
+\quad\text{with}\quad M =
+\left(\begin{array}{cccc}
+1 & 1 & ⋯ & 1 \\
+(v_1)_1 & (v_2)_1  & ⋯ & (v_N)_1  \\
+⋮  & ⋮   & ⋯ & ⋮   \\
+(v_1)_D & (v_2)_D  & ⋯ & (v_N)_D  \\
+\end{array}\right)^{-1}
+```
+where the inverse exists because ``T`` is non-degenerate [1], cf. functions
+`_cart_to_bary` and `_compute_cart_to_bary_matrix`. Additionally, we have
+``∂_{x_i} λ_j(\bm{x}) = M_{j,i+1}``, so
+```math
+∇ λ_j = M_{2:N, j}.
+```
+The matrix ``M`` is all we need that depends on ``T`` in order to compute
+Bernstein polynomials and their derivatives, it is stored in the field
+`cart_to_bary_matrix` of [`BernsteinBasisOnSimplex`](@ref), when ``T`` is not
+the reference simplex.
+
+On the reference simplex defined by the vertices `get_vertex_coordinates(SEGMENT / TRI / TET⋯)`:
+```math
+\begin{aligned}
+v_1     & = (0\ 0\ ⋯\ 0), \\
+v_2     & = (1\ 0\ ⋯\ 0), \\
+⋮  &         \\
+v_N     & = (0\ ⋯\ 0\ 1),
+\end{aligned}
+```
+the matrix ``M`` is not stored because
+```math
+λ(\bm{x}) = \Big(1-∑_{1≤ i≤ D} x_i, x_1, x_2, ⋯, x_D\Big)
+\quad\text{and}\quad
+∂_{x_i} λ_j = δ_{i+1,j} - δ_{1j} = M_{j,i+1}.
+```
+
+### Bernstein polynomials definition
+
+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`` relative to a
+simplex ``T`` are defined by
+```math
+B^{D,K}_α(\bm{x}) = \binom{K}{α} λ(\bm{x})^α\qquad\text{for all }α ∈\mathcal{I}_K^D
+```
+where
+- ``\mathcal{I}_K^D = \{\ α∈(\mathbb{Z}_+)^{D+1} \quad|\quad |α|=K\ \}``
+- ``|α|=∑_{1≤ i≤ N} α_i``
+- ``\binom{K}{α} = \frac{K!}{α_1 !α_2 !… α_N!}``
+- ``λ`` are the barycentric coordinates relative to ``T`` (defined above)
+
+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).
+
+### Bernstein indices and indexing
+
+Working with Bernstein polynomials requires dealing with several quantities
+indexed by some ``α ∈ \mathcal{I}_K^D``, the polynomials themselves but also the
+coefficients ``c_α`` of a polynomial in the basis, the domain points
+``{\bm{x}_α = \underset{1≤i≤N}{∑} α_i v_i}`` and the intermediate
+coefficients used in the de Casteljau algorithm.
+
+These indices are returned by [`bernstein_terms(K,D)`](@ref bernstein_terms).
+When storing such quantities in arrays, ``∙_α`` is stored at index
+[`bernstein_term_id(α)`](@ref bernstein_term_id), which is the index of `α`
+in `bernstein_terms(sum(α),length(α)-1)`.
+
+We adopt the convention that a quantity indexed by a ``α ∉ ℤ_+^N`` is equal to
+zero (to simplify the definition of algorithms where ``α=β-e_i`` appears).
+
+### The de Casteljau algorithms
+
+A polynomial ``p ∈ ℙ^D_K`` in Bernstein form ``p = ∑_{α∈\mathcal{I}^D_K}\, p_α
+B_α`` can be evaluated at ``\bm{x}`` using the de Casteljau algorithms
+[1, Algo. 2.9] by iteratively computing
+```math
+\qquad p_β^{(l)} = \underset{1 ≤ i ≤ N}{∑} λ_i\, p_{β+e_i}^{(l-1)} \qquad ∀β ∈ \mathcal{I}^D_{K-l},
+```
+for ``l=1, 2, …, K`` where ``p_α^{(0)}=p_α``, ``λ=λ(\bm{x})`` and the
+result is ``p(\bm{x})=p_𝟎^{(K)}``. This algorithm is implemented (in
+place) by [`_de_Casteljau_nD!`](@ref).
+
+But Gridap implements the polynomial bases themselves instead of individual
+polynomials in a basis. To compute all ``B_α`` at ``\bm{x}``, one can
+use the de Casteljau algorithm going "downwards" (from the tip of the pyramid
+to the base). The idea is to use the relation
+```math
+B_α = ∑_{1 ≤ i ≤ N} λ_i B_{α-e_i}\qquad ∀α ∈ ℤ_+^N,\ |α|≥1.
+```
+
+Starting from ``b_𝟎^{(0)}=B_𝟎(\bm{x})=1``, compute iteratively
+```math
+\qquad b_β^{(l)} = \underset{1 ≤ i ≤ N}{∑} λ_i\, b_{β-e_i}^{(l-1)} \qquad ∀β ∈ \mathcal{I}^D_{l},
+```
+for ``l=1,2, …, K``, where again ``λ=λ(\bm{x})`` and the result is
+``B_α(\bm{x})=b_α^{(K)}`` for all ``α`` in ``\mathcal{I}^D_K``. This
+algorithm is implemented (in place) by [`_downwards_de_Casteljau_nD!`](@ref).
+The implementation is a bit tricky, because the iterations must be done in
+reverse order to avoid erasing coefficients needed later, and a lot of summands
+disappear (when ``(β-e_i)_i < 0``).
+
+The gradient and hessian of the `BernsteinBasisOnSimplex` are also implemented.
+They rely on the following
+```math
+∂_q B_α(\bm{x}) = K\!∑_{1 ≤ i ≤ N} ∂_qλ_i\, B_{α-e_i}(\bm{x}),\qquad
+∂_t ∂_q B_α(\bm{x}) = K\!∑_{1 ≤ i,j ≤ N} ∂_tλ_j\, ∂_qλ_i\, B_{α-e_i-e_j}(\bm{x}).
+```
+The gradient formula comes from [1, Eq. (2.28)], and the second is derived from
+the first using the fact that ``∂_qλ`` is homogeneous. The implementation of
+the gradient and hessian compute the ``B_β`` using
+`_downwards_de_Casteljau_nD!` up to order ``K-1`` and ``K-2`` respectively, and
+then the results are assembled by [`_grad_Bα_from_Bαm!`](@ref) and
+[`_hess_Bα_from_Bαmm!`](@ref) respectively. The implementation makes sure to
+only access each relevant ``B_β`` once per ``(∇/H)B_α`` computed. Also, on the
+reference simplex, the barycentric coordinates derivatives are computed at
+compile time using ``∂_qλ_i = δ_{i q}-δ_{i N}``.
+
+## References
+
+[1] [M.J. Lai & L.L. Schumaker, Spline Functions on Triangulations, Chapter 2 - Bernstein–Bézier Methods for Bivariate Polynomials, pp. 18 - 61.](https://doi.org/10.1017/CBO9780511721588.003)

src/Polynomials/UniformPolyBases.jl → src/Polynomials/CartProdPolyBases.jl

@@ -3,34 +3,36 @@
 #################################
 
 """
-    struct UniformPolyBasis{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
+    struct CartProdPolyBasis{D,V,PT} <: PolynomialBasis{D,V,PT}
 
-Type representing a uniform basis of (an)isotropic `D`-multivariate `V`-valued
-polynomial space
+"Cartesian product of a scalar tensor polynomial basis"
+
+Type representing a basis of a (an)isotropic `D`-multivariate `V`-valued
+cartesian product polynomial space
 
 `V`(𝕊, 𝕊, ..., 𝕊)
 
-where 𝕊 is a scalar multivariate polynomial space. So each (independant)
-component of `V` holds the same space (hence the name 'uniform').
+where the scalar space 𝕊 is a (subspace of a) tensor product space of an
+univariate polynomial basis.
 
 The scalar polynomial basis spanning 𝕊 is defined as
 
   { x ⟶ bα`ᴷ`(x) = bα₁`ᴷ`(x₁) × bα₂`ᴷ`(x₂) × ... × bα`Dᴷ`(x`D`) |  α ∈ `terms` }
 
 where bαᵢ`ᴷ`(xᵢ) is the αᵢth 1D basis polynomial of the basis `PT` of order `K`
 evaluated at xᵢ (iᵗʰ comp. of x), and where α = (α₁, α₂, ..., α`D`) is a
-multi-index in `terms`, a subset of ⟦0,`K`⟧`ᴰ`. `terms` is a field that can be
+multi-index in `terms`, a subset of {0:`K`}`ᴰ`. `terms` is a field that can be
 passed in a constructor.
 
-The fields of this `struct` are not public.
 This type fully implements the [`Field`](@ref) interface, with up to second
 order derivatives.
 """
-struct UniformPolyBasis{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
+struct CartProdPolyBasis{D,V,PT} <: PolynomialBasis{D,V,PT}
+  max_order::Int
   orders::NTuple{D,Int}
   terms::Vector{CartesianIndex{D}}
 
-  function UniformPolyBasis{D}(
+  function CartProdPolyBasis{D}(
     ::Type{PT},
     ::Type{V},
     orders::NTuple{D,Int},
@@ -41,24 +43,25 @@ struct UniformPolyBasis{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
     K = maximum(orders; init=0)
     msg =  "Some term contain a higher index than the maximum degree + 1."
     @check all( term -> (maximum(Tuple(term), init=0) <= K+1), terms) msg
-    new{D,V,K,PT}(orders,terms)
+    new{D,V,PT}(K,orders,terms)
   end
 end
 
-@inline Base.size(a::UniformPolyBasis{D,V}) where {D,V} = (length(a.terms)*num_indep_components(V),)
+@inline Base.size(a::CartProdPolyBasis{D,V}) where {D,V} = (length(a.terms)*num_indep_components(V),)
+@inline get_order(b::CartProdPolyBasis) = b.max_order
 
-function UniformPolyBasis(
+function CartProdPolyBasis(
    ::Type{PT},
    ::Val{D},
    ::Type{V},
    orders::NTuple{D,Int},
    terms::Vector{CartesianIndex{D}}) where {PT<:Polynomial,D,V}
 
-  UniformPolyBasis{D}(PT,V,orders,terms)
+  CartProdPolyBasis{D}(PT,V,orders,terms)
 end
 
 """
-    UniformPolyBasis(::Type{PT}, ::Val{D}, ::Type{V}, orders::Tuple [, filter=_q_filter])
+    CartProdPolyBasis(::Type{PT}, ::Val{D}, ::Type{V}, orders::Tuple [, filter=_q_filter])
 
 This constructor allows to pass a tuple `orders` containing the maximal
 polynomial order to be used in each of the `D` spatial dimensions in order to
@@ -67,35 +70,35 @@ construct a tensorial anisotropic `D`-multivariate space 𝕊.
 If a filter is provided, it is applied on the cartesian product terms
 CartesianIndices(`orders`), with maximum(`orders`) as order argument.
 """
-function UniformPolyBasis(
+function CartProdPolyBasis(
   ::Type{PT}, ::Val{D}, ::Type{V}, orders::NTuple{D,Int}, filter::Function=_q_filter
   ) where {PT,D,V}
 
   terms = _define_terms(filter, orders)
-  UniformPolyBasis{D}(PT,V,orders,terms)
+  CartProdPolyBasis{D}(PT,V,orders,terms)
 end
 
 """
-    UniformPolyBasis(::Type{PT}, ::Type{V}, ::Val{D}, order::Int [, filter=_q_filter])
+    CartProdPolyBasis(::Type{PT}, ::Type{V}, ::Val{D}, order::Int [, filter=_q_filter])
 
-Return a `UniformPolyBasis{D,V,order,PT}` where 𝕊 is defined by the terms
+Return a `CartProdPolyBasis{D,V,order,PT}` where 𝕊 is defined by the terms
 filtered by
 
     term -> filter(term, order).
 
 See the [Filter functions](@ref) section of the documentation for more details.
 """
-function UniformPolyBasis(
+function CartProdPolyBasis(
   ::Type{PT}, VD::Val{D}, ::Type{V}, order::Int, filter::Function=_q_filter) where {PT,D,V}
 
   orders = tfill(order,VD)
-  UniformPolyBasis(PT,Val(D),V,orders,filter)
+  CartProdPolyBasis(PT,Val(D),V,orders,filter)
 end
 
 # API
 
 """
-    get_exponents(b::UniformPolyBasis)
+    get_exponents(b::CartProdPolyBasis)
 
 Get a vector of tuples with the exponents of all the terms in the basis of 𝕊,
 the components scalar space of `b`.
@@ -115,17 +118,17 @@ println(exponents)
 Tuple{Int,Int}[(0, 0), (1, 0), (2, 0), (0, 1), (1, 1), (2, 1), (0, 2), (1, 2), (2, 2)]
 ```
 """
-function get_exponents(b::UniformPolyBasis)
+function get_exponents(b::CartProdPolyBasis)
   indexbase = 1
-  [Tuple(t) .- indexbase for t in b.terms]
+  Tuple(Tuple(t) .- indexbase for t in b.terms)
 end
 
 """
-    get_orders(b::UniformPolyBasis)
+    get_orders(b::CartProdPolyBasis)
 
 Return the D-tuple of polynomial orders in each spatial dimension
 """
-function get_orders(b::UniformPolyBasis)
+function get_orders(b::CartProdPolyBasis)
   b.orders
 end
 
@@ -134,42 +137,42 @@ end
 #################################
 
 function _evaluate_nd!(
-  b::UniformPolyBasis{D,V,K,PT}, x,
+  b::CartProdPolyBasis{D,V,PT}, x,
   r::AbstractMatrix{V}, i,
-  c::AbstractMatrix{T}) where {D,V,K,PT,T}
-
-  terms  = b.terms
-  orders = b.orders
+  c::AbstractMatrix{T}, VK::Val) where {D,V,PT,T}
 
   for d in 1:D
-    Kd = Val(orders[d])
-    _evaluate_1d!(PT,Kd,c,x,d)
+    # The optimization below of fine tuning Kd for `orders` is a bottlneck if
+    # `orders` are not passed in `Val`s, due to runtime dispatch depending on
+    # none inferable Val(b.orders[d])
+    # Kd = Val(b.orders[d])
+    _evaluate_1d!(PT,VK,c,x,d)
   end
 
   k = 1
-  for ci in terms
+  for ci in b.terms
 
     s = one(T)
     for d in 1:D
       @inbounds s *= c[d,ci[d]]
     end
 
-    k = _uniform_set_value!(r,i,s,k)
+    k = _cartprod_set_value!(r,i,s,k)
   end
 end
 
 """
-    _uniform_set_value!(r::AbstractMatrix{<:Real},i,s,k)
+    _cartprod_set_value!(r::AbstractMatrix{<:Real},i,s,k)
 
 r[i,k] = s; return k+1
 """
-function _uniform_set_value!(r::AbstractMatrix{<:Real},i,s,k)
+function _cartprod_set_value!(r::AbstractMatrix{<:Real},i,s,k)
   @inbounds r[i,k] = s
   k+1
 end
 
 """
-    _uniform_set_value!(r::AbstractMatrix{V},i,s::T,k,l)
+    _cartprod_set_value!(r::AbstractMatrix{V},i,s::T,k,l)
 
 ```
 r[i,k]     = V(s, 0,    ..., 0)
@@ -181,7 +184,7 @@ return k+N
 
 where `N = num_indep_components(V)`.
 """
-function _uniform_set_value!(r::AbstractMatrix{V},i,s::T,k) where {V,T}
+function _cartprod_set_value!(r::AbstractMatrix{V},i,s::T,k) where {V,T}
   ncomp = num_indep_components(V)
   z = zero(T)
   @inbounds for j in 1:ncomp
@@ -192,22 +195,18 @@ function _uniform_set_value!(r::AbstractMatrix{V},i,s::T,k) where {V,T}
 end
 
 function _gradient_nd!(
-  b::UniformPolyBasis{D,V,K,PT}, x,
+  b::CartProdPolyBasis{D,V,PT}, x,
   r::AbstractMatrix{G}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
-  s::MVector{D,T}) where {D,V,K,PT,G,T}
-
-  terms  = b.terms
-  orders = b.orders
+  s::MVector{D,T}, VK::Val) where {D,V,PT,G,T}
 
   for d in 1:D
-    Kd = Val(orders[d])
-    _derivatives_1d!(PT,Kd,(c,g),x,d)
+    _derivatives_1d!(PT,VK,(c,g),x,d)
   end
 
   k = 1
-  for ci in terms
+  for ci in b.terms
 
     for i in eachindex(s)
       @inbounds s[i] = one(T)
@@ -223,12 +222,12 @@ function _gradient_nd!(
       end
     end
 
-    k = _uniform_set_derivative!(r,i,s,k,V)
+    k = _cartprod_set_derivative!(r,i,s,k,V)
   end
 end
 
 """
-    _uniform_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{<:Real})
+    _cartprod_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{<:Real})
 
 ```
 r[i,k] = s = (∇bᵏ)(xi); return k+1
@@ -237,15 +236,15 @@ r[i,k] = s = (∇bᵏ)(xi); return k+1
 where bᵏ is the kᵗʰ basis polynomial. Note that `r[i,k]` is a `VectorValue` or
 `TensorValue` and `s` a `MVector` or `MMatrix` respectively, of same size.
 """
-function _uniform_set_derivative!(
+function _cartprod_set_derivative!(
   r::AbstractMatrix{G},i,s,k,::Type{<:Real}) where G
 
   @inbounds r[i,k] = s
   k+1
 end
 
 """
-    _uniform_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{V})
+    _cartprod_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{V})
 
 ```
 z = zero(s)
@@ -259,7 +258,7 @@ return k+n
 Note that `r[i,k]` is a `TensorValue` or `ThirdOrderTensorValue` and `s` a
 `MVector` or `MMatrix`.
 """
-@generated function _uniform_set_derivative!(
+@generated function _cartprod_set_derivative!(
   r::AbstractMatrix{G},i,s,k,::Type{V}) where {G,V}
   # Git blame me for readable non-generated version
 
@@ -294,7 +293,7 @@ end
 # necessary as long as outer(Point, V<:AbstractSymTensorValue)::G does not
 # return a tensor type G that implements the appropriate symmetries of the
 # gradient (and hessian)
-@generated function _uniform_set_derivative!(
+@generated function _cartprod_set_derivative!(
   r::AbstractMatrix{G},i,s,k,::Type{V}) where {G,V<:AbstractSymTensorValue{D}} where D
   # Git blame me for readable non-generated version
 
@@ -333,24 +332,20 @@ end
 end
 
 function _hessian_nd!(
-  b::UniformPolyBasis{D,V,K,PT}, x,
+  b::CartProdPolyBasis{D,V,PT}, x,
   r::AbstractMatrix{G}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
   h::AbstractMatrix{T},
-  s::MMatrix{D,D,T}) where {D,V,K,PT,G,T}
-
-  terms  = b.terms
-  orders = b.orders
+  s::MMatrix{D,D,T}, VK::Val) where {D,V,PT,G,T}
 
   for d in 1:D
-    Kd = Val(orders[d])
-    _derivatives_1d!(PT,Kd,(c,g,h),x,d)
+    _derivatives_1d!(PT,VK,(c,g,h),x,d)
   end
 
   k = 1
 
-  for ci in terms
+  for ci in b.terms
 
     for i in eachindex(s)
       @inbounds s[i] = one(T)
@@ -370,7 +365,7 @@ function _hessian_nd!(
       end
     end
 
-    k = _uniform_set_derivative!(r,i,s,k,V)
+    k = _cartprod_set_derivative!(r,i,s,k,V)
   end
 end
 

src/Polynomials/ChebyshevBases.jl

@@ -11,11 +11,11 @@ struct Chebyshev{kind} <: Polynomial end
 isHierarchical(::Type{<:Chebyshev}) = true
 
 """
-    ChebyshevBasis{D,V,kind,K} = UniformPolyBasis{D,V,K,Chebyshev{kind}}
+    ChebyshevBasis{D,V,kind} = CartProdPolyBasis{D,V,Chebyshev{kind}}
 
-Alias for Chebyshev multivariate scalar' or `Multivalue`'d basis.
+Alias for cartesian product Chebyshev basis, scalar valued or multivalued.
 """
-const ChebyshevBasis{D,V,kind,K} = UniformPolyBasis{D,V,K,Chebyshev{kind}}
+const ChebyshevBasis{D,V,kind} = CartProdPolyBasis{D,V,Chebyshev{kind}}
 
 """
     ChebyshevBasis(::Val{D}, ::Type{V}, order::Int, terms::Vector; kind=:T)
@@ -24,9 +24,9 @@ const ChebyshevBasis{D,V,kind,K} = UniformPolyBasis{D,V,K,Chebyshev{kind}}
 
 High level constructors of [`ChebyshevBasis`](@ref).
 """
-ChebyshevBasis(args...; kind=:T) = UniformPolyBasis(Chebyshev{kind}, args...)
+ChebyshevBasis(args...; kind=:T) = CartProdPolyBasis(Chebyshev{kind}, args...)
 
-function UniformPolyBasis(
+function CartProdPolyBasis(
   ::Type{Chebyshev{:U}}, ::Val{D}, ::Type{V}, ::Int) where {D, V}
 
   @notimplemented "1D evaluation for second kind need to be implemented here"

src/Polynomials/CompWiseTensorPolyBases.jl

@@ -1,5 +1,5 @@
 """
-    CompWiseTensorPolyBasis{D,V,K,PT,L} <: PolynomialBasis{D,V,K,PT}
+    CompWiseTensorPolyBasis{D,V,PT,L} <: PolynomialBasis{D,V,PT}
 
 "Polynomial basis of component wise tensor product polynomial spaces"
 
@@ -20,37 +20,39 @@ with `L`>1, where the scalar `D`-multivariate spaces 𝕊ˡ (for 1 ≤ l ≤ `L`
 The `L`×`D` matrix of orders α is given in the constructor, and `K` is the
 maximum of α. Any 1D polynomial family `PT<:Polynomial` is usable.
 """
-struct CompWiseTensorPolyBasis{D,V,K,PT,L} <: PolynomialBasis{D,V,K,PT}
+struct CompWiseTensorPolyBasis{D,V,PT,L} <: PolynomialBasis{D,V,PT}
+  max_order::Int
   orders::SMatrix{L,D,Int}
 
   function CompWiseTensorPolyBasis{D}(
     ::Type{PT}, ::Type{V}, orders::SMatrix{L,D,Int}) where {D,PT<:Polynomial,V,L}
 
     msg1 = "The orders matrix rows number must match the number of independent components of V"
     @check L == num_indep_components(V) msg1
-    msg2 = "The Component Wise construction is useless for one component, use UniformPolyBasis instead"
+    msg2 = "The Component Wise construction is useless for one component, use CartProdPolyBasis instead"
     @check L > 1 msg2
     @check D > 0
     @check isconcretetype(PT) "PT needs to be a concrete <:Polynomial type"
     K = maximum(orders)
 
-    new{D,V,K,PT,L}(orders)
+    new{D,V,PT,L}(K,orders)
   end
 end
 
 Base.size(a::CompWiseTensorPolyBasis) = ( sum(prod.(eachrow(a.orders .+ 1))), )
+get_order(b::CompWiseTensorPolyBasis) = b.max_order
 
 """
     get_comp_terms(f::CompWiseTensorPolyBasis{D,V})
 
 Return a tuple (terms\\_1, ..., terms\\_l, ..., terms\\_L) containing, for each
 component of V, the Cartesian indices iterator over the terms that define 𝕊ˡ,
-that is all elements of ⟦1,`o`(l,1)+1⟧ × ⟦1,`o`(l,2)+1⟧ × … × ⟦1,`o`(l,D)+1⟧.
+that is all elements of {1 : `o`(l,1)+1} × {1 : `o`(l,2)+1} × … × {1 : `o`(l,D)+1}.
 
 E.g., if `orders=[ 0 1; 1 0]`, then the `comp_terms` are
 `( CartesianIndices{2}((1,2)), CartesianIndices{2}((2,1)) )`.
 """
-function get_comp_terms(f::CompWiseTensorPolyBasis{D,V,K,PT,L}) where {D,V,K,PT,L}
+function get_comp_terms(f::CompWiseTensorPolyBasis{D,V,PT,L}) where {D,V,PT,L}
   _terms(l) = CartesianIndices( Tuple(f.orders[l,:] .+ 1) )
   comp_terms = ntuple(l -> _terms(l), Val(L))
   comp_terms::NTuple{L,CartesianIndices{D}}
@@ -62,24 +64,21 @@ end
 #################################
 
 function _evaluate_nd!(
-  b::CompWiseTensorPolyBasis{D,V,K,PT,L}, x,
+  b::CompWiseTensorPolyBasis{D,V,PT,L}, x,
   r::AbstractMatrix{V}, i,
-  c::AbstractMatrix{T}) where {D,V,K,PT,L,T}
-
-  orders = b.orders
-  comp_terms = get_comp_terms(b)
+  c::AbstractMatrix{T}, VK::Val) where {D,V,PT,L,T}
 
   for d in 1:D
-    # for each coordinate d, the order at which the basis should be evaluated is
-    # the maximum d-order for any component l
-    Kd = Val(maximum(orders[:,d]))
-    _evaluate_1d!(PT,Kd,c,x,d)
+    # The optimization below of fine tuning Kd is a bottlneck if not put in a
+    # function due to runtime dispatch and creation of Val(Kd)
+    #  # for each coordinate d, the order at which the basis should be evaluated is
+    #  # the maximum d-order for any component l
+    #  Kd = Val(maximum(b.orders[:,d]))
+    _evaluate_1d!(PT,VK,c,x,d)
   end
 
-  m = zero(Mutable(V))
   k = 1
-
-  for (l,terms) in enumerate(comp_terms)
+  for (l,terms) in enumerate(get_comp_terms(b))
     for ci in terms
 
       s = one(T)
@@ -109,25 +108,18 @@ function _comp_wize_set_value!(r::AbstractMatrix{V},i,s::T,k,l) where {V,T}
 end
 
 function _gradient_nd!(
-  b::CompWiseTensorPolyBasis{D,V,K,PT,L}, x,
+  b::CompWiseTensorPolyBasis{D,V,PT,L}, x,
   r::AbstractMatrix{G}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
-  s::MVector{D,T}) where {D,V,K,PT,L,G,T}
-
-  orders = b.orders
-  comp_terms = get_comp_terms(b)
+  s::MVector{D,T}, VK::Val) where {D,V,PT,L,G,T}
 
   for d in 1:D
-    # for each spatial coordinate d, the order at which the basis should be
-    # evaluated is the maximum d-order for any component l
-    Kd = Val(maximum(orders[:,d]))
-    _derivatives_1d!(PT,Kd,(c,g),x,d)
+    _derivatives_1d!(PT,VK,(c,g),x,d)
   end
 
   k = 1
-
-  for (l,terms) in enumerate(comp_terms)
+  for (l,terms) in enumerate(get_comp_terms(b))
     for ci in terms
 
       for i in eachindex(s)
@@ -150,7 +142,7 @@ function _gradient_nd!(
 end
 
 """
-    _comp_wize_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{V})
+    _comp_wize_set_derivative!(r::AbstractMatrix{G},i,s,k,::Val{l},::Type{V})
 
 ```
 z = zero(s)
@@ -180,34 +172,27 @@ the `k`ᵗʰ basis polynomial, whose nonzero component in `V` is the `l`ᵗʰ.
   return Expr(:block, body ,:(return k+1))
 end
 
-# See _uniform_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{V}) where {G,V<:AbstractSymTensorValue{D}} where D
+# See _cartprod_set_derivative!(r::AbstractMatrix{G},i,s,k,::Type{V}) where {G,V<:AbstractSymTensorValue{D}} where D
 @generated function _comp_wize_set_derivative!(
   r::AbstractMatrix{G},i,s,k,::Type{V}) where {G,V<:AbstractSymTensorValue{D}} where D
 
   @notimplemented
 end
 
 function _hessian_nd!(
-  b::CompWiseTensorPolyBasis{D,V,K,PT,L}, x,
+  b::CompWiseTensorPolyBasis{D,V,PT,L}, x,
   r::AbstractMatrix{H}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
   h::AbstractMatrix{T},
-  s::MMatrix{D,D,T}) where {D,V,K,PT,L,H,T}
-
-  orders = b.orders
-  comp_terms = get_comp_terms(b)
+  s::MMatrix{D,D,T}, VK::Val) where {D,V,PT,L,H,T}
 
   for d in 1:D
-    # for each spatial coordinate d, the order at which the basis should be
-    # evaluated is the maximum d-order for any component l
-    Kd = Val(maximum(orders[:,d]))
-    _derivatives_1d!(PT,Kd,(c,g,h),x,d)
+    _derivatives_1d!(PT,VK,(c,g,h),x,d)
   end
 
   k = 1
-
-  for (l,terms) in enumerate(comp_terms)
+  for (l,terms) in enumerate(get_comp_terms(b))
     for ci in terms
 
       for i in eachindex(s)
@@ -263,7 +248,7 @@ b = QGradBasis(Monomial, Val(3), Float64, 2)
 For more details, see [`CompWiseTensorPolyBasis`](@ref), as `QGradBasis` returns
 an instance of\\
 `CompWiseTensorPolyBasis{D, VectorValue{D,T}, order+1, PT}` for `D`>1, or\\
-`UniformPolyBasis{1, VectorValue{1,T}, order+1, PT}` for `D`=1.
+`CartProdPolyBasis{1, VectorValue{1,T}, order+1, PT}` for `D`=1.
 """
 function QGradBasis(::Type{PT},::Val{D},::Type{T},order::Int) where {PT,D,T}
   @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
@@ -278,7 +263,7 @@ function QGradBasis(::Type{PT},::Val{1},::Type{T},order::Int) where {PT,T}
   @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
 
   V = VectorValue{1,T}
-  UniformPolyBasis(PT, Val(1), V, order+1)
+  CartProdPolyBasis(PT, Val(1), V, order+1)
 end
 
 
@@ -311,7 +296,7 @@ b = QCurlGradBasis(Bernstein, Val(2), Float64, 3)
 For more details, see [`CompWiseTensorPolyBasis`](@ref), as `QCurlGradBasis`
 returns an instance of\\
 `CompWiseTensorPolyBasis{D, VectorValue{D,T}, order+1, PT}` for `D`>1, or\\
-`UniformPolyBasis{1, VectorValue{1,T}, order+1, PT}` for `D`=1.
+`CartProdPolyBasis{1, VectorValue{1,T}, order+1, PT}` for `D`=1.
 """
 function QCurlGradBasis(::Type{PT},::Val{D},::Type{T},order::Int) where {PT,D,T}
   @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
@@ -326,5 +311,5 @@ function QCurlGradBasis(::Type{PT},::Val{1},::Type{T},order::Int) where {PT,T}
   @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
 
   V = VectorValue{1,T}
-  UniformPolyBasis(PT, Val(1), V, order+1)
+  CartProdPolyBasis(PT, Val(1), V, order+1)
 end

src/Polynomials/LegendreBases.jl

@@ -8,11 +8,11 @@ struct Legendre <: Polynomial   end
 isHierarchical(::Type{Legendre}) = true
 
 """
-    LegendreBasis{D,V,K} = UniformPolyBasis{D,V,K,Legendre}
+    LegendreBasis{D,V} = CartProdPolyBasis{D,V,Legendre}
 
-Alias for Legendre multivariate scalar' or `Multivalue`'d basis.
+Alias for cartesian product Legendre basis, scalar valued or multi-valued.
 """
-const LegendreBasis{D,V,K} = UniformPolyBasis{D,V,K,Legendre}
+const LegendreBasis{D,V} = CartProdPolyBasis{D,V,Legendre}
 
 """
     LegendreBasis(::Val{D}, ::Type{V}, order::Int, terms::Vector)
@@ -21,7 +21,7 @@ const LegendreBasis{D,V,K} = UniformPolyBasis{D,V,K,Legendre}
 
 High level constructors of [`LegendreBasis`](@ref).
 """
-LegendreBasis(args...) = UniformPolyBasis(Legendre, args...)
+LegendreBasis(args...) = CartProdPolyBasis(Legendre, args...)
 
 
 # 1D evaluation implementation

src/Polynomials/ModalC0Bases.jl

@@ -8,13 +8,14 @@ Reference: Eq. (17) in https://doi.org/10.1016/j.camwa.2022.09.027
 struct ModalC0 <: Polynomial end
 
 """
-    ModalC0Basis{D,V,T,K} <: PolynomialBasis{D,V,K,ModalC0}
+    ModalC0Basis{D,V,T} <: PolynomialBasis{D,V,ModalC0}
 
 Tensor product basis of generalised modal C0 1D basis from section 5.2 in
 https://doi.org/10.1016/j.camwa.2022.09.027.
 See also [ModalC0 polynomials](@ref) section of the documentation.
 """
-struct ModalC0Basis{D,V,T,K} <: PolynomialBasis{D,V,K,ModalC0}
+struct ModalC0Basis{D,V,T} <: PolynomialBasis{D,V,ModalC0}
+  max_order::Int
   orders::NTuple{D,Int}
   terms::Vector{CartesianIndex{D}}
   a::Vector{Point{D,T}}
@@ -32,7 +33,7 @@ struct ModalC0Basis{D,V,T,K} <: PolynomialBasis{D,V,K,ModalC0}
     @check T == eltype(V) "Point and polynomial values should have the same scalar body"
     K = maximum(orders, init=0)
 
-    new{D,V,T,K}(orders,terms,a,b)
+    new{D,V,T}(K,orders,terms,a,b)
   end
 end
 
@@ -50,8 +51,8 @@ At last, all scalar basis polynomial will have its bounding box `(a[i],b[i])`,
 but they are assumed iddentical if only two points `a` and `b` are provided,
 and default to `a=Point{D}(0...)`, `b=Point{D}(1...)` if not provided.
 
-The basis is uniform, isotropic if one `order` is provided, or anisotropic if a
-`D` tuple `orders` is provided.
+The basis is a cartesian product when multi-valued, isotropic if one `order` is
+provided, or anisotropic if a `D` tuple `orders` is provided.
 """
 function ModalC0Basis() end
 
@@ -118,7 +119,8 @@ end
 
 # API
 
-@inline Base.size(a::ModalC0Basis{D,V}) where {D,V} = (length(a.terms)*num_indep_components(V),)
+@inline Base.size(b::ModalC0Basis{D,V}) where {D,V} = (length(b.terms)*num_indep_components(V),)
+@inline get_order(b::ModalC0Basis) = b.max_order
 
 function get_orders(b::ModalC0Basis)
   b.orders
@@ -177,10 +179,12 @@ end
 # nD evaluations implementation #
 #################################
 
+_get_static_parameters(::ModalC0Basis) = nothing
+
 function _evaluate_nd!(
-  basis::ModalC0Basis{D,V,T,K}, x,
+  basis::ModalC0Basis{D,V,T}, x,
   r::AbstractMatrix{V}, i,
-  c::AbstractMatrix{T}) where {D,V,T,K}
+  c::AbstractMatrix{T}, ::Nothing) where {D,V,T}
 
   terms  = basis.terms
   orders = basis.orders
@@ -221,11 +225,11 @@ end
 end
 
 function _gradient_nd!(
-  basis::ModalC0Basis{D,V,T,K}, x,
+  basis::ModalC0Basis{D,V,T}, x,
   r::AbstractMatrix{G}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
-  s::MVector{D,T}) where {D,V,T,K,G}
+  s::MVector{D,T}, ::Nothing) where {D,V,T,G}
 
   terms  = basis.terms
   orders = basis.orders
@@ -302,12 +306,12 @@ end
 end
 
 function _hessian_nd!(
-  basis::ModalC0Basis{D,V,T,K}, x,
+  basis::ModalC0Basis{D,V,T}, x,
   r::AbstractMatrix{G}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
   h::AbstractMatrix{T},
-  s::MMatrix{D,D,T}) where {D,V,T,K,G}
+  s::MMatrix{D,D,T}, ::Nothing) where {D,V,T,G}
 
   terms  = basis.terms
   orders = basis.orders
@@ -414,7 +418,7 @@ end
 # Generic 1D internal polynomial APIs #
 #######################################
 
-# For possible use with UniformPolyBasis etc.
+# For possible use with CartProdPolyBasis etc.
 # Make it for x∈[0,1] like the other 1D bases.
 
 function _evaluate_1d!(::Type{ModalC0},::Val{K},c::AbstractMatrix{T},x,d) where {K,T<:Number}

src/Polynomials/MonomialBases.jl

@@ -8,11 +8,11 @@ struct Monomial <: Polynomial   end
 isHierarchical(::Type{Monomial}) = true
 
 """
-    MonomialBasis{D,V,K} = UniformPolyBasis{D,V,K,Monomial}
+    MonomialBasis{D,V} = CartProdPolyBasis{D,V,Monomial}
 
-Alias for monomial Multivariate scalar' or `Multivalue`'d basis.
+Alias for cartesian product monomial basis, scalar valued or multi-valued.
 """
-const MonomialBasis{D,V,K} = UniformPolyBasis{D,V,K,Monomial}
+const MonomialBasis{D,V} = CartProdPolyBasis{D,V,Monomial}
 
 """
     MonomialBasis(::Val{D}, ::Type{V}, order::Int, terms::Vector)
@@ -21,15 +21,15 @@ const MonomialBasis{D,V,K} = UniformPolyBasis{D,V,K,Monomial}
 
 High level constructors of [`MonomialBasis`](@ref).
 """
-MonomialBasis(args...) = UniformPolyBasis(Monomial, args...)
+MonomialBasis(args...) = CartProdPolyBasis(Monomial, args...)
 
 function PGradBasis(::Type{Monomial},::Val{D},::Type{T},order::Int) where {D,T}
   NedelecPolyBasisOnSimplex{D}(Monomial,T,order)
 end
 function PGradBasis(::Type{Monomial},::Val{1},::Type{T},order::Int) where T
   @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
   V = VectorValue{1,T}
-  UniformPolyBasis(Monomial, Val(1), V, order+1)
+  CartProdPolyBasis(Monomial, Val(1), V, order+1)
 end
 
 # 1D evaluation implementation
@@ -78,3 +78,63 @@ function _hessian_1d!(::Type{Monomial},::Val{K},h::AbstractMatrix{T},x,d) where
   end
 end
 
+
+# Optimizations for 0 to 1/2 derivatives at once
+
+function _derivatives_1d!(::Type{Monomial},v::Val_01,t::NTuple{2},x,d)
+  @inline _evaluate_1d!(Monomial, v, t[1], x, d)
+  @inline _gradient_1d!(Monomial, v, t[2], x, d)
+end
+
+function _derivatives_1d!(::Type{Monomial},::Val{K},t::NTuple{2},x,d) where K
+  @inbounds begin
+    n = K + 1 # n > 2
+    v, g = t
+    T = eltype(v)
+
+    z = zero(T)
+    xn = one(T)
+    xd = x[d]
+
+    v[d,1] = xn
+    g[d,1] = z
+    for i in 2:n
+      g[d,i] = (i-1)*xn
+      xn *= xd
+      v[d,i] = xn
+    end
+  end
+end
+
+
+function _derivatives_1d!(::Type{Monomial},v::Val_012,t::NTuple{3},x,d)
+  @inline _evaluate_1d!(Monomial, v, t[1], x, d)
+  @inline _gradient_1d!(Monomial, v, t[2], x, d)
+  @inline _hessian_1d!( Monomial, v, t[3], x, d)
+end
+
+function _derivatives_1d!(::Type{Monomial},::Val{K},t::NTuple{3},x,d) where K
+  @inbounds begin
+    n = K + 1 # n > 2
+    v, g, h = t
+    T = eltype(v)
+
+    z = zero(T)
+    o = one(T)
+    xd = x[d]
+
+    v[d,1] = o;  g[d,1] = z; h[d,1] = z
+    v[d,2] = xd; g[d,2] = o; h[d,2] = z
+
+    xn  = xd
+    xnn = o
+    for i in 3:n
+      h[d,i] = (i-1)*xnn
+      xnn  = (i-1)*xn
+      g[d,i] = xnn
+      xn *= xd
+      v[d,i] = xn
+    end
+  end
+end
+

src/Polynomials/NedelecPolyBases.jl

@@ -1,38 +1,41 @@
 """
-    NedelecPolyBasisOnSimplex{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
+    NedelecPolyBasisOnSimplex{D,V,PT} <: PolynomialBasis{D,V,PT}
 
 Basis of the vector valued (`V<:VectorValue{D}`) space ℕ𝔻ᴰₙ(△) for `D`=2,3.
 This space is the polynomial space for Nedelec elements on simplices with
 curl in (ℙᴰₙ)ᴰ. Its maximum degree is n+1 = `K`. `get_order` on it returns `K`.
 
-Currently, the basis is implemented as the union of a UniformPolyBasis{...,PT}
+Currently, the basis is implemented as the union of a CartProdPolyBasis{...,PT}
 for ℙᴰₙ and a monomial basis for x × (ℙᴰₙ \\ ℙᴰₙ₋₁)ᴰ.
 """
-struct NedelecPolyBasisOnSimplex{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
+struct NedelecPolyBasisOnSimplex{D,V,PT} <: PolynomialBasis{D,V,PT}
   order::Int
   function NedelecPolyBasisOnSimplex{D}(::Type{PT},::Type{T},order::Integer) where {D,PT<:Polynomial,T}
     @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
     @check isconcretetype(PT) "PT needs to be a concrete <:Polynomial type"
     @notimplementedif !(D in (2,3))
     K = Int(order)+1
     V = VectorValue{D,T}
-    new{D,V,K,PT}(Int(order))
+    new{D,V,PT}(Int(order))
   end
 end
 
-function Base.size(a::NedelecPolyBasisOnSimplex{D}) where D
-  k = a.order+1
-  n = div(k*prod(i->(k+i),2:D),factorial(D-1))
+get_order(f::NedelecPolyBasisOnSimplex) = f.order + 1 # Return actual maximum poly order
+
+function Base.size(f::NedelecPolyBasisOnSimplex{D}) where D
+  K = get_order(f)
+  n = div(K*prod(i->(K+i),2:D),factorial(D-1))
   (n,)
 end
 
 function return_cache(
-  f::NedelecPolyBasisOnSimplex{D,V,K,PT},x::AbstractVector{<:Point}) where {D,V,K,PT}
+  f::NedelecPolyBasisOnSimplex{D,V,PT},x::AbstractVector{<:Point}) where {D,V,PT}
 
+  K = get_order(f)
   np = length(x)
   ndofs = length(f)
   a = zeros(V,(np,ndofs))
-  P = UniformPolyBasis(PT,Val(D),V,K-1,_p_filter)
+  P = CartProdPolyBasis(PT,Val(D),V,K-1,_p_filter)
   cP = return_cache(P,x)
   CachedArray(a), cP, P
 end
@@ -115,15 +118,16 @@ function evaluate!(
 end
 
 function return_cache(
-  g::FieldGradientArray{1,<:NedelecPolyBasisOnSimplex{D,V,K,PT}},
-  x::AbstractVector{<:Point}) where {D,V,K,PT}
+  g::FieldGradientArray{1,<:NedelecPolyBasisOnSimplex{D,V,PT}},
+  x::AbstractVector{<:Point}) where {D,V,PT}
   f = g.fa
+  K = get_order(f)
   np = length(x)
   ndofs = length(f)
   xi = testitem(x)
   G = gradient_type(V,xi)
   a = zeros(G,(np,ndofs))
-  mb = UniformPolyBasis(PT,Val(D),V,K-1,_p_filter)
+  mb = CartProdPolyBasis(PT,Val(D),V,K-1,_p_filter)
   P = Broadcasting(∇)(mb)
   cP = return_cache(P,x)
   CachedArray(a), cP, P

src/Polynomials/PolynomialInterfaces.jl

@@ -7,10 +7,6 @@
 
 Abstract type for polynomial bases families/types. It has trait
 [`isHierarchical`](@ref).
-
-The currently implemented families are [Monomial](@ref), [Legendre](@ref),
-[Chebyshev](@ref), [ModalC0](@ref) and [Bernstein](@ref). Only Bernstein is not
-hierarchical.
 """
 abstract type Polynomial  <: Field end
 
@@ -20,6 +16,10 @@ abstract type Polynomial  <: Field end
 Return `true` if the 1D basis of order `K` of the given [`Polynomial`](@ref)
 basis family is the union of the basis of order `K-1` and an other order `K`
 polynomial. Equivalently, if the iᵗʰ basis polynomial is of order i-1.
+
+The currently implemented families are [Monomial](@ref), [Legendre](@ref),
+[Chebyshev](@ref), [ModalC0](@ref) and [Bernstein](@ref). Only Bernstein is not
+hierarchical.
 """
 isHierarchical(::Type{<:Polynomial}) = @abstractmethod
 
@@ -44,38 +44,39 @@ isHierarchical(::Type{<:Polynomial}) = @abstractmethod
 # ndof_1d: maximum of 1D polynomial vector in any spatial dimension
 
 """
-    PolynomialBasis{D,V,K,PT<:Polynomial} <: AbstractVector{PT}
+    PolynomialBasis{D,V,PT<:Polynomial} <: AbstractVector{PT}
 
 Abstract type representing a generic multivariate polynomial basis.
 The parameters are:
 - `D`: the spatial dimension
 - `V`: the image values type, a concrete type `<:Real` or `<:MultiValue`
-- `K`: the maximum order of a basis polynomial in a spatial component
 - `PT <: Polynomial`: the family of the basis polynomials (must be a concrete type).
+
+The implementations also stores `K`: the maximum order of a basis polynomial in a spatial component
 """
-abstract type PolynomialBasis{D,V,K,PT<:Polynomial} <: AbstractVector{PT}  end
+abstract type PolynomialBasis{D,V,PT<:Polynomial} <: AbstractVector{PT}  end
 
 @inline Base.size(::PolynomialBasis{D,V}) where {D,V} = @abstractmethod
-@inline Base.getindex(::PolynomialBasis{D,V,K,PT}, i::Integer) where {D,V,K,PT} = PT()
+@inline Base.getindex(::PolynomialBasis{D,V,PT}, i::Integer) where {D,V,PT} = PT()
 @inline Base.IndexStyle(::PolynomialBasis) = IndexLinear()
 @inline return_type(::PolynomialBasis{D,V}) where {D,V} = V
 
 """
-    get_order(b::PolynomialBasis{D,V,K}) = K
+    get_order(b::PolynomialBasis)
 
-Return the maximum polynomial order in a dimension, is should be `0` in 0D.
+Return the maximum polynomial order in a dimension, or `0` in 0D.
 """
-@inline get_order(::PolynomialBasis{D,V,K}) where {D,V,K} = K
+@inline get_order(::PolynomialBasis) = @abstractmethod
 
 
 ###########
 # Helpers #
 ###########
 
 _q_filter( e,order)  = (maximum(e,init=0) <= order) # ℚₙ
-_qs_filter(e,order)  = (maximum(e,init=0) == order) # ℚₙ\ℚ₍ₙ₋₁₎
+_qh_filter(e,order)  = (maximum(e,init=0) == order) # ℚₙ\ℚ₍ₙ₋₁₎
 _p_filter( e,order)  = (sum(e) <= order)            # ℙₙ
-_ps_filter(e,order)  = (sum(e) == order)            # ℙₙ\ℙ₍ₙ₋₁₎
+_ph_filter(e,order)  = (sum(e) == order)            # ℙₙ\ℙ₍ₙ₋₁₎
 _ser_filter(e,order) = (sum( [ i for i in e if i>1 ] ) <= order) # Serendipity
 
 function _define_terms(filter,orders)
@@ -102,10 +103,12 @@ function _return_cache(
   ndof_1d = get_order(f) + 1
   # Cache for the returned array
   r = CachedArray(zeros(G,(np,ndof)))
+  # Mutable cache for one N_deriv's derivative of a T-valued scalar polynomial
+  s = MArray{Tuple{Vararg{D,N_deriv}},T}(undef)
   # Cache for the 1D basis function values in each dimension (to be
   # tensor-producted), and of their N_deriv'th 1D derivatives
   t = ntuple( _ -> CachedArray(zeros(T,(D,ndof_1d ))), Val(N_deriv+1))
-  (r, t...)
+  (r, s, t...)
 end
 
 function return_cache(f::PolynomialBasis{D,V}, x::AbstractVector{<:Point}) where {D,V}
@@ -135,16 +138,25 @@ function _setsize!(f::PolynomialBasis{D}, np, r, t...) where D
   end
 end
 
+"""
+    _get_static_parameters(::PolynomialBasis)
+
+Return a (tuple of) static parameter(s) appended to low level `[...]_nd!` evaluation
+calls, default is `Val(get_order(b))`.
+"""
+_get_static_parameters(b::PolynomialBasis) = Val(get_order(b))
+
 function evaluate!(cache,
   f::PolynomialBasis,
   x::AbstractVector{<:Point})
 
-  r, c = cache
+  r, _, c = cache
   np = length(x)
   _setsize!(f,np,r,c)
+  params = _get_static_parameters(f)
   for i in 1:np
     @inbounds xi = x[i]
-    _evaluate_nd!(f,xi,r,i,c)
+    _evaluate_nd!(f,xi,r,i,c,params)
   end
   r.array
 end
@@ -154,13 +166,13 @@ function evaluate!(cache,
   x::AbstractVector{<:Point}) where {D,V}
 
   f = fg.fa
-  r, c, g = cache
+  r, s, c, g = cache
   np = length(x)
   _setsize!(f,np,r,c,g)
-  s = zero(Mutable(VectorValue{D,eltype(V)}))
+  params = _get_static_parameters(f)
   for i in 1:np
     @inbounds xi = x[i]
-    _gradient_nd!(f,xi,r,i,c,g,s)
+    _gradient_nd!(f,xi,r,i,c,g,s,params)
   end
   r.array
 end
@@ -170,13 +182,13 @@ function evaluate!(cache,
   x::AbstractVector{<:Point}) where {D,V}
 
   f = fg.fa
-  r, c, g, h = cache
+  r, s, c, g, h = cache
   np = length(x)
   _setsize!(f,np,r,c,g,h)
-  s = zero(Mutable(TensorValue{D,D,eltype(V)}))
+  params = _get_static_parameters(f)
   for i in 1:np
     @inbounds xi = x[i]
-    _hessian_nd!(f,xi,r,i,c,g,h,s)
+    _hessian_nd!(f,xi,r,i,c,g,h,s,params)
   end
   r.array
 end
@@ -234,7 +246,7 @@ end
 ###############################
 
 """
-    _evaluate_nd!(b,xi,r,i,c)
+    _evaluate_nd!(b,xi,r,i,c,params)
 
 Compute and assign: `r`[`i`] = `b`(`xi`) = (`b`₁(`xi`), ..., `b`ₙ(`xi`))
 
@@ -245,13 +257,13 @@ row of `r`.
 function _evaluate_nd!(
   b::PolynomialBasis, xi,
   r::AbstractMatrix, i,
-  c::AbstractMatrix)
+  c::AbstractMatrix, params)
 
   @abstractmethod
 end
 
 """
-    _gradient_nd!(b,xi,r,i,c,g,s)
+    _gradient_nd!(b,xi,r,i,c,g,s,params)
 
 Compute and assign: `r`[`i`] = ∇`b`(`xi`) = (∇`b`₁(`xi`), ..., ∇`b`ₙ(`xi`))
 
@@ -260,19 +272,20 @@ for gradients of `b`ₖ(`xi`), and
 
 - `g` is a mutable `D`×`K` cache (for the 1D polynomials first derivatives).
 - `s` is a mutable length `D` cache for ∇`b`ₖ(`xi`).
+- `params` is an optional (tuple of) parameter(s) returned by [`_get_static_parameters(b)`](@ref _get_static_parameters)
 """
 function _gradient_nd!(
   b::PolynomialBasis, xi,
   r::AbstractMatrix, i,
   c::AbstractMatrix,
   g::AbstractMatrix,
-  s::MVector)
+  s::MVector, params)
 
   @abstractmethod
 end
 
 """
-    _hessian_nd!(b,xi,r,i,c,g,h,s)
+    _hessian_nd!(b,xi,r,i,c,g,h,s,params)
 
 Compute and assign: `r`[`i`] = H`b`(`xi`) = (H`b`₁(`xi`), ..., H`b`ₙ(`xi`))
 
@@ -288,7 +301,7 @@ function _hessian_nd!(
   c::AbstractMatrix,
   g::AbstractMatrix,
   h::AbstractMatrix,
-  s::MMatrix)
+  s::MMatrix, params)
 
   @abstractmethod
 end

src/Polynomials/Polynomials.jl

@@ -13,8 +13,8 @@ or second derivatives.
 Constructors for commonly used bases (see the documentation for the spaces definitions):
 - ℚ spaces: `[Polynomial]Basis(Val(D), V, order)`
 - ℙ spaces: `[Polynomial]Basis(..., Polynomials._p_filter)`
-- ℚₙ\\ℚₙ₋₁: `[Polynomial]Basis(..., Polynomials._qs_filter)`
-- ℙₙ\\ℙₙ₋₁: `[Polynomial]Basis(..., Polynomials._ps_filter)`
+- ℚₙ\\ℚₙ₋₁: `[Polynomial]Basis(..., Polynomials._qh_filter)`
+- ℙₙ\\ℙₙ₋₁: `[Polynomial]Basis(..., Polynomials._ph_filter)`
 - ℕ𝔻(△): [`PGradBasis`](@ref)`(Val(D), T, order)`
 - ℕ𝔻(□): [`QGradBasis`](@ref)`(...)`
 - ℝ𝕋(△): [`PCurlGradBasis`](@ref)`(...)`
@@ -80,13 +80,16 @@ module Polynomials
 
 using DocStringExtensions
 using LinearAlgebra: mul!
+using LinearAlgebra: I
 using StaticArrays
 using Gridap.Helpers
 using Gridap.Arrays
 using Gridap.TensorValues
 using Gridap.Fields
 
 using PolynomialBases: jacobi, jacobi_and_derivative
+using Combinatorics: multiexponents, multinomial, combinations
+using Base.Iterators: take
 
 import Gridap.Fields: evaluate!
 import Gridap.Fields: return_cache
@@ -103,14 +106,18 @@ export Bernstein
 export PolynomialBasis
 export get_order
 
-export UniformPolyBasis
+export CartProdPolyBasis
 export get_exponents
 export get_orders
 export MonomialBasis
 export LegendreBasis
 export ChebyshevBasis
 export BernsteinBasis
 
+export BernsteinBasisOnSimplex
+export bernstein_terms
+export bernstein_term_id
+
 export CompWiseTensorPolyBasis
 export QGradBasis
 export QCurlGradBasis
@@ -132,7 +139,7 @@ export PCurlGradMonomialBasis
 
 include("PolynomialInterfaces.jl")
 
-include("UniformPolyBases.jl")
+include("CartProdPolyBases.jl")
 
 include("CompWiseTensorPolyBases.jl")
 

src/Polynomials/RaviartThomasPolyBases.jl

@@ -1,5 +1,5 @@
 """
-    RaviartThomasPolyBasis{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
+    RaviartThomasPolyBasis{D,V,PT} <: PolynomialBasis{D,V,PT}
 
 Basis of the vector valued (`V<:VectorValue{D}`) space
 
@@ -14,11 +14,12 @@ The space 𝕊ₙ, typically ℙᴰₙ or ℚᴰₙ, does not need to have a ten
 structure of 1D scalar spaces. Thus, the ℝ𝕋ᴰₙ component's scalar spaces are not
 tensor products either.
 
-𝕊ₙ is defined like a scalar valued [`UniformPolyBasis`](@ref) via the `_filter`
+𝕊ₙ is defined like a scalar valued [`CartProdPolyBasis`](@ref) via the `_filter`
 argument of the constructor, by default `_p_filter` for ℙᴰₙ.
 As a consequence, `PT` must be hierarchical, see [`isHierarchical`](@ref).
 """
-struct RaviartThomasPolyBasis{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
+struct RaviartThomasPolyBasis{D,V,PT} <: PolynomialBasis{D,V,PT}
+  max_order::Int
   pterms::Vector{CartesianIndex{D}}
   sterms::Vector{CartesianIndex{D}}
 
@@ -51,36 +52,33 @@ struct RaviartThomasPolyBasis{D,V,K,PT} <: PolynomialBasis{D,V,K,PT}
     _minus_one_order_filter = term -> _filter(Tuple(term) .- indexbase, order-1)
     sterms = filter(!_minus_one_order_filter, pterms)
 
-    new{D,V,order+1,PT}(pterms,sterms)
+    new{D,V,PT}(order+1,pterms,sterms)
   end
 end
 
-Base.size(a::RaviartThomasPolyBasis{D}) where {D} = (D*length(a.pterms) + length(a.sterms), )
+Base.size(b::RaviartThomasPolyBasis{D}) where {D} = (D*length(b.pterms) + length(b.sterms), )
+get_order(b::RaviartThomasPolyBasis) = b.max_order
 
 
 #################################
 # nD evaluations implementation #
 #################################
 
 function _evaluate_nd!(
-  b::RaviartThomasPolyBasis{D,V,K,PT}, x,
+  b::RaviartThomasPolyBasis{D,V,PT}, x,
   r::AbstractMatrix{V}, i,
-  c::AbstractMatrix{T}) where {D,V,K,PT,T}
-
-  pterms = b.pterms
-  sterms = b.sterms
+  c::AbstractMatrix{T}, VK::Val) where {D,V,PT,T}
 
   for d in 1:D
-    Kv = Val(K)
-    _evaluate_1d!(PT,Kv,c,x,d)
+    _evaluate_1d!(PT,VK,c,x,d)
   end
 
   m = zero(Mutable(V))
   k = 1
 
   @inbounds begin
     for l in 1:D
-      for ci in pterms
+      for ci in b.pterms
 
         s = one(T)
         for d in 1:D
@@ -91,7 +89,7 @@ function _evaluate_nd!(
       end
     end
 
-    for ci in sterms
+    for ci in b.sterms
       for i in 1:D
         m[i] = zero(T)
       end
@@ -114,26 +112,22 @@ function _evaluate_nd!(
 end
 
 function _gradient_nd!(
-  b::RaviartThomasPolyBasis{D,V,K,PT}, x,
+  b::RaviartThomasPolyBasis{D,V,PT}, x,
   r::AbstractMatrix{G}, i,
   c::AbstractMatrix{T},
   g::AbstractMatrix{T},
-  s::MVector{D,T}) where {D,V,K,PT,G,T}
-
-  pterms = b.pterms
-  sterms = b.sterms
+  s::MVector{D,T}, VK::Val) where {D,V,PT,G,T}
 
   for d in 1:D
-    Kv = Val(K)
-    _derivatives_1d!(PT,Kv,(c,g),x,d)
+    _derivatives_1d!(PT,VK,(c,g),x,d)
   end
 
   m = zero(Mutable(G))
   k = 1
 
   @inbounds begin
     for l in 1:D
-      for ci in pterms
+      for ci in b.pterms
 
         for i in eachindex(s)
           s[i] = one(T)
@@ -153,7 +147,7 @@ function _gradient_nd!(
       end
     end
 
-    for ci in sterms
+    for ci in b.sterms
 
       for i in eachindex(m)
         m[i] = zero(T)
@@ -215,7 +209,7 @@ b = PCurlGradBasis(Monomial, Val(3), Float64, 2)
 For more details, see [`RaviartThomasPolyBasis`](@ref), as `PCurlGradBasis` returns
 an instance of\\
 `RaviartThomasPolyBasis{D, VectorValue{D,T}, order+1, PT}`  for `D`>1, or\\
-`UniformPolyBasis{1, VectorValue{1,T}, order+1, PT}` for `D`=1.
+`CartProdPolyBasis{1, VectorValue{1,T}, order+1, PT}` for `D`=1.
 """
 function PCurlGradBasis(::Type{PT},::Val{D},::Type{T},order::Int) where {PT,D,T}
   RaviartThomasPolyBasis{D}(PT, T, order)
@@ -225,5 +219,5 @@ function PCurlGradBasis(::Type{PT},::Val{1},::Type{T},order::Int) where {PT,T}
   @check T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
 
   V = VectorValue{1,T}
-  UniformPolyBasis(PT, Val(1), V, order+1)
+  CartProdPolyBasis(PT, Val(1), V, order+1)
 end

src/ReferenceFEs/Dofs.jl

@@ -56,19 +56,19 @@ end
 """
     struct MappedDofBasis{T<:Dof,MT,BT} <: AbstractVector{T}
       F :: MT
-      σ :: BT
+      Σ :: BT
       args
     end
 
-Represents η = σ∘F, evaluated as η(φ) = σ(F(φ,args...))
+Represents { η = σ∘F : σ ∈ Σ }, evaluated as η(φ) = σ(F(φ,args...)) where
 
-  - σ : V* -> R is a dof basis
-  - F : W  -> V is a map between function spaces
+- σ : V -> R  are dofs on V
+- F : W  -> V is a map between function spaces
 
-Intended combinations would be: 
+Intended combinations would be:
 
-- σ : V* -> R dof basis in the physical domain and F* : V̂ -> V is a pushforward map.
-- ̂σ : V̂* -> R dof basis in the reference domain and (F*)^-1 : V -> V̂ is an inverse pushforward map.
+- Σ ⊂ V* a dof basis in the physical domain and F* : V̂ -> V is a pushforward map.
+- ̂Σ ⊂ V̂* a dof basis in the reference domain and (F*)⁻¹ : V -> V̂ is an inverse pushforward map.
 
 """
 struct MappedDofBasis{T<:Dof,MT,BT,A} <: AbstractVector{T}

src/TensorValues/MultiValueTypes.jl

@@ -80,11 +80,11 @@ num_components(T::Type{<:MultiValue}) = @unreachable "$T type is too abstract to
     num_indep_components(::Type{<:Number})
     num_indep_components(a::Number)
 
-Number of independant components of a `Number`, that is `num_components`
+Number of independent components of a `Number`, that is `num_components`
 minus the number of components determined from others by symmetries or constraints.
 
-For example, a `TensorValue{3,3}` has 9 independant components, a `SymTensorValue{3}`
-has 6 and a `SymTracelessTensorValue{3}` has 5. But they all have 9 (non independant) components.
+For example, a `TensorValue{3,3}` has 9 independent components, a `SymTensorValue{3}`
+has 6 and a `SymTracelessTensorValue{3}` has 5. But they all have 9 (non independent) components.
 """
 num_indep_components(::Type{T}) where T<:Number = num_components(T)
 num_indep_components(::T) where T<:Number = num_indep_components(T)
@@ -110,7 +110,7 @@ end
     indep_comp_getindex(a::Number,i)
 
 Get the `i`th independent component of `a`. It only differs from `getindex(a,i)`
-when the components of `a` are interdependant, see [`num_indep_components`](@ref).
+when the components of `a` are interdependent, see [`num_indep_components`](@ref).
 `i` should be in `1:num_indep_components(a)`.
 """
 function indep_comp_getindex(a::Number,i)

test/PolynomialsTests/BernsteinBasesTests.jl

@@ -3,23 +3,25 @@ module BernsteinBasisTests
 using Test
 using Gridap.TensorValues
 using Gridap.Fields
+using Gridap.Arrays
 using Gridap.Polynomials
-using Gridap.Polynomials: binoms
 using ForwardDiff
+using StaticArrays
 
 @test isHierarchical(Bernstein) == false
 
-np = 3
 x = [Point(0.),Point(1.),Point(.4)]
 x1 = x[1]
 
-# Only test 1D evaluations as tensor product structure is tested in monomial tests
-#
+#####################################
+# Tests for 1D Bernstein polynomial #
+#####################################
+
 V = Float64
 G = gradient_type(V,x1)
 H = gradient_type(G,x1)
 
-function test_internals(order,x,bx,∇bg,Hbx)
+function test_internals(order,x,bx,∇bx,Hbx)
   sz = (1,order+1)
   for (i,xi) in enumerate(x)
     v2 = zeros(sz)
@@ -42,7 +44,7 @@ end
 order = 0
 b = BernsteinBasis(Val(1),V,order)
 @test get_order(b) == 0
-@test get_orders(b) == (0,)
+@test get_exponents(b) == ((0,),)
 
 bx =   [ 1.; 1.; 1.;; ]
 ∇bx = G[ 0.; 0.; 0.;; ]
@@ -59,8 +61,8 @@ test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 order = 1
 b = BernsteinBasis(Val(1),V,order)
 
-bx = [ 1.0  0.0
-       0.0  1.0
+bx = [ 1.0  0.
+       0.  1.0
        0.6  0.4]
 
 ∇bx = G[ -1. 1.
@@ -82,8 +84,8 @@ test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 order = 2
 b = BernsteinBasis(Val(1),V,order)
 
-bx  = [ 1.0   0.0   0.0
-        0.0   0.0   1.0
+bx  = [ 1.0   0.    0.
+        0.    0.    1.0
         0.36  0.48  0.16]
 
 
@@ -104,25 +106,26 @@ test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 # Order 3
 
 function bernstein(K,N)
-  b = binoms(Val(K))
+  b = ntuple( i -> binomial(K,i-1), Val(K+1)) # all binomial(i,K) for 0≤i≤K
   t -> b[N+1]*(t^N)*((1-t)^(K-N))
 end
 _∇(b) = t -> ForwardDiff.derivative(b, t)
 _H(b) = t -> ForwardDiff.derivative(y -> ForwardDiff.derivative(b, y), t)
 
+_bx_1D( order,x)   = [      bernstein(order,n)( xi[1])  for xi in x,  n in 0:order]
+_∇bx_1D(order,x,G) = [ G(_∇(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
+_Hbx_1D(order,x,H) = [ H(_H(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
+
 order = 3
 b = BernsteinBasis(Val(1),V,order)
 
 # x=x^1; x2 = x^2; x3 = x^3
 #    -x3+3x2-3x+1  3x3-6x2+3x -3x3+3x2 x3
-bx  = [      bernstein(order,n)( xi[1])  for xi in x,  n in 0:order]
-
+bx  = _bx_1D( order,x)
 #   -3x2+6x-3  9x2-12x+3  -9x2+6x 3x2
-∇bx = [ G(_∇(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
-
+∇bx = _∇bx_1D(order,x,G)
 #    -6x+6 18x-12 -18x+6  6x
-Hbx = [ H(_H(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
-
+Hbx = _Hbx_1D(order,x,H)
 test_internals(order,x,bx,∇bx,Hbx)
 
 test_field_array(b,x,bx,≈, grad=∇bx,gradgrad=Hbx)
@@ -134,14 +137,209 @@ test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 order = 4
 b = BernsteinBasis(Val(1),V,order)
 
-bx  = [      bernstein(order,n)( xi[1])  for xi in x,  n in 0:order]
-∇bx = [ G(_∇(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
-Hbx = [ H(_H(bernstein(order,n))(xi[1])) for xi in x,  n in 0:order]
+bx  = _bx_1D( order,x)
+∇bx = _∇bx_1D(order,x,G)
+Hbx = _Hbx_1D(order,x,H)
 
 test_internals(order,x,bx,∇bx,Hbx)
 
 test_field_array(b,x,bx,≈, grad=∇bx,gradgrad=Hbx)
 test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
 
+#####################################
+# Tests for ND Bernstein polynomial #
+#####################################
+
+function bernstein_nD(D,K,x2λ=nothing)
+  terms = Polynomials.bernstein_terms(K,D)
+  N = length(terms)
+  x -> begin
+    @assert length(x) == D
+    vals = zeros(eltype(x),N)
+    # change to barycentric coords of reference simplex
+    if isnothing(x2λ)
+      λ = zeros(eltype(x), D+1)
+      copyto!(λ, x)
+      λ[D+1] = 1 - sum(x)
+    else
+      λ = x2λ*SVector(1, x...)
+    end
+    for i in 1:N
+      vals[i] = Polynomials.multinomial(terms[i]...)
+      for (λi,ei) in zip(λ,terms[i])
+        vals[i] *= λi^ei
+      end
+    end
+    return vals
+  end
+end
+
+_∇(b) = x -> ForwardDiff.jacobian(b, get_array(x))
+_H(b) = x -> ForwardDiff.jacobian(y -> ForwardDiff.jacobian(b, y), get_array(x))
+
+_bx( D,order,x,  x2λ=nothing) = transpose(reduce(hcat, (                                    bernstein_nD(D,order,x2λ )(xi)           for xi in x)))
+_∇bx(D,order,x,G,x2λ=nothing) = transpose(reduce(hcat, ( map(G,          eachrow(        _∇(bernstein_nD(D,order,x2λ))(xi)))         for xi in x)))
+_Hbx(D,order,x,H,x2λ=nothing) = transpose(reduce(hcat, ( map(x->H(x...), eachrow(reshape(_H(bernstein_nD(D,order,x2λ))(xi), :,D*D))) for xi in x)))
+
+D = 2
+x = [Point(1.,0.), Point(.0,.5), Point(1.,.5), Point(.2,.3)]
+x3 = x[3]
+
+# scalar valued in 2D
+V = Float64
+G = gradient_type(V,x3)
+H = gradient_type(G,x3)
+
+order = 0
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+@test get_order(b) == 0
+bx  = _bx( D,order,x)
+∇bx = _∇bx(D,order,x,G)
+Hbx = _Hbx(D,order,x,H)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+
+order = 1
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+bx  = _bx( D,order,x)
+∇bx = _∇bx(D,order,x,G)
+Hbx = _Hbx(D,order,x,H)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+
+order = 2
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+bx  = _bx( D,order,x)
+∇bx = _∇bx(D,order,x,G)
+Hbx = _Hbx(D,order,x,H)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+
+order = 3
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+bx  = _bx( D,order,x)
+∇bx = _∇bx(D,order,x,G)
+Hbx = _Hbx(D,order,x,H)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+
+order = 4
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+bx  = _bx( D,order,x)
+∇bx = _∇bx(D,order,x,G)
+Hbx = _Hbx(D,order,x,H)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+test_field_array(b,x[1],bx[1,:],≈,grad=∇bx[1,:],gradgrad=Hbx[1,:])
+
+b_terms = bernstein_terms(order,D)
+λ = Polynomials._cart_to_bary(x3, b.cart_to_bary_matrix)
+for j in 1:length(b)
+  α = b_terms[j]
+  id_α = bernstein_term_id(α)
+  @test id_α == j
+  c = Float64[ Int(i==id_α) for i in 1:length(b) ] # Bernstein coefficients of Bα
+  Polynomials._de_Casteljau_nD!(c, λ, Val(order), Val(D))
+  Bα_x3 = c[1]
+  @test Bα_x3 == bx[3,id_α]
+end
+
+# Vector valued in 2D
+V = VectorValue{3,Float64}
+G = gradient_type(V,x3)
+H = gradient_type(G,x3)
+
+order = 2
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+bx  = V[(1., 0., 0.) (0., 1., 0.) (0., 0., 1.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.);
+        (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0., 0., 0.) (0.25, 0., 0.) (0., 0.25, 0.) (0., 0., 0.25) (0.5, 0., 0.) (0., 0.5, 0.) (0., 0., 0.5) (0.25, 0., 0.) (0., 0.25, 0.) (0., 0., 0.25);
+        (1., 0., 0.) (0., 1., 0.) (0., 0., 1.) (1., 0., 0.) (0., 1., 0.) (0., 0., 1.) (-1., 0., 0.) (0., -1., 0.) (0., 0., -1.) (0.25, 0., 0.) (0., 0.25, 0.) (0., 0., 0.25) (-0.5, 0., 0.) (0., -0.5, 0.) (0., 0., -0.5) (0.25, 0., 0.) (0., 0.25, 0.) (0., 0., 0.25);
+        (0.04, 0., 0.) (0., 0.04, 0.) (0., 0., 0.04) (0.12, 0., 0.) (0., 0.12, 0.) (0., 0., 0.12) (0.2, 0., 0.) (0., 0.2, 0.) (0., 0., 0.2) (0.09, 0., 0.) (0., 0.09, 0.) (0., 0., 0.09) (0.3, 0., 0.) (0., 0.3, 0.) (0., 0., 0.3) (0.25, 0., 0.) (0., 0.25, 0.) (0., 0., 0.25)]
+∇bx = G[(2., 0., 0., 0., 0., 0.) (0., 0., 2., 0., 0., 0.) (0., 0., 0., 0., 2., 0.) (0., 2., 0., 0., 0., 0.) (0., 0., 0., 2., 0., 0.) (0., 0., 0., 0., 0., 2.) (-2., -2., 0., 0., 0., 0.) (0., 0., -2., -2., 0., 0.) (0., 0., 0., 0., -2., -2.) (0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0.) (-0., 0., 0., 0., 0., 0.) (0., 0., -0., 0., 0., 0.) (0., 0., 0., 0., -0., 0.) (0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0.);
+        (0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0.) (1., 0., 0., 0., 0., 0.) (0., 0., 1., 0., 0., 0.) (0., 0., 0., 0., 1., 0.) (1., -0., 0., 0., 0., 0.) (0., 0., 1., -0., 0., 0.) (0., 0., 0., 0., 1., -0.) (0., 1., 0., 0., 0., 0.) (0., 0., 0., 1., 0., 0.) (0., 0., 0., 0., 0., 1.) (-1., 0., 0., 0., 0., 0.) (0., 0., -1., 0., 0., 0.) (0., 0., 0., 0., -1., 0.) (-1., -1., 0., 0., 0., 0.) (0., 0., -1., -1., 0., 0.) (0., 0., 0., 0., -1., -1.);
+        (2., 0., 0., 0., 0., 0.) (0., 0., 2., 0., 0., 0.) (0., 0., 0., 0., 2., 0.) (1., 2., 0., 0., 0., 0.) (0., 0., 1., 2., 0., 0.) (0., 0., 0., 0., 1., 2.) (-3., -2., 0., 0., 0., 0.) (0., 0., -3., -2., 0., 0.) (0., 0., 0., 0., -3., -2.) (0., 1., 0., 0., 0., 0.) (0., 0., 0., 1., 0., 0.) (0., 0., 0., 0., 0., 1.) (-1., -2., 0., 0., 0., 0.) (0., 0., -1., -2., 0., 0.) (0., 0., 0., 0., -1., -2.) (1., 1., 0., 0., 0., 0.) (0., 0., 1., 1., 0., 0.) (0., 0., 0., 0., 1., 1.);
+        (0.4, 0., 0., 0., 0., 0.) (0., 0., 0.4, 0., 0., 0.) (0., 0., 0., 0., 0.4, 0.) (0.6, 0.4, 0., 0., 0., 0.) (0., 0., 0.6, 0.4, 0., 0.) (0., 0., 0., 0., 0.6, 0.4) (0.6, -0.4, 0., 0., 0., 0.) (0., 0., 0.6, -0.4, 0., 0.) (0., 0., 0., 0., 0.6, -0.4) (0., 0.6, 0., 0., 0., 0.) (0., 0., 0., 0.6, 0., 0.) (0., 0., 0., 0., 0., 0.6) (-0.6, 0.4, 0., 0., 0., 0.) (0., 0., -0.6, 0.4, 0., 0.) (0., 0., 0., 0., -0.6, 0.4) (-1., -1., 0., 0., 0., 0.) (0., 0., -1., -1., 0., 0.) (0., 0., 0., 0., -1., -1.)]
+Hbx = H[(2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0.) (0., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 2., 2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 0.) (-4., -2., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., -4., -2., -2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., -4., -2., -2., 0.) (0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.) (0., -2., -2., -4., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., -2., -2., -4., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., -2., -2., -4.) (2., 2., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 2., 2., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2.);
+        (2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0.) (0., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 2., 2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 0.) (-4., -2., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., -4., -2., -2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., -4., -2., -2., 0.) (0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.) (0., -2., -2., -4., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., -2., -2., -4., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., -2., -2., -4.) (2., 2., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 2., 2., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2.);
+        (2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0.) (0., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 2., 2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 0.) (-4., -2., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., -4., -2., -2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., -4., -2., -2., 0.) (0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.) (0., -2., -2., -4., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., -2., -2., -4., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., -2., -2., -4.) (2., 2., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 2., 2., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2.);
+        (2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0.) (0., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 2., 2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 0.) (-4., -2., -2., 0., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., -4., -2., -2., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., -4., -2., -2., 0.) (0., 0., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.) (0., -2., -2., -4., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 0., -2., -2., -4., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 0., -2., -2., -4.) (2., 2., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0.) (0., 0., 0., 0., 2., 2., 2., 2., 0., 0., 0., 0.) (0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2.)]
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
+
+# scalar valued in 3D
+
+x = [Point(0.,0.,1.), Point(.5,.5,.5), Point(1.,.2,.4), Point(.2,.4,.3)]
+x1 = x[1]
+
+V = Float64
+G = gradient_type(V,x1)
+H = gradient_type(G,x1)
+
+order = 4
+D = 3
+b = BernsteinBasisOnSimplex(Val(D),V,order)
+bx  = _bx( D,order,x)
+∇bx = _∇bx(D,order,x,G)
+Hbx = _Hbx(D,order,x,H)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
+
+
+############################################################################
+# Tests for ND Bernstein polynomial with arbitrary barycentric coordinates #
+############################################################################
+
+D = 2
+T = Float64
+
+vertices = (Point(0.,1.), Point(0.,2.), Point(0.,3.))
+@test_throws DomainError Polynomials._compute_cart_to_bary_matrix(vertices, Val(D+1))
+
+vertices = (Point(5.,0.), Point(7.,2.), Point(0.,3.))
+
+b = BernsteinBasisOnSimplex(Val(2), Float64, 3, vertices)
+x = [Point(.0,.5), Point(1.,.5), Point(.2,.3), Point(5.,0.), Point(7.,2.), Point(0.,3.)]
+x1 = x[1]
+
+for xi in x
+  λi = Polynomials._cart_to_bary(xi, b.cart_to_bary_matrix)
+  @test sum(λi) ≈ 1.
+  @test xi ≈ sum(λi .* vertices)
+end
+
+
+# Scalar value in 2D
+D = 2
+vertices = (Point(5.,0.), Point(7.,2.), Point(0.,3.))
+x2λ = Polynomials._compute_cart_to_bary_matrix(vertices, Val(D+1))
+
+T = Float64
+G = gradient_type(V,x1)
+H = gradient_type(G,x1)
+
+order = 4
+b = BernsteinBasisOnSimplex(Val(D),V,order,vertices)
+bx  = _bx( D,order,x,  x2λ)
+∇bx = _∇bx(D,order,x,G,x2λ)
+Hbx = _Hbx(D,order,x,H,x2λ)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+test_field_array(b,x1,bx[1,:],≈,grad=∇bx[1,:],gradgrad=Hbx[1,:])
+
+
+# scalar valued in 3D
+D = 3
+vertices = (Point(5.,0.,0.), Point(7.,0.,2.), Point(0.,3.,3.), Point(3.,0.,3.))
+x2λ = Polynomials._compute_cart_to_bary_matrix(vertices, Val(D+1))
+
+x = [Point(0.,0.,1.), Point(.5,.5,.5), Point(1.,.2,.4), Point(.2,.4,.3)]
+x1 = x[1]
+
+V = Float64
+G = gradient_type(V,x1)
+H = gradient_type(G,x1)
+
+order = 4
+b = BernsteinBasisOnSimplex(Val(D),V,order,vertices)
+bx  = _bx( D,order,x,  x2λ)
+∇bx = _∇bx(D,order,x,G,x2λ)
+Hbx = _Hbx(D,order,x,H,x2λ)
+test_field_array(b,x,bx,≈, grad=∇bx, gradgrad=Hbx)
+test_field_array(b,x1,bx[1,:],≈,grad=∇bx[1,:],gradgrad=Hbx[1,:])
+
 end # module

test/PolynomialsTests/ChebyshevBasesTests.jl

@@ -4,7 +4,6 @@ using Test
 using Gridap.TensorValues
 using Gridap.Fields
 using Gridap.Polynomials
-using Gridap.Polynomials: binoms
 using ForwardDiff
 
 @test isHierarchical(Chebyshev) == true

test/PolynomialsTests/ModalC0BasesTests.jl

@@ -39,13 +39,13 @@ b1 = ModalC0Basis{1}(V,order,a,b)
                                      (0.0,) (0.0,) (3.4641016151377544,) (-1.4907119849998598,);
                                      (0.0,) (0.0,) (3.4641016151377544,) (2.9814239699997196,)]
 
-# Validate generic 1D implem using UniformPolyBasis
+# Validate generic 1D implem using CartProdPolyBasis
 
 order = 3
 a = fill(Point(0.),order+1)
 b = fill(Point(1.),order+1)
 b1 = ModalC0Basis{1}(V,order,a,b)
-b1u= UniformPolyBasis(ModalC0,Val(1),V,order)
+b1u= CartProdPolyBasis(ModalC0,Val(1),V,order)
 
 ∇b1  = Broadcasting(∇)(b1)
 ∇b1u = Broadcasting(∇)(b1u)
@@ -83,15 +83,15 @@ H = gradient_type(G,x1)
                                             (0.0, 0.5590169943749475, 0.5590169943749475, 1.118033988749895);
                                             (0.0, -2.23606797749979, -2.23606797749979, 0.0) ]
 
-# Validate generic 2D implem using UniformPolyBasis
+# Validate generic 2D implem using CartProdPolyBasis
 
 order = 3
 len_b2 = (order+1)^2
 a = fill(Point(0.,0.), len_b2)
 b = fill(Point(1.,1.), len_b2)
 
 b2 = ModalC0Basis{2}(V,order,a,b)
-b2u= UniformPolyBasis(ModalC0,Val(2),V,order)
+b2u= CartProdPolyBasis(ModalC0,Val(2),V,order)
 ∇b2  = Broadcasting(∇)(b2)
 ∇b2u = Broadcasting(∇)(b2u)
 

test/PolynomialsTests/MonomialBasesTests.jl

@@ -5,7 +5,7 @@ using Gridap.TensorValues
 using Gridap.Fields
 using Gridap.Polynomials
 
-using Gridap.Polynomials: _q_filter, _qs_filter, _p_filter, _ps_filter
+using Gridap.Polynomials: _q_filter, _qh_filter, _p_filter, _ph_filter
 
 @test isHierarchical(Monomial) == true
 
@@ -53,18 +53,21 @@ test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
 # Real-valued Q space with an isotropic order
 
-orders = (1,2)
+orders = (1,3)
 V = Float64
 G = gradient_type(V,xi)
+H = gradient_type(G,xi)
 b = MonomialBasis(Val(2),V,orders)
 
-v = V[1.0, 2.0, 3.0, 6.0, 9.0, 18.0]
-g = G[(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (3.0, 2.0), (0.0, 6.0), (9.0, 12.0)]
+v = V[1.0, 2.0, 3.0, 6.0, 9.0, 18.0, 27.0, 54.0]
+g = G[(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (3.0, 2.0), (0.0, 6.0), (9.0, 12.0), (0., 27.0), (27.0, 54.0)]
+h = H[(0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (0.0, 1.0, 1.0, 0.0), (0.0, 0.0, 0.0, 2.0), (0.0, 6.0, 6.0, 4.0), (0.0, 0.0, 0.0, 18.0), (0.0, 27.0, 27.0, 36.0)]
 
 bx = repeat(permutedims(v),np)
 ∇bx = repeat(permutedims(g),np)
-test_field_array(b,x,bx,grad=∇bx)
-test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:])
+Hbx = repeat(permutedims(h),np)
+test_field_array(b,x,bx,grad=∇bx,gradgrad=Hbx)
+test_field_array(b,x[1],bx[1,:],grad=∇bx[1,:],gradgrad=Hbx[1,:])
 
 # Vector-valued Q space with isotropic order
 
@@ -254,11 +257,11 @@ order = 2
 @test _q_filter( (1,1) ,order) == true
 @test _q_filter( (3,1) ,order) == false
 
-@test _qs_filter( (1,2) ,order) == true
-@test _qs_filter( (2,0) ,order) == true
-@test _qs_filter( (2,2) ,order) == true
-@test _qs_filter( (1,1) ,order) == false
-@test _qs_filter( (3,1) ,order) == false
+@test _qh_filter( (1,2) ,order) == true
+@test _qh_filter( (2,0) ,order) == true
+@test _qh_filter( (2,2) ,order) == true
+@test _qh_filter( (1,1) ,order) == false
+@test _qh_filter( (3,1) ,order) == false
 
 @test _p_filter( (1,2) ,order) == false
 @test _p_filter( (2,0) ,order) == true
@@ -267,10 +270,10 @@ order = 2
 @test _p_filter( (3,1) ,order) == false
 @test _p_filter( (0,1) ,order) == true
 
-@test _ps_filter( (1,2) ,order) == false
-@test _ps_filter( (2,0) ,order) == true
-@test _ps_filter( (2,2) ,order) == false
-@test _ps_filter( (1,1) ,order) == true
-@test _ps_filter( (3,1) ,order) == false
+@test _ph_filter( (1,2) ,order) == false
+@test _ph_filter( (2,0) ,order) == true
+@test _ph_filter( (2,2) ,order) == false
+@test _ph_filter( (1,1) ,order) == true
+@test _ph_filter( (3,1) ,order) == false
 
 end # module

test/PolynomialsTests/PCurlGradBasesTests.jl

@@ -73,7 +73,7 @@ T = Float64
 V = VectorValue{D,T}
 b = PCurlGradBasis(PT,Val(D),T,order)
 
-@test b isa UniformPolyBasis{D,V,order+1,PT}
+@test b isa CartProdPolyBasis{D,V,PT}
 
 @test_throws AssertionError PCurlGradBasis(PT,Val(D),V,order)
 

test/PolynomialsTests/PGradBasesTests.jl

@@ -60,7 +60,7 @@ T = Float64
 V = VectorValue{D,T}
 b = PGradBasis(PT,Val(D),T,order)
 
-@test b isa UniformPolyBasis{D,V,order+1,PT}
+@test b isa CartProdPolyBasis{D,V,PT}
 
 @test_throws AssertionError PGradBasis(PT,Val(D),V,order)
 

test/PolynomialsTests/PolynomialInterfacesTests.jl

@@ -1,4 +1,4 @@
-module ChangeBasisTests
+module PolynomialInterfacesTests
 
 using Test
 using Gridap.Fields
@@ -42,17 +42,17 @@ Polynomials._derivatives_1d!(MockPolynomial, Val(1), (c,), xi, 1)
 
 T = Float64
 D = 1
-struct MockPolyBasis <: PolynomialBasis{D,T,0,MockPolynomial} end
+struct MockPolyBasis <: PolynomialBasis{D,T,MockPolynomial} end
 
 mb = MockPolyBasis()
 
-
 # Implemented interfaces
 @test IndexStyle(mb) == IndexLinear()
 @test return_type(mb) == T
-@test get_order(mb) == 0
 @test mb[1] == MockPolynomial()
+@test_throws ErrorException get_order(mb)
 
+Polynomials.get_order(b::MockPolyBasis) = 0
 
 # Interfaces to implement
 @test_throws ErrorException size(mb)
@@ -61,21 +61,16 @@ Base.size(::MockPolyBasis) = (1,)
 @test length(mb) == 1
 
 
-r, c = return_cache(mb,x)
-
-@test_throws ErrorException Polynomials._evaluate_nd!(mb, xi, r, 1, c)
+r, _, c = return_cache(mb,x)
+@test_throws ErrorException Polynomials._evaluate_nd!(mb, xi, r, 1, c, nothing)
 
 ∇mb = FieldGradientArray{1}(mb)
-r, c, g = return_cache(∇mb,x)
-s = MVector{D,T}(0.)
-
-@test_throws ErrorException Polynomials._gradient_nd!(mb, xi, r, 1, c, g, s)
+r, s, c, g = return_cache(∇mb,x)
+@test_throws ErrorException Polynomials._gradient_nd!(mb, xi, r, 1, c, g, s, nothing)
 
 Hmb = FieldGradientArray{2}(mb)
-r, c, g, h = return_cache(Hmb,x)
-s = MMatrix{D,D,T}(0.)
-
-@test_throws ErrorException Polynomials._hessian_nd!(mb, xi, r, 1, c, g, h, s)
+r, s, c, g, h = return_cache(Hmb,x)
+@test_throws ErrorException Polynomials._hessian_nd!(mb, xi, r, 1, c, g, h, s, nothing)
 
 
 end

test/PolynomialsTests/QCurlGradBasesTests.jl

@@ -94,7 +94,7 @@ T = Float64
 V = VectorValue{D,T}
 b = QCurlGradBasis(PT,Val(D),T,order)
 
-@test b isa UniformPolyBasis{D,V,order+1,PT}
+@test b isa CartProdPolyBasis{D,V,PT}
 
 @test_throws AssertionError QCurlGradBasis(PT,Val(D),V,order)
 

test/PolynomialsTests/QGradBasesTests.jl

@@ -148,7 +148,7 @@ T = Float64
 V = VectorValue{D,T}
 b = QGradBasis(PT,Val(D),T,order)
 
-@test b isa UniformPolyBasis{D,V,order+1,PT}
+@test b isa CartProdPolyBasis{D,V,PT}
 
 @test_throws AssertionError QGradBasis(PT,Val(D),V,order)
 

test/ReferenceFEsTests/CLagrangianRefFEsTests.jl

@@ -12,22 +12,22 @@ using Gridap.ReferenceFEs: monomial_basis
 
 orders = (2,3)
 b = monomial_basis(Float64,QUAD,orders)
-r = [(0,0), (1,0), (2,0), (0,1), (1,1), (2,1), (0,2), (1,2), (2,2), (0,3), (1,3), (2,3)]
+r = ((0,0), (1,0), (2,0), (0,1), (1,1), (2,1), (0,2), (1,2), (2,2), (0,3), (1,3), (2,3))
 @test get_exponents(b) == r
 
 orders = (1,1,2)
 b = monomial_basis(Float64,WEDGE,orders)
-r = [(0,0,0), (1,0,0), (0,1,0), (0,0,1), (1,0,1), (0,1,1), (0,0,2), (1,0,2), (0,1,2)]
+r = ((0,0,0), (1,0,0), (0,1,0), (0,0,1), (1,0,1), (0,1,1), (0,0,2), (1,0,2), (0,1,2))
 @test get_exponents(b) == r
 
 orders = (1,1,1)
 b = monomial_basis(Float64,PYRAMID,orders)
-r = [(0,0,0), (1,0,0), (0,1,0), (1,1, 0), (0,0,1)]
+r = ((0,0,0), (1,0,0), (0,1,0), (1,1, 0), (0,0,1))
 @test get_exponents(b) == r
 
 orders = (1,1,1)
 b = monomial_basis(Float64,TET,orders)
-r = [(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
+r = ((0,0,0), (1,0,0), (0,1,0), (0,0,1))
 @test get_exponents(b) == r
 
 orders = (2,2)