From 7ab734164927bddd9221edcec2cb7a5a6334754f Mon Sep 17 00:00:00 2001 From: Antoine Marteau Date: Fri, 10 Jan 2025 18:15:12 +1100 Subject: [PATCH] make ModalC0 compatible with all basis parametrised by Polynomial --- docs/src/Polynomials.md | 26 ++++---- src/Polynomials/ModalC0Bases.jl | 69 ++++++++++++++++++---- test/PolynomialsTests/ModalC0BasesTests.jl | 45 +++++++++++++- 3 files changed, 117 insertions(+), 23 deletions(-) diff --git a/docs/src/Polynomials.md b/docs/src/Polynomials.md index 58143c701..c1f116dea 100644 --- a/docs/src/Polynomials.md +++ b/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 diff --git a/src/Polynomials/ModalC0Bases.jl b/src/Polynomials/ModalC0Bases.jl index 77d88ac8a..42cc40a0c 100644 --- a/src/Polynomials/ModalC0Bases.jl +++ b/src/Polynomials/ModalC0Bases.jl @@ -27,6 +27,8 @@ 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) @@ -34,6 +36,25 @@ struct ModalC0Basis{D,V,T,K} <: PolynomialBasis{D,V,K,ModalC0} 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 diff --git a/test/PolynomialsTests/ModalC0BasesTests.jl b/test/PolynomialsTests/ModalC0BasesTests.jl index de293238e..04c4f839b 100644 --- a/test/PolynomialsTests/ModalC0BasesTests.jl +++ b/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