make ModalC0 compatible with all basis parametrised by Polynomial

docs/src/Polynomials.md

@@ -96,29 +96,35 @@ polynomials ``\phi_K`` for which ``\phi_K(0) = \delta_{0K}`` and ``\phi_K(1) =
 \delta_{1K}``. This module implements the generalised version introduced in Eq.
 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
+bounding box `(a=Point{D}(0...), b=Point{D}(1...))` is assumed, which
+coincides with the usual definition of the ModalC0 bases.
+
 
 ### Polynomial spaces in FEM
 
 #### 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``
 ```math
-\mathbb{P}_K(x) = \text{Span}\big\{\quad x\rightarrow x^i \quad\big|\quad 0\leq i\leq K \quad\big\}
+\mathbb{P}_K(x) = \text{Span}\big\{\quad x\rightarrow x^i \quad\big|\quad 0\leq i\leq K \quad\big\}.
 ```
 
+Then, ``\mathbb{Q}^D`` and ``\mathbb{P}^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\}
+    \alpha_1, \alpha_2, \dots, \alpha_D \leq 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
-    K \quad\big\}
+    K \quad\big\}.
 ```
 
-```math
-\mathbb{P}_K = \mathbb{P}^1_K = \mathbb{Q}^1_K
-```
+To note, there is ``\mathbb{P}_K = \mathbb{P}^1_K = \mathbb{Q}^1_K``.
 
 #### Serendipity space Sr
 
@@ -263,10 +269,6 @@ BernsteinBasis(args...)
 
 ```@docs
 BernsteinBasis
-ModalC0Basis
-```
-
-```@docs
 PGradBasis
 QGradBasis
 PCurlGradBasis
@@ -285,6 +287,8 @@ get_exponents
 CompWiseTensorPolyBasis
 NedelecPolyBasisOnSimplex
 RaviartThomasPolyBasis
