diff --git a/NEWS.md b/NEWS.md index be4ea7e64..d8788309c 100644 --- a/NEWS.md +++ b/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` diff --git a/docs/generate_diagrams.jl b/docs/generate_diagrams.jl index ba0e7fc3f..1916cfebf 100644 --- a/docs/generate_diagrams.jl +++ b/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 diff --git a/docs/make.jl b/docs/make.jl index eaa82166a..1dd1e3678 100644 --- a/docs/make.jl +++ b/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", ], ] diff --git a/docs/src/Polynomials.md b/docs/src/Polynomials.md index c1f116dea..d65fcf5b6 100644 --- a/docs/src/Polynomials.md +++ b/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 diff --git a/docs/src/assets/poly_2.svg b/docs/src/assets/poly_2.svg index 0299b85dc..0973a1c45 100644 --- a/docs/src/assets/poly_2.svg +++ b/docs/src/assets/poly_2.svg @@ -1 +1 @@ -AbstractArray{<:Polynomial}PolynomialBasisget_orderreturn_typeUniformPolyBasisget_exponentsget_ordersCompWiseTensorPolyBasisRaviartThomasPolyBasisNedelecPolyBasisOnSimplexModalC0Basisget_orders(<:Polynomial)BasisQGrad[<:Polynomial]BasisQCurlGrad[<:Polynomial]BasisPCurlGrad[<:Polynomial]BasisPGradMonomialBasis \ No newline at end of file +AbstractArray{<:Polynomial}PolynomialBasisget_orderreturn_typeCartProdPolyBasisget_exponentsget_ordersCompWiseTensorPolyBasisRaviartThomasPolyBasisNedelecPolyBasisOnSimplexBernsteinBasisOnSimplexModalC0Basisget_orders(<:Polynomial)BasisQGrad[<:Polynomial]BasisQCurlGrad[<:Polynomial]BasisPCurlGrad[<:Polynomial]BasisPGradMonomialBasis diff --git a/docs/src/dev-notes/bernstein.md b/docs/src/dev-notes/bernstein.md new file mode 100644 index 000000000..7f62746bb --- /dev/null +++ b/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) diff --git a/src/Polynomials/BernsteinBases.jl b/src/Polynomials/BernsteinBases.jl index 655f2b988..d8c2a3cd9 100644 --- a/src/Polynomials/BernsteinBases.jl +++ b/src/Polynomials/BernsteinBases.jl @@ -7,12 +7,18 @@ struct Bernstein <: Polynomial end isHierarchical(::Type{Bernstein}) = false + +##################################### +# Cartesian product Bernstein bases # +##################################### + """ - BernsteinBasis{D,V,K} = UniformPolyBasis{D,V,K,Bernstein} + BernsteinBasis{D,V} = CartProdPolyBasis{D,V,Bernstein} -Alias for Bernstein multivariate scalar' or `Multivalue`'d basis. +Alias for cartesian product of a scalar tensor 1D-Bernstein basis, scalar valued +or multivalued. """ -const BernsteinBasis{D,V,K} = UniformPolyBasis{D,V,K,Bernstein} +const BernsteinBasis{D,V} = CartProdPolyBasis{D,V,Bernstein} """ BernsteinBasis(::Val{D}, ::Type{V}, order::Int, terms::Vector) @@ -21,24 +27,18 @@ const BernsteinBasis{D,V,K} = UniformPolyBasis{D,V,K,Bernstein} High level constructors of [`BernsteinBasis`](@ref). """ -BernsteinBasis(args...) = UniformPolyBasis(Bernstein, args...) - - -# 1D evaluation implementation - -""" - binoms(::Val{K}) +BernsteinBasis(args...) = CartProdPolyBasis(Bernstein, args...) -Returns the tuple of binomials ( C₍ₖ₀₎, C₍ₖ₁₎, ..., C₍ₖₖ₎ ). -""" -binoms(::Val{K}) where K = ntuple( i -> binomial(K,i-1), Val(K+1)) +################################ +# 1D evaluation implementation # +################################ function _evaluate_1d!(::Type{Bernstein},::Val{0},v::AbstractMatrix{T},x,d) where {T<:Number} @inbounds v[d,1] = one(T) end -@inline function _De_Casteljau_step_1D!(v,d,i,λ1,λ2) +@inline function _de_Casteljau_step_1D!(v,d,i,λ1,λ2) # i = k+1 # vₖ <- xvₖ₋₁ # Bᵏₖ(x) = x*Bᵏ⁻¹ₖ₋₁(x) @@ -64,19 +64,11 @@ function _evaluate_1d!(::Type{Bernstein},::Val{K},v::AbstractMatrix{T},x,d) wher v[d,2] = λ2 for i in 3:n - _De_Casteljau_step_1D!(v,d,i,λ1,λ2) - ## vₖ <- xvₖ₋₁ # Bᵏₖ(x) = x*Bᵏ⁻¹ₖ₋₁(x) - #v[d,i] = λ2*v[d,i-1] - ## vⱼ <- xvⱼ₋₁ + (1-x)vⱼ # Bᵏⱼ(x) = x*Bᵏ⁻¹ⱼ₋₁(x) + (1-x)*Bᵏ⁻¹ⱼ(x) for j = k-1, k-2, ..., 1 - #for l in i-1:-1:2 - # v[d,l] = λ2*v[d,l-1] + λ1*v[d,l] - #end - ## v₀ <- (1-x)v₀ # Bᵏ₀(x) = (1-x)*Bᵏ⁻¹₀(x) - #v[d,1] = λ1*v[d,1] + _de_Casteljau_step_1D!(v,d,i,λ1,λ2) end end # still optimisable for K > 2/3: - # - compute bj = binoms(k,j) at compile time (binoms(Val(K)) function) + # - compute bj = binomials(k,j) ∀j at compile time # - compute vj = xʲ*(1-x)ᴷ⁻ʲ recursively in place like De Casteljau (saving half the redundant multiplications) # - do it in a stack allocated cache (MVector, Bumber.jl) # - @simd affect bj * vj in v[d,i] for all j @@ -94,7 +86,7 @@ end # First derivative of the jth Bernstein poly of order K at x: # (Bᵏⱼ)'(x) = K * ( Bᵏ⁻¹ⱼ₋₁(x) - Bᵏ⁻¹ⱼ(x) ) # = K * x^(j-1) * (1-x)^(K-j-1) * ((1-x)*binom(K-1,j-1) - x*binom(K-1,j)) -function _gradient_1d!(::Type{Bernstein},::Val{K}, g::AbstractMatrix{T},x,d) where {K,T<:Number} +function _gradient_1d!(::Type{Bernstein},::Val{K},g::AbstractMatrix{T},x,d) where {K,T<:Number} @inbounds begin n = K + 1 # n > 2 @@ -181,7 +173,7 @@ function _derivatives_1d!(::Type{Bernstein},::Val{K},t::NTuple{2},x,d) where K g[d,1] = -K*v[d,1] # Last step of De Casteljau for _evaluate_1d! - _De_Casteljau_step_1D!(v,d,n,λ1,λ2) + _de_Casteljau_step_1D!(v,d,n,λ1,λ2) end end @@ -196,14 +188,14 @@ function _derivatives_1d!(::Type{Bernstein},::Val{K},t::NTuple{3},x,d) where K n = K + 1 # n > 3 v, g, h = t - KK = K*(K-1) λ2 = x[d] λ1 = one(eltype(v)) - λ2 # De Casteljau until Bᵏ⁻²ⱼ ∀j _evaluate_1d!(Bernstein,Val(K-2),v,x,d) - # Compute hessians as _hessian_1d! + # Compute hessians as in _hessian_1d! + KK = K*(K-1) h[d,n] = KK*v[d,n-2] h[d,n-1] = KK*( v[d,n-3] -2*v[d,n-2] ) @simd for l in n-2:-1:3 @@ -213,9 +205,9 @@ function _derivatives_1d!(::Type{Bernstein},::Val{K},t::NTuple{3},x,d) where K h[d,1] = KK*v[d,1] # One step of De Casteljau to get Bᵏ⁻¹ⱼ ∀j - _De_Casteljau_step_1D!(v,d,n-1,λ1,λ2) + _de_Casteljau_step_1D!(v,d,n-1,λ1,λ2) - # Compute gradients as _gradient_1d! + # Compute gradients as in _gradient_1d! g[d,n] = K*v[d,n-1] @simd for l in n-1:-1:2 g[d,l] = K*(v[d,l-1] - v[d,l]) @@ -223,6 +215,668 @@ function _derivatives_1d!(::Type{Bernstein},::Val{K},t::NTuple{3},x,d) where K g[d,1] = -K*v[d,1] # Last step of De Casteljau for _evaluate_1d! - _De_Casteljau_step_1D!(v,d,n,λ1,λ2) + _de_Casteljau_step_1D!(v,d,n,λ1,λ2) + end +end + + +################################################### +# Bernstein bases on simplices using de Casteljau # +################################################### + +""" + BernsteinBasisOnSimplex{D,V,M} <: PolynomialBasis{D,V,Bernstein} + +Type for the multivariate Bernstein basis in barycentric coordinates, +c.f. [Bernstein polynomials](@ref) section of the documentation. If `V` is not +scalar, a Cartesian product of the Bernstein scalar basis is made for each +independent component of `V`. + +The index of `B_α` in the basis is [`bernstein_term_id(α)`](@ref bernstein_term_id). + +`M` is Nothing for the reference tetrahedra barycentric coordinates or +`SMatrix{D+1,D+1}` if some simplex (triangle, tetrahedra, ...) vertices +coordinates are given. +""" +struct BernsteinBasisOnSimplex{D,V,M} <: PolynomialBasis{D,V,Bernstein} + max_order::Int + cart_to_bary_matrix::M # Nothing or SMatrix{D+1,D+1} + + function BernsteinBasisOnSimplex{D}(::Type{V},order::Int,vertices=nothing) where {D,V} + msg = "A D simplex defined by D+1 (linearly independent) vertices of type <:Point{D} is required" + @check (isnothing(vertices) || length(vertices) == D+1 && eltype(vertices) <: Point{D}) msg + + K = Int(order) + cart_to_bary_matrix = _compute_cart_to_bary_matrix(vertices, Val(D+1)) + M = typeof(cart_to_bary_matrix) # Nothing or SMatrix + new{D,V,M}(K,cart_to_bary_matrix) + end +end + +""" + BernsteinBasisOnSimplex(::Val{D},::Type{V},order::Int) + BernsteinBasisOnSimplex(::Val{D},::Type{V},order::Int,vertices) + +Constructor for [`BernsteinBasisOnSimplex`](@ref). + +If specified, `vertices` is a collection of `D+1` `Point{D}` defining a simplex +used to compute the barycentric coordinates from, it must be non-degenerated +(have nonzero volume). +""" +function BernsteinBasisOnSimplex(::Val{D},::Type{V},order::Int,vertices=nothing) where {D,V} + BernsteinBasisOnSimplex{D}(V,order,vertices) +end + +Base.size(b::BernsteinBasisOnSimplex{D,V}) where {D,V} = (num_indep_components(V)*binomial(D+get_order(b),D),) +get_order(b::BernsteinBasisOnSimplex) = b.max_order + + +##################### +# Bernstein Helpers # +##################### + +""" + _compute_cart_to_bary_matrix(vertices, ::Val{N}) + _compute_cart_to_bary_matrix(::Nothing,::Val) = nothing + +For the given the vertices of a `D`-simplex (`D` = `N`-1, computes the change +of coordinate matrix `x_to_λ` from cartesian to barycentric, that is +`λ` = `x_to_λ` * `x` such that `sum(λ) == 1` and `x == sum(λ .* vertices)`. +""" +function _compute_cart_to_bary_matrix(vertices, ::Val{N}) where N + T = eltype(eltype(vertices)) + λ_to_x = MMatrix{N,N,T}(undef) + for (i,v) in enumerate(vertices) + λ_to_x[:,i] .= tuple(one(T), v...) + end + + x_to_λ = inv(λ_to_x) # Doesn't throw if singular because this is a StaticArrays Matrix + msg = "The simplex defined by the given vertices is degenerated (is flat / has zero volume)." + !all(isfinite, x_to_λ) && throw(DomainError(vertices,msg)) + + return SMatrix{N,N,T}(x_to_λ) +end +_compute_cart_to_bary_matrix(::Nothing, ::Val) = nothing + +""" + _cart_to_bary(x::Point{D,T}, ::Nothing) + +Converts the cartesian coordinates `x` into the barycentric coordinates with +respect to the reference simplex, that is `λ`=(x1, ..., xD, 1-x1-x2-...-xD). +""" +@inline function _cart_to_bary(x::Point{D,T}, ::Nothing) where {D,T} + return SVector(x..., 1-sum(x)) +end + +""" + _cart_to_bary(x::Point{D,T}, x_to_λ) + +Converts the cartesian coordinates `x` into the barycentric coordinates using +the `x_to_λ` change of coordinate matrix, see [`_compute_cart_to_bary_matrix`](@ref). +""" +@inline function _cart_to_bary(x::Point{D,T}, x_to_λ) where {D,T} + x_1 = SVector{D+1,T}(one(T), x...) + return x_to_λ*x_1 +end + +""" + bernstein_terms(K,D) + +Return the vector of multi-indices for the `D`-dimensional Bernstein basis of +order `K`, that is + +{ α ∈ {0:K}ᴰ⁺¹ | |α| = K } + +ordered in decreasing lexicographic order, e.g. {200, 110, 101, 020, 011, 002} for K=2, D=2. +""" +function bernstein_terms(K,D) + terms = collect(multiexponents(D+1,K)) + terms = convert(Vector{Vector{Int}}, terms) +end + + +################################ +# nD evaluation implementation # +################################ + +# Overload _return_cache and _setsize for in place D-dimensional de Casteljau algorithm +function _return_cache( + b::BernsteinBasisOnSimplex{D}, x,::Type{G},::Val{N_deriv}) where {D,G,N_deriv} + + @assert D == length(eltype(x)) "Incorrect number of point components" + T = eltype(G) + K = get_order(b) + np = length(x) + ndof = length(b) + ndof_scalar = binomial(K+D,D) + + r = CachedArray(zeros(G,(np,ndof))) + s = MArray{Tuple{Vararg{D,N_deriv}},T}(undef) + c = CachedVector(zeros(T,ndof_scalar)) + # The cache c here holds all scalar nD-Bernstein polynomials, no other caches needed for derivatives + t = ntuple( _ -> nothing, Val(N_deriv)) + (r, s, c, t...) +end + +function _setsize!(b::BernsteinBasisOnSimplex{D}, np, r, t...) where D + K = get_order(b) + ndof = length(b) + ndof_scalar = binomial(K+D,D) + setsize!(r,(np,ndof)) + setsize!(t[1],(ndof_scalar,)) +end + +function _evaluate_nd!( + b::BernsteinBasisOnSimplex{D,V}, x, + r::AbstractMatrix{V}, i, + c::AbstractVector{T}, VK::Val) where {D,V,T} + + λ = _cart_to_bary(x, b.cart_to_bary_matrix) + c[1] = one(T) + _downwards_de_Casteljau_nD!(c,λ,VK,Val(D)) + + k = 1 + for s in c + k = _cartprod_set_value!(r,i,s,k) + end +end + +function _gradient_nd!( + b::BernsteinBasisOnSimplex{D,V}, x, + r::AbstractMatrix{G}, i, + c::AbstractVector{T}, + g::Nothing, + s::MVector{D,T}, ::Val{K}) where {D,V,G,T,K} + + x_to_λ = b.cart_to_bary_matrix + λ = _cart_to_bary(x, x_to_λ) + + c[1] = one(T) + _downwards_de_Casteljau_nD!(c,λ,Val(K-1),Val(D)) + + _grad_Bα_from_Bαm!(r,i,c,s,Val(K),Val(D),V,x_to_λ) +end + +function _hessian_nd!( + b::BernsteinBasisOnSimplex{D,V}, x, + r::AbstractMatrix{G}, i, + c::AbstractVector{T}, + g::Nothing, + h::Nothing, + s::MMatrix{D,D,T}, ::Val{K}) where {D,V,G,T,K} + + x_to_λ = b.cart_to_bary_matrix + λ = _cart_to_bary(x, x_to_λ) + + c[1] = one(T) + _downwards_de_Casteljau_nD!(c,λ,Val(K-2),Val(D)) + + _hess_Bα_from_Bαmm!(r,i,c,s,Val(K),Val(D),V,x_to_λ) +end + +# @generated functions as otherwise the time and allocation for +# computing the indices are the bottlneck... +""" + _downwards_de_Casteljau_nD!(c, λ,::Val{K},::Val{D},::Val{K0}=Val(1)) + +Iteratively applies de Casteljau algorithm in reverse in place using `λ`s as +coefficients. + +If `K0 = 1`, `λ` are the barycentric coordinates of some point `x` and `c[1] = 1`, +this computes all order `K` basis Bernstein polynomials at `x`: + +``c[α_id] = B_α(x)  ∀α in bernstein_terms(K,D)`` + +where `α_id` = [`bernstein_term_id`](@ref)(α). +""" +@generated function _downwards_de_Casteljau_nD!(c, λ,::Val{K},::Val{D},::Val{K0}=Val(1)) where {K,D,K0} + z = zero(eltype(c)) + ex_v = Vector{Expr}() + sub_ids = MVector{D+1,Tuple{Int,Int}}(undef) + + for Ki in K0:K + # Iterations are in reverse lexicographic order (left to right), because α-ei is + # always stored on the left of α (as α-ei < α in lexicographic order), so the + # erased B_β replaced by B_α won't be used to compute the remainings B_γ for |γ|=`K` + # with γ>α in lexicographic order. + terms = bernstein_terms(Ki,D) + for (id,α) in Iterators.reverse(enumerate(terms)) # For all |α| = Ki + # s = 0. + push!(ex_v, :(s = $z)) + + # For all |β| = |α|-1; β ≥ 0 + nb_sα = _sub_multi_indices!(sub_ids, α, Val(D+1)) + for (id_β, d) in take(sub_ids, nb_sα) + # s += λ_d * B_β + push!(ex_v, :(@inbounds s += λ[$d]*c[$id_β])) + end + + # c[id] = B_α + push!(ex_v, :(@inbounds c[$id] = s)) + end + end + return Expr(:block, ex_v...) +end + +""" + _de_Casteljau_nD!(c, λ,::Val{K},::Val{D},::Val{Kf}=Val(0)) + +Iteratively applies de Casteljau algorithm in place using `λ`s as +coefficients. + +If `Kf = 0`, `λ` are the barycentric coordinates of some point `x` and `c` contains +the Bernstein coefficients ``c_α`` of a polynomial ``p`` (that is ``p(x) = ∑_α c_α B_α(x)`` for +``α`` in [`bernstein_terms`](@ref)(`K`,`D`) ), this computes + +``c[1] = p(x)`` + +where the ``c_α`` must be initially stored in `c`[`α_id`], where +`α_id` = [`bernstein_term_id`](@ref)(α). +""" +@generated function _de_Casteljau_nD!(c, λ,::Val{K},::Val{D},::Val{Kf}=Val(0)) where {K,D,Kf} + z = zero(eltype(c)) + ex_v = Vector{Expr}() + sup_ids = MVector{D+1,Tuple{Int,Int}}(undef) + + for Ki in (K-1):-1:Kf + # Iterations are in lexicographic order (right to left), because α+ei is + # always stored on the right of α (as α+ei > α in lexicographic order), so the + # erased B_β replaced by B_α won't be used to compute the remainings B_γ for |γ|=`K` + # with γ<α in lexicographic order. + terms = bernstein_terms(Ki,D) + for (id,α) in enumerate(terms) # For all |α| = Ki + # s = 0. + push!(ex_v, :(s = $z)) + + # For all |β| = |α|+1 + _sup_multi_indices!(sup_ids, α, Val(D+1)) + for (id_β, d) in sup_ids + # s += λ_d * B_β + push!(ex_v, :(@inbounds s += λ[$d]*c[$id_β])) + end + + # c[id] = B_α (= s) + push!(ex_v, :(@inbounds c[$id] = s)) + end + end + return Expr(:block, ex_v...) +end + + +# ∂q(B_α) = K ∑_{1 ≤ i ≤ N} ∂q(λi) B_β +# for 1 ≤ q ≤ D and β = α-ei +@generated function _grad_Bα_from_Bαm!( + r,i,c,s,::Val{K},::Val{D},::Type{V},x_to_λ=nothing) where {K,D,V} + + ex_v = Vector{Expr}() + ncomp = num_indep_components(V) + z = zero(eltype(c)) + N = D+1 + δ(i,j) = Int(i==j) + sub_ids = MVector{D+1,Tuple{Int,Int}}(undef) + + for (id_α,α) in enumerate(bernstein_terms(K,D)) + push!(ex_v, :(@inbounds s .= $z)) # s = 0 + + nb_sα = _sub_multi_indices!(sub_ids, α, Val(D+1)) + for (id_β, i) in take(sub_ids, nb_sα) # β = α - ei + push!(ex_v, :(@inbounds B_β = c[$id_β])) + # s[q] = Σ_β ∂q(λi) B_β + for q in 1:D + if x_to_λ == Nothing + # ∂q(λi) = δiq - δiN + Cqi = δ(i,q) - δ(i,N) + iszero(Cqi) || push!(ex_v, :(@inbounds s[$q] += $Cqi*B_β)) + else + # ∂q(λi) = ei (x_to_λ*(e1 - e{q+1}) - x_to_λ*(e1)) = ei * x_to_λ * e{q+1} + # ∂q(λi) = x_to_λ[i,q+1] + push!(ex_v, :(@inbounds s[$q] += x_to_λ[$i,$(q+1)]*B_β)) + end + end + end + push!(ex_v, :(@inbounds s .*= $K)) # s = Ks. + + k = ncomp*(id_α-1) + 1 + push!(ex_v, :(_cartprod_set_derivative!(r,i,s,$k,V))) end + + return Expr(:block, ex_v...) end + +# ∂t(∂q(B_α)) = K(K-1) ∑_{1 ≤ i,j ≤ N} ∂t(λj) ∂q(λi) B_β +# for 1 ≤ t,q ≤ D and β = α - ei - ej +@generated function _hess_Bα_from_Bαmm!( + r,i,c,s,::Val{K},::Val{D},::Type{V},x_to_λ=nothing) where {K,D,V} + + ex_v = Vector{Expr}() + ncomp = num_indep_components(V) + z = zero(eltype(c)) + N = D+1 + δ(i,j) = Int(i==j) + C(q,t,i,j) = (δ(i,q)-δ(i,N))*(δ(j,t)-δ(j,N)) + N_max_ssα = binomial(D+2,2) + sub_sub_α_ids = MVector{N_max_ssα, NTuple{3,Int}}(undef) + + for (id_α,α) in enumerate(bernstein_terms(K,D)) + push!(ex_v, :(@inbounds s .= $z)) # s = 0 + + nb_subsub = _sub_sub_multi_indices!(sub_sub_α_ids, α, Val(D+1)) + for (id_β, i, j) in take(sub_sub_α_ids, nb_subsub) + + push!(ex_v, :(@inbounds B_β = c[$id_β])) + for t in 1:D + for q in 1:D + if x_to_λ == Nothing + # s[t,q] = Σ_β B_β (δ_iq - δ_iN)*(δ_jt - δ_jN) + Cβ = C(q,t,i,j) + if i≠j Cβ += C(q,t,j,i) end + iszero(Cβ) || push!(ex_v, :(@inbounds s[$t,$q] += $Cβ*B_β)) + else + push!(ex_v, :(@inbounds C = x_to_λ[$i,$(q+1)]*x_to_λ[$j,$(t+1)])) + if i≠j push!(ex_v, :(@inbounds C += x_to_λ[$j,$(q+1)]*x_to_λ[$i,$(t+1)])) end + push!(ex_v, :(@inbounds s[$t,$q] += C*B_β)) + end + end + end + end + push!(ex_v, :(@inbounds s .*= $(K*(K-1))) ) # s = K(K-1)s + + k = ncomp*(id_α-1) + 1 + push!(ex_v, :(_cartprod_set_derivative!(r,i,s,$k,V))) + end + + return Expr(:block, ex_v...) +end + + +######################## +# de Casteljau helpers # +######################## + +""" + bernstein_term_id(α) + +For a given Bernstein multi-index `α` (vector or tuple), return the associated +linear index of `α` ordered in decreasing lexicographic order, that is the `i` +such that + + (i,α) ∈ enumerate(bernstein_terms(K,D)) + +where K = sum(`α`), see also [`bernstein_terms`](@ref). +""" +function bernstein_term_id(α) + D = length(α)-1 + @check D ≥ 0 + @inbounds i = sum(_L_slices_size(L, D, _L_slice_2(L,α,D)) for L in 1:D; init=0) + 1 + return i +end + + +# For a given Bernstein term `α`, return the index (starting from 0) of the +# (D-`L`)-slice to which `α` belongs within the (D-`L`-1)-slice of +# the D-multiexponent simplex (D = `N`-1). +# +# In a D-multiexponent simplex of elements `α`, ordered in a vector in +# decreasing lexicographic order, the (D-`L`)-slices are the consecutive `α`s +# having iddentical first `L` indices `α`. +""" + _L_slice(L, α, D) + +where `L` ∈ 1:N, returns `sum( α[i] for i in (L+1):(D+1) )`. +""" +Base.@propagate_inbounds _L_slice_2(L, α, D) = sum( α[i] for i in (L+1):(D+1) ) + +# Return the length of the `l`-1 first (`D`-`L`)-slices in the +# `D`-multiexponent simplex (ordered in decreasing lexicographic order). +# Those numbers are the "(`D`-`L`)-simplex numbers". +""" + _L_slices_size(L,D,l) = binomial(D-L+l, D-L+1) +""" +_L_slices_size(L,D,l) = binomial(D-L+l, D-L+1) + +""" + _sub_multi_indices!(sub_ids, α) + +Given a positive multi-index `α`, sets in place in `sub_ids` the couples +(`id`, `d`) with `d` in 1:`N` for which the multi-index `αd⁻` = `α`-e`d` is +positive (that is `α`[`d`]>0), and `id` is the linear index of `αd⁻` +(see [`bernstein_term_id`](@ref)). + +The function returns the number of sub indices set. +""" +function _sub_multi_indices!(sub_ids, α, ::Val{N}) where N + @check length(sub_ids) >= N + nb_sα = 0 + for i in 1:N + α⁻ = ntuple(k -> α[k]-Int(k==i), Val(N)) + if all(α⁻ .≥ 0) + nb_sα += 1 + id⁻ = bernstein_term_id(α⁻) + sub_ids[nb_sα] = (id⁻, i) + end + end + return nb_sα +end + +""" + _sub_sub_multi_indices!(sub_ids, α, ::Val{N}) + +Like [`_sub_multi_indices`](@ref), but sets the triples (`id`, `t`, `q`) in `sub_ids`, +with `t,q` in 1:`N` for which the multi-index `αd⁻⁻` = `α`-e`t`-e`q` is positive, +and returns the number of triples set. +""" +function _sub_sub_multi_indices!(sub_ids, α, ::Val{N}) where N + @check length(sub_ids) >= binomial(N+1,2) + nb_ssα = 0 + for i in 1:N + for j in i:N + α⁻⁻ = ntuple(k -> α[k]-Int(k==i)-Int(k==j), Val(N)) + if all(α⁻⁻ .≥ 0) + nb_ssα += 1 + id⁻⁻ = bernstein_term_id(α⁻⁻) + sub_ids[nb_ssα] = (id⁻⁻, i, j) + end + end + end + return nb_ssα +end + +""" + _sup_multi_indices!(sup_ids, α, ::Val{N}) + +Like [`_sub_multi_indices!`](@ref), but sets the indices for the `N` multi-indices +`αd⁺` = `α`+e`d` for 1≤d≤`N`, and returns `N` +""" +function _sup_multi_indices!(sup_ids, α, ::Val{N}) where N + @check length(sup_ids) >= N + for i in 1:N + α⁺ = ntuple(k -> α[k]+Int(k==i), Val(N)) + id⁺ = bernstein_term_id(α⁺) + sup_ids[i] = (id⁺, i) + end + return N +end + + +####################################################### +#### Bernstein bases on simplices Naive implementation# +####################################################### +# +# """ +# BernsteinBasisOnSimplex{D,V} <: PolynomialBasis{D,V,Bernstein} +# +# This basis uses barycentric coordinates defined by the vertices of the +# reference `D`-simplex. +# """ +# struct BernsteinBasisOnSimplex{D,V} <: PolynomialBasis{D,V,Bernstein} +# function BernsteinBasisOnSimplex{D}(::Type{V},order::Int) where {D,V} +# K = Int(order) +# new{D,V}() +# end +# end +# +# function BernsteinBasisOnSimplex(::Val{D},::Type{V},order::Int) where {D,V} +# BernsteinBasisOnSimplex{D}(V,order) +# end +# +# Base.size(::BernsteinBasisOnSimplex{D,V}) where {D,V} = (num_indep_components(V)*binomial(D+K,D),) +# get_exponents(::BernsteinBasisOnSimplex{D,V}) where {D,V} = bernstein_terms(K,D) +# +# ################################ +# # nD evaluation implementation # +# ################################ +# +# # Overload _return_cache and _setsize to add +1 coordinate cache in t +# function _return_cache( +# f::BernsteinBasisOnSimplex{D}, x,::Type{G},::Val{N_deriv}) where {D,G,N_deriv} +# +# @assert D == length(eltype(x)) "Incorrect number of point components" +# T = eltype(G) +# np = length(x) +# ndof = length(f) +# ndof_1d = get_order(f) + 1 +# r = CachedArray(zeros(G,(np,ndof))) +# s = MArray{Tuple{Vararg{D,N_deriv}},T}(undef) +# bernstein_D = D+1 # There are D+1 barycentric coordinates +# t = ntuple( _ -> CachedArray(zeros(T,(bernstein_D,ndof_1d))), Val(N_deriv+1)) +# (r, s, t...) +# end +# function _setsize!(f::BernsteinBasisOnSimplex{D}, np, r, t...) where D +# ndof = length(f) +# ndof_1d = get_order(f) + 1 +# setsize!(r,(np,ndof)) +# bernstein_D = D+1 # There are D+1 barycentric coordinates +# for c in t +# setsize!(c,(bernstein_D,ndof_1d)) +# end +# end +# +# +# function _evaluate_nd!( +# b::BernsteinBasisOnSimplex{D,V}, x, +# r::AbstractMatrix{V}, i, +# c::AbstractMatrix{T}) where {D,V,T} +# +# K = get_order(b) +# terms = _get_terms(b) +# coefs = multinoms(Val(K),Val(D)) +# +# λ = _cart_to_bary(x,nothing) +# +# for d in 1:(D+1) +# _evaluate_1d!(Monomial,Val(K),c,λ,d) # compute powers 0:K of all bary. coords. +# end +# +# k = 1 +# for (ci,m) in zip(terms,coefs) +# +# for d in 1:(D+1) +# @inbounds m *= c[d,ci[d]] +# end +# +# k = _cartprod_set_value!(r,i,m,k) +# end +# end +# +# function _gradient_nd!( +# b::BernsteinBasisOnSimplex{D,V}, x, +# r::AbstractMatrix{G}, i, +# c::AbstractMatrix{T}, +# g::AbstractMatrix{T}, +# s::MVector{D,T}) where {D,V,G,T} +# +# K = get_order(b) +# N = D+1 +# terms = _get_terms(b) +# coefs = multinoms(Val(K),Val(D)) +# +# λ = _cart_to_bary(x,nothing) +# +# for d in 1:N +# _derivatives_1d!(Monomial,Val(K),(c,g),λ,d) +# end +# +# k = 1 +# @inbounds for (ci,m) in zip(terms,coefs) +# +# for i in eachindex(s) +# s[i] = m +# end +# +# for q in 1:D +# for d in 1:D +# if d != q +# s[q] *= c[d,ci[d]] +# else +# s[q] *= g[q,ci[q]]*c[N,ci[N]] - g[N,ci[N]]*c[q,ci[q]] +# end +# end +# end +# +# k = _cartprod_set_derivative!(r,i,s,k,V) +# end +# end +# +# function _hessian_nd!( +# b::BernsteinBasisOnSimplex{D,V}, x, +# r::AbstractMatrix{G}, i, +# c::AbstractMatrix{T}, +# g::AbstractMatrix{T}, +# h::AbstractMatrix{T}, +# s::MMatrix{D,D,T}) where {D,V,G,T} +# +# N = D+1 +# terms = _get_terms(b) +# coefs = multinoms(Val(K),Val(D)) +# +# λ = _cart_to_bary(x,nothing) +# +# for d in 1:N +# _derivatives_1d!(Monomial,Val(K),(c,g,h),λ,d) +# end +# +# k = 1 +# @inbounds for (ci,m) in zip(terms,coefs) +# +# for i in eachindex(s) +# s[i] = m +# end +# +# for t in 1:D +# for q in 1:D +# for d in 1:D +# if d != q && d != t +# # if q == t, D-1 factors +# # else, D-2 factors +# s[t,q] *= c[d,ci[d]] +# elseif q == t # == d +# # +2 factors -> D+1 +# s[t,q] *= (h[d,ci[d]]*c[N,ci[N]] -2g[d,ci[d]]*g[N,ci[N]] + c[d,ci[d]]*h[N,ci[N]]) +# elseif d == q # q ≠ t, we multiply once with the factors with q and t derivative terms +# # +3 factors -> D+1 +# s[t,q] *=( g[t,ci[t]]*g[q,ci[q]]*c[N,ci[N]] +# - g[t,ci[t]]*c[q,ci[q]]*g[N,ci[N]] +# - c[t,ci[t]]*g[q,ci[q]]*g[N,ci[N]] +# + c[t,ci[t]]*c[q,ci[q]]*h[N,ci[N]]) +# end +# end +# end +# end +# +# k = _cartprod_set_derivative!(r,i,s,k,V) +# end +# end +# +# """ +# multinoms(::Val{K}, ::Val{D}) +# +# Returns the tuple of multinomial coefficients for each term in +# [`bernstein_terms`](@ref)(`K`,`D`). For e.g. a term `t`, the +# multinomial can be computed by `factorial(sum(t)) ÷ prod(factorial.(t)` +# """ +# @generated function multinoms(::Val{K},::Val{D}) where {K,D} +# terms = bernstein_terms(K,D) +# multinomials = tuple( (multinomial(α...) for α in terms)... ) +# Meta.parse("return $multinomials") +# end +# diff --git a/src/Polynomials/UniformPolyBases.jl b/src/Polynomials/CartProdPolyBases.jl similarity index 71% rename from src/Polynomials/UniformPolyBases.jl rename to src/Polynomials/CartProdPolyBases.jl index 8068a9c48..81555d0fd 100644 --- a/src/Polynomials/UniformPolyBases.jl +++ b/src/Polynomials/CartProdPolyBases.jl @@ -3,15 +3,17 @@ ################################# """ - 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 @@ -19,18 +21,18 @@ The scalar polynomial basis spanning 𝕊 is defined as 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,7 +236,7 @@ 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 @@ -245,7 +244,7 @@ function _uniform_set_derivative!( 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 diff --git a/src/Polynomials/ChebyshevBases.jl b/src/Polynomials/ChebyshevBases.jl index 913b58e21..f5e02c3f9 100644 --- a/src/Polynomials/ChebyshevBases.jl +++ b/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" diff --git a/src/Polynomials/CompWiseTensorPolyBases.jl b/src/Polynomials/CompWiseTensorPolyBases.jl index f0206e125..9895160c5 100644 --- a/src/Polynomials/CompWiseTensorPolyBases.jl +++ b/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,7 +20,8 @@ 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}( @@ -28,29 +29,30 @@ struct CompWiseTensorPolyBasis{D,V,K,PT,L} <: PolynomialBasis{D,V,K,PT} 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,7 +172,7 @@ 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 @@ -188,26 +180,19 @@ end 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 diff --git a/src/Polynomials/LegendreBases.jl b/src/Polynomials/LegendreBases.jl index 9b5109d47..197295baf 100644 --- a/src/Polynomials/LegendreBases.jl +++ b/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 diff --git a/src/Polynomials/ModalC0Bases.jl b/src/Polynomials/ModalC0Bases.jl index 42cc40a0c..cd29f1d4a 100644 --- a/src/Polynomials/ModalC0Bases.jl +++ b/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} diff --git a/src/Polynomials/MonomialBases.jl b/src/Polynomials/MonomialBases.jl index f37a7ab55..d451e75ce 100644 --- a/src/Polynomials/MonomialBases.jl +++ b/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,7 +21,7 @@ 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) @@ -29,7 +29,7 @@ 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 + diff --git a/src/Polynomials/NedelecPolyBases.jl b/src/Polynomials/NedelecPolyBases.jl index a5d49af9d..6f2847ca1 100644 --- a/src/Polynomials/NedelecPolyBases.jl +++ b/src/Polynomials/NedelecPolyBases.jl @@ -1,14 +1,14 @@ """ - 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" @@ -16,23 +16,26 @@ struct NedelecPolyBasisOnSimplex{D,V,K,PT} <: PolynomialBasis{D,V,K,PT} @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 diff --git a/src/Polynomials/PolynomialInterfaces.jl b/src/Polynomials/PolynomialInterfaces.jl index bdabfde1a..74615363d 100644 --- a/src/Polynomials/PolynomialInterfaces.jl +++ b/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,28 +44,29 @@ 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 ########### @@ -73,9 +74,9 @@ Return the maximum polynomial order in a dimension, is should be `0` in 0D. ########### _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 diff --git a/src/Polynomials/Polynomials.jl b/src/Polynomials/Polynomials.jl index b4fb645f2..00d2eadd5 100644 --- a/src/Polynomials/Polynomials.jl +++ b/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,6 +80,7 @@ module Polynomials using DocStringExtensions using LinearAlgebra: mul! +using LinearAlgebra: I using StaticArrays using Gridap.Helpers using Gridap.Arrays @@ -87,6 +88,8 @@ 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,7 +106,7 @@ export Bernstein export PolynomialBasis export get_order -export UniformPolyBasis +export CartProdPolyBasis export get_exponents export get_orders export MonomialBasis @@ -111,6 +114,10 @@ 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") diff --git a/src/Polynomials/RaviartThomasPolyBases.jl b/src/Polynomials/RaviartThomasPolyBases.jl index bc810f6ca..9178002e9 100644 --- a/src/Polynomials/RaviartThomasPolyBases.jl +++ b/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,11 +52,12 @@ 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 ################################# @@ -63,16 +65,12 @@ Base.size(a::RaviartThomasPolyBasis{D}) where {D} = (D*length(a.pterms) + length ################################# 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)) @@ -80,7 +78,7 @@ function _evaluate_nd!( @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,18 +112,14 @@ 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)) @@ -133,7 +127,7 @@ function _gradient_nd!( @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 diff --git a/src/ReferenceFEs/Dofs.jl b/src/ReferenceFEs/Dofs.jl index 546afc1aa..913311b90 100644 --- a/src/ReferenceFEs/Dofs.jl +++ b/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} diff --git a/src/TensorValues/MultiValueTypes.jl b/src/TensorValues/MultiValueTypes.jl index cc73d31aa..cb36a5d49 100644 --- a/src/TensorValues/MultiValueTypes.jl +++ b/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) diff --git a/test/PolynomialsTests/BernsteinBasesTests.jl b/test/PolynomialsTests/BernsteinBasesTests.jl index 60abe1ccd..75d7dbe3a 100644 --- a/test/PolynomialsTests/BernsteinBasesTests.jl +++ b/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,9 +137,9 @@ 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) @@ -144,4 +147,199 @@ 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 diff --git a/test/PolynomialsTests/ChebyshevBasesTests.jl b/test/PolynomialsTests/ChebyshevBasesTests.jl index b32ddfd90..d7a93062b 100644 --- a/test/PolynomialsTests/ChebyshevBasesTests.jl +++ b/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 diff --git a/test/PolynomialsTests/ModalC0BasesTests.jl b/test/PolynomialsTests/ModalC0BasesTests.jl index 3ba41a78e..51dc58b8c 100644 --- a/test/PolynomialsTests/ModalC0BasesTests.jl +++ b/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,7 +83,7 @@ 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 @@ -91,7 +91,7 @@ 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) diff --git a/test/PolynomialsTests/MonomialBasesTests.jl b/test/PolynomialsTests/MonomialBasesTests.jl index d3a734982..9604b1949 100644 --- a/test/PolynomialsTests/MonomialBasesTests.jl +++ b/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 diff --git a/test/PolynomialsTests/PCurlGradBasesTests.jl b/test/PolynomialsTests/PCurlGradBasesTests.jl index 07d171796..5e31085e4 100644 --- a/test/PolynomialsTests/PCurlGradBasesTests.jl +++ b/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) diff --git a/test/PolynomialsTests/PGradBasesTests.jl b/test/PolynomialsTests/PGradBasesTests.jl index 01138775d..b7e8569cf 100644 --- a/test/PolynomialsTests/PGradBasesTests.jl +++ b/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) diff --git a/test/PolynomialsTests/PolynomialInterfacesTests.jl b/test/PolynomialsTests/PolynomialInterfacesTests.jl index 7d4660761..3e0e02079 100644 --- a/test/PolynomialsTests/PolynomialInterfacesTests.jl +++ b/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 diff --git a/test/PolynomialsTests/QCurlGradBasesTests.jl b/test/PolynomialsTests/QCurlGradBasesTests.jl index 30a7d3207..8acf3feab 100644 --- a/test/PolynomialsTests/QCurlGradBasesTests.jl +++ b/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) diff --git a/test/PolynomialsTests/QGradBasesTests.jl b/test/PolynomialsTests/QGradBasesTests.jl index 1a7d0efd1..130c6e64b 100644 --- a/test/PolynomialsTests/QGradBasesTests.jl +++ b/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) diff --git a/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl b/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl index 4dc8dee38..1adf51774 100644 --- a/test/ReferenceFEsTests/CLagrangianRefFEsTests.jl +++ b/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)