Merge pull request #1294 from gridap/bugfix-basis-type-instability

Bugfix on basis type instability

Author: JordiManyer

NEWS.md

@@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 
 - Fixed type instability for tensor operations between `MultiValue` and scalars. Since PR[#1293](https://github.com/gridap/Gridap.jl/pull/1293).
+- Fixed type instability in basis construction when user gives a non-concrete output type. Since PR[#1294](https://github.com/gridap/Gridap.jl/pull/1294).
 
 ## [0.20.5] - 2026-04-28
 

src/Polynomials/BernsteinBases.jl

@@ -247,10 +247,11 @@ struct BernsteinBasisOnSimplex{D,V,M,K} <: PolynomialBasis{D,V,Bernstein}
   function BernsteinBasisOnSimplex{D}(::Type{V},order::Int,vertices=nothing) where {D,V}
     _simplex_vertices_checks(Val(D), vertices)
 
+    VV = make_concretetype(V)
     cart_to_bary_matrix = _compute_cart_to_bary_matrix(vertices, Val(D+1))
     M = typeof(cart_to_bary_matrix) # Nothing or SMatrix
     K = order
-    new{D,V,M,K}(cart_to_bary_matrix)
+    new{D,VV,M,K}(cart_to_bary_matrix)
   end
 end
 

src/Polynomials/CartProdPolyBases.jl

@@ -39,12 +39,13 @@ struct CartProdPolyBasis{D,V,PT} <: PolynomialBasis{D,V,PT}
     orders::NTuple{D,Int},
     terms::Vector{CartesianIndex{D}}) where {D,V,PT<:Polynomial}
 
+    VV = make_concretetype(V)
     @check isconcretetype(PT) "PT needs to be a concrete <:Polynomial type"
 
     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,PT}(K,orders,terms)
+    new{D,VV,PT}(K,orders,terms)
   end
 end
 

src/Polynomials/CompWiseTensorPolyBases.jl

@@ -59,12 +59,13 @@ struct CompWiseTensorPolyBasis{D,V,PT} <: PolynomialBasis{D,V,PT}
     @check size(orders,2) == D msg3
     @check isconcretetype(PT) "PT needs to be a concrete <:Polynomial type"
 
+    VV = make_concretetype(V)
     K = maximum(orders)
     num_poly = mapreduce(length, +, comp_terms)
 
     #TODO check orders in `orders` greater or equal than max index in terms
 
-    new{D,V,PT}(num_poly,K,orders,comp_terms)
+    new{D,VV,PT}(num_poly,K,orders,comp_terms)
   end
 end
 

src/Polynomials/ModalC0Bases.jl

@@ -43,9 +43,10 @@ struct ModalC0Basis{D,V,T} <: PolynomialBasis{D,V,ModalC0}
     _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"
+    VV = make_concretetype(V)
     K = maximum(orders, init=0)
 
-    new{D,V,T}(K,orders,terms,a,b)
+    new{D,VV,T}(K,orders,terms,a,b)
   end
 end
 

src/TensorValues/MultiValueTypes.jl

@@ -97,6 +97,11 @@ function rand(rng::AbstractRNG,::Random.SamplerType{V}) where V<:MultiValue{D,T}
   V(Tuple(vrand))
 end
 
+function make_concretetype(::Type{T}) where T <: Number
+  TT = ifelse(isconcretetype(T),T,typeof(zero(T)))
+  @check isconcretetype(TT) "Type $(T) cannot be made concrete."
+  return TT
+end
 
 ## ATM it is not possible to implement array like axes because lazy_mapping
 ## operations / broadcast rely on axes(::MultiValue) adopting the Number convention to return ().

src/TensorValues/TensorValues.jl

@@ -66,6 +66,7 @@ export skew_symmetric_part
 export num_components
 export num_indep_components
 export change_eltype
+export make_concretetype
 export diagonal_tensor
 export ⊙
 export ⊗

test/PolynomialsTests/BernsteinBasesTests.jl

@@ -345,4 +345,11 @@ Hbx = _Hbx(D,order,x,H,x2λ)
 test_field_array(b,x,bx,≈, grad=Gbx, gradgrad=Hbx)
 test_field_array(b,x1,bx[1,:],≈,grad=Gbx[1,:],gradgrad=Hbx[1,:])
 
+# value_type must be concrete even when user passes a non-concrete tensor type
+for V in (TensorValue{2,2,Float64}, SymTensorValue{2,Float64},
+          SkewSymTensorValue{2,Float64}, SymFourthOrderTensorValue{2,Float64})
+  @test isconcretetype(value_type(BernsteinBasis(Val(2), V, 1)))
+  @test isconcretetype(value_type(BernsteinBasisOnSimplex{2}(V, 1)))
+end
+
 end # module

test/PolynomialsTests/CompWiseTensorPolyBasesTests.jl

@@ -59,5 +59,10 @@ Gbx = hcat( [ getindex.(bi,1) .* VectorValue(1.,0)⊗V(1.,0) for bi in evaluate(
 
 test_field_array(b,x,bx,≈, grad=Gbx)
 
+# value_type must be concrete even when user passes a non-concrete tensor type
+for V in (TensorValue{2,2,Float64}, SymTensorValue{2,Float64}, SymFourthOrderTensorValue{2,Float64})
+  L = num_indep_components(V)
+  @test isconcretetype(value_type(CompWiseTensorPolyBasis{2}(Monomial, V, ones(Int, L, 2))))
+end
 
 end # module

test/PolynomialsTests/ModalC0BasesTests.jl

@@ -141,5 +141,10 @@ G = gradient_type(V,x1)
 r = zeros(G, (1,1))
 @test_throws ErrorException Polynomials._set_derivative_mc0!(r,1,s,0,0,V)
 
+# value_type must be concrete even when user passes a non-concrete tensor type
+for V in (TensorValue{2,2,Float64}, SymTensorValue{2,Float64},
+          SkewSymTensorValue{2,Float64}, SymFourthOrderTensorValue{2,Float64})
+  @test isconcretetype(value_type(ModalC0Basis{2}(V, 1)))
+end
 
 end # module