+ModalC0Basis
+ModalC0Basis()
 ```
 
 ### Deprecated APIs

src/Polynomials/ModalC0Bases.jl

@@ -27,13 +27,34 @@ struct ModalC0Basis{D,V,T,K} <: PolynomialBasis{D,V,K,ModalC0}
     a::Vector{Point{D,T}},
     b::Vector{Point{D,T}}) where {D,V,T}
 
+    _msg =  "The number of bounding box points in a and b should match the number of terms"
+    @check length(terms) == length(a) == length(b) _msg
     @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)
   end
 end
 
+"""
+    ModalC0Basis{D}(::Type{V},order::Int;                        [, filter][, sort!:])
+    ModalC0Basis{D}(::Type{V},order::Int,   a::Point ,b::Point ; [, filter][, sort!:])
+    ModalC0Basis{D}(::Type{V},orders::Tuple;                     [, filter][, sort!:])
+    ModalC0Basis{D}(::Type{V},orders::Tuple,a::Point ,b::Point ; [, filter][, sort!:])
+    ModalC0Basis{D}(::Type{V},orders::Tuple,a::Vector,b::Vector; [, filter][, sort!:])
+
+where `filter` is a `Function` defaulting to `_q_filter`, and `sort!` is a
+`Function` defaulting to `_sort_by_nfaces!`.
+
+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.
+"""
+function ModalC0Basis() end
+
 function ModalC0Basis{D}(
   ::Type{V},
   orders::NTuple{D,Int},
@@ -94,6 +115,7 @@ function ModalC0Basis{D}(
   ModalC0Basis{D}(V,orders; filter=filter, sort! =sort!)
 end
 
+
 # API
 
 @inline Base.size(a::ModalC0Basis{D,V}) where {D,V} = (length(a.terms)*num_indep_components(V),)
@@ -352,12 +374,12 @@ end
 function _gradient_1d_mc0!(g::AbstractMatrix{T},x,a,b,order,d) where T
   @assert order > 0
   n = order + 1
-  z = one(T)
-  @inbounds g[d,1] = -z
-  @inbounds g[d,2] = z
+  o = one(T)
+  @inbounds g[d,1] = -o
+  @inbounds g[d,2] = o
   if n > 2
     ξ = ( 2*x[d] - ( a[d] + b[d] ) ) / ( b[d] - a[d] )
-    v1 = z - x[d]
+    v1 = o - x[d]
     v2 = x[d]
     for i in 3:n
       j, dj = jacobi_and_derivative(ξ,i-3,1,1)
@@ -369,16 +391,16 @@ end
 function _hessian_1d_mc0!(h::AbstractMatrix{T},x,a,b,order,d) where T
   @assert order > 0
   n = order + 1
-  y = zero(T)
-  z = one(T)
-  @inbounds h[d,1] = y
-  @inbounds h[d,2] = y
+  z = zero(T)
+  o = one(T)
+  @inbounds h[d,1] = z
+  @inbounds h[d,2] = z
   if n > 2
     ξ = ( 2*x[d] - ( a[d] + b[d] ) ) / ( b[d] - a[d] )
-    v1 = z - x[d]
+    v1 = o - x[d]
     v2 = x[d]
-    dv1 = -z
-    dv2 = z
+    dv1 = -o
+    dv2 = o
     for i in 3:n
       j, dj = jacobi_and_derivative(ξ,i-3,1,1)
       _, d2j = jacobi_and_derivative(ξ,i-4,2,2)
@@ -387,3 +409,28 @@ function _hessian_1d_mc0!(h::AbstractMatrix{T},x,a,b,order,d) where T
   end
 end
 
+
+#######################################
+# Generic 1D internal polynomial APIs #
+#######################################
+
+# For possible use with UniformPolyBasis 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}
+  a = zero(x)
+  b = zero(x) .+ one(T)
+  @inline _evaluate_1d_mc0!(c,x,a,b,K,d)
+end
+
+function _gradient_1d!(::Type{ModalC0},::Val{K},g::AbstractMatrix{T},x,d) where {K,T<:Number}
+  a = zero(x)
+  b = zero(x) .+ one(T)
+  @inline _gradient_1d_mc0!(g,x,a,b,K,d)
+end
+
+function _hessian_1d!(::Type{ModalC0},::Val{K},h::AbstractMatrix{T},x,d) where {K,T<:Number}
+  a = zero(x)
+  b = zero(x) .+ one(T)
+  @inline _hessian_1d_mc0!(h,x,a,b,K,d)
+end

test/PolynomialsTests/ModalC0BasesTests.jl

@@ -39,6 +39,24 @@ 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
+
+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)
+
+∇b1  = Broadcasting(∇)(b1)
+∇b1u = Broadcasting(∇)(b1u)
+∇∇b1 = Broadcasting(∇)(∇b1)
+∇∇b1u= Broadcasting(∇)(∇b1u)
+
+@test evaluate(b1,  [x1,x2,x3,]) ≈ evaluate(b1u,  [x1,x2,x3,])
+@test evaluate(∇b1, [x1,x2,x3,]) ≈ evaluate(∇b1u, [x1,x2,x3,])
+@test evaluate(∇∇b1,[x1,x2,x3,]) ≈ evaluate(∇∇b1u,[x1,x2,x3,])
+
+
 x1 = Point(0.0,0.0)
 x2 = Point(0.5,0.5)
 x3 = Point(1.0,1.0)
@@ -51,7 +69,6 @@ b = [ Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0),
       Point(-1.0,1.0),Point(-1.0,1.0),Point(-1.0,1.25),Point(-1.0,1.25),
       Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0),Point(1.0,1.0) ]
 b2 = ModalC0Basis{2}(V,order,a,b)
-#b2 = ModalC0Basis(Val(2),V,order,a,b)
 ∇b2 = Broadcasting(∇)(b2)
 ∇∇b2 = Broadcasting(∇)(∇b2)
 
@@ -65,6 +82,32 @@ H = gradient_type(G,x1)
 @test evaluate(∇∇b2,[x1,x2,x3,])[:,10] ≈ H[ (0.0, 0.0, 0.0, -4.47213595499958);
                                             (0.0, 0.5590169943749475, 0.5590169943749475, 1.118033988749895);
                                             (0.0, -2.23606797749979, -2.23606797749979, 0.0) ]
+
+# Validate generic 2D implem using UniformPolyBasis
+
+order = 3
+len_b2 = (order+1)^2
+a = fill(Point(0.,0.), len_b2)
+b = fill(Point(1.,1.), len_b2)
+
+b2 = ModalC0Basis{2}(V,order,a,b)
+b2u= UniformPolyBasis(ModalC0,Val(2),V,order)
+∇b2  = Broadcasting(∇)(b2)
+∇b2u = Broadcasting(∇)(b2u)
+
+b2x   = eachcol(evaluate(b2,    [x1,x2,x3,]))
+b2xu  = eachcol(evaluate(b2u,   [x1,x2,x3,]))
+∇b2x  = eachcol(evaluate(∇b2,   [x1,x2,x3,]))
+∇b2xu = eachcol(evaluate(∇b2u,  [x1,x2,x3,]))
+
+# re order basis polynomials as each basis has different ordering ...
+b2x_perm   = b2x[  sortperm(b2x)[  invperm(sortperm(b2xu))]]
+∇b2x_perm  = ∇b2x[ sortperm(∇b2x)[ invperm(sortperm(∇b2xu))]]
+
+@test b2xu  == b2x_perm
+@test ∇b2xu == ∇b2x_perm
+
+
 # Misc
 
 # Derivatives not implemented for symetric tensor types