[wip] added indep_comp_getindex and fix FESpaces

of value with linked components. Indeed, the code previously made the assumption that all the components of the Number unknown of the FESpaces are unlinked to the other, which is false for the symmetric tensor types. A new getter for the independent components of the Number is added and used in the MultiValue'd FESpaces machinery. Also, the number of DoFs of the ReferenceFEs are given by num_indep_components(::Number). fixes gridap#923 gridap#908
Antoinemarteau · Aug 8, 2024 · 1543f7e · 1543f7e
1 parent 7769a3e
commit 1543f7e
Show file tree

Hide file tree

Showing 13 changed files with 188 additions and 57 deletions.
diff --git a/src/Adaptivity/FineToCoarseReferenceFEs.jl b/src/Adaptivity/FineToCoarseReferenceFEs.jl
@@ -45,8 +45,8 @@ function Arrays.evaluate!(cache,s::FineToCoarseDofBasis{T,<:LagrangianDofBasis},
   vals  = evaluate!(cf,field,b.nodes,s.child_ids)
   ndofs = length(b.dof_to_node)
   T2    = eltype(vals)
-  ncomps = num_components(T2)
-  @check ncomps == num_components(eltype(b.node_and_comp_to_dof)) """\n
+  ncomps = num_indep_components(T2) # use num_indep_components ?
+  @check ncomps == num_indep_components(eltype(b.node_and_comp_to_dof)) """\n
   Unable to evaluate LagrangianDofBasis. The number of components of the
   given Field does not match with the LagrangianDofBasis.
 
@@ -80,8 +80,8 @@ end
 
 
 """
-  Wrapper for a ReferenceFE which is specialised for 
-  efficiently evaluating FineToCoarseFields. 
+  Wrapper for a ReferenceFE which is specialised for
+  efficiently evaluating FineToCoarseFields.
 """
 struct FineToCoarseRefFE{T,D,A} <: ReferenceFE{D}
   reffe     :: T
@@ -132,4 +132,4 @@ function FESpaces.TestFESpace(model::DiscreteModel,rrules::AbstractVector{<:Refi
   basis, reffe_args, reffe_kwargs = reffe
   reffes = lazy_map(rr -> ReferenceFE(get_polytope(rr),rr,basis,reffe_args...;reffe_kwargs...),rrules)
   return TestFESpace(model,reffes;kwargs...)
-end
+end
diff --git a/src/Exports.jl b/src/Exports.jl
@@ -41,6 +41,7 @@ using Gridap.Arrays: ∑; export ∑
 @publish TensorValues outer
 @publish TensorValues diagonal_tensor
 @publish TensorValues num_components
+@publish TensorValues num_indep_components
 using Gridap.TensorValues: ⊙; export ⊙
 using Gridap.TensorValues: ⊗; export ⊗
 

diff --git a/src/FESpaces/CLagrangianFESpaces.jl b/src/FESpaces/CLagrangianFESpaces.jl
@@ -119,7 +119,7 @@ end
 # Helpers
 
 _default_mask(::Type) = true
-_default_mask(::Type{T}) where T <: MultiValue = ntuple(i->true,Val{length(T)}())
+_default_mask(::Type{T}) where T <: MultiValue = ntuple(i->true,Val{num_indep_components(T)}())
 
 _dof_type(::Type{T}) where T = T
 _dof_type(::Type{T}) where T<:MultiValue = eltype(T)
@@ -193,7 +193,7 @@ function _generate_node_to_dof_glue_component_major(
   z::MultiValue,node_to_tag,tag_to_masks)
   nfree_dofs = 0
   ndiri_dofs = 0
-  ncomps = length(z)
+  ncomps = num_indep_components(z)
   @check length(testitem(tag_to_masks)) == ncomps
   for (node,tag) in enumerate(node_to_tag)
     if tag == UNSET
@@ -218,7 +218,7 @@ function _generate_node_to_dof_glue_component_major(
   node_and_comp_to_dof = zeros(T,nnodes)
   nfree_dofs = 0
   ndiri_dofs = 0
-  m = zero(Mutable(T))
+  m = zeros(Int32, ncomps)
   for (node,tag) in enumerate(node_to_tag)
     if tag == UNSET
       for comp in 1:ncomps
@@ -245,7 +245,7 @@ function _generate_node_to_dof_glue_component_major(
         end
       end
     end
-    node_and_comp_to_dof[node] = m
+    node_and_comp_to_dof[node] = T(m...)
   end
   glue = NodeToDofGlue(
    free_dof_to_node,
@@ -301,7 +301,7 @@ function _generate_cell_dofs_clagrangian(
   cell_to_ctype,
   node_and_comp_to_dof)
 
-  ncomps = num_components(z)
+  ncomps = num_indep_components(z)
 
   ctype_to_lnode_to_comp_to_ldof = map(get_node_and_comp_to_dof,ctype_to_reffe)
   ctype_to_num_ldofs = map(num_dofs,ctype_to_reffe)
@@ -353,8 +353,8 @@ function _fill_cell_dofs_clagrangian!(
     p = cell_to_dofs.ptrs[cell]-1
     for (lnode, node) in enumerate(nodes)
       for comp in 1:ncomps
-        ldof = lnode_and_comp_to_ldof[lnode][comp]
-        dof = node_and_comp_to_dof[node][comp]
+        ldof = indep_comp_getindex(lnode_and_comp_to_ldof[lnode], comp)
+        dof =  indep_comp_getindex(node_and_comp_to_dof[node], comp)
         cell_to_dofs.data[p+ldof] = dof
       end
     end

diff --git a/src/FESpaces/ConstantFESpaces.jl b/src/FESpaces/ConstantFESpaces.jl
@@ -32,7 +32,7 @@ struct ConstantFESpace{V,T,A,B,C} <: SingleFieldFESpace
     cell_dof_basis_array = lazy_map(get_dof_basis,cell_reffe)
     cell_dof_basis = CellDof(cell_dof_basis_array,Triangulation(model),ReferenceDomain())
 
-    cell_dof_ids = Fill(Int32(1):Int32(num_components(field_type)),num_cells(model))
+    cell_dof_ids = Fill(Int32(1):Int32(num_indep_components(field_type)),num_cells(model))
     A = typeof(cell_basis)
     B = typeof(cell_dof_basis)
     C = typeof(cell_dof_ids)

diff --git a/src/Polynomials/JacobiPolynomialBases.jl b/src/Polynomials/JacobiPolynomialBases.jl
@@ -9,7 +9,7 @@ struct JacobiPolynomialBasis{D,T} <: AbstractVector{JacobiPolynomial}
   end
 end
 
-@inline Base.size(a::JacobiPolynomialBasis{D,T}) where {D,T} = (length(a.terms)*num_components(T),)
+@inline Base.size(a::JacobiPolynomialBasis{D,T}) where {D,T} = (length(a.terms)*num_indep_components(T),)
 @inline Base.getindex(a::JacobiPolynomialBasis,i::Integer) = JacobiPolynomial()
 @inline Base.IndexStyle(::JacobiPolynomialBasis) = IndexLinear()
 
@@ -49,7 +49,7 @@ return_type(::JacobiPolynomialBasis{D,T}) where {D,T} = T
 function return_cache(f::JacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(T,(np,ndof)))
   v = CachedArray(zeros(T,(ndof,)))
@@ -60,7 +60,7 @@ end
 function evaluate!(cache,f::JacobiPolynomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   r, v, c = cache
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -82,7 +82,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(V,xi)
   n = 1 + _maximum(f.orders)
@@ -101,7 +101,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -124,7 +124,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(gradient_type(V,xi),xi)
   n = 1 + _maximum(f.orders)
@@ -144,7 +144,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g, h = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -164,7 +164,7 @@ end
 # Optimizing evaluation at a single point
 
 function return_cache(f::JacobiPolynomialBasis{D,T},x::Point) where {D,T}
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   r = CachedArray(zeros(T,(ndof,)))
   xs = [x]
   cf = return_cache(f,xs)

diff --git a/src/Polynomials/ModalC0Bases.jl b/src/Polynomials/ModalC0Bases.jl
@@ -11,11 +11,12 @@ struct ModalC0Basis{D,T,V} <: AbstractVector{ModalC0BasisFunction}
     terms::Vector{CartesianIndex{D}},
     a::Vector{Point{D,V}},
     b::Vector{Point{D,V}}) where {D,T,V}
+
     new{D,T,V}(orders,terms,a,b)
   end
 end
 
-@inline Base.size(a::ModalC0Basis{D,T,V}) where {D,T,V} = (length(a.terms)*num_components(T),)
+@inline Base.size(a::ModalC0Basis{D,T,V}) where {D,T,V} = (length(a.terms)*num_indep_components(T),)
 @inline Base.getindex(a::ModalC0Basis,i::Integer) = ModalC0BasisFunction()
 @inline Base.IndexStyle(::ModalC0Basis) = IndexLinear()
 
@@ -101,7 +102,7 @@ return_type(::ModalC0Basis{D,T,V}) where {D,T,V} = T
 function return_cache(f::ModalC0Basis{D,T,V},x::AbstractVector{<:Point}) where {D,T,V}
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(T,(np,ndof)))
   v = CachedArray(zeros(T,(ndof,)))
@@ -112,7 +113,7 @@ end
 function evaluate!(cache,f::ModalC0Basis{D,T,V},x::AbstractVector{<:Point}) where {D,T,V}
   r, v, c = cache
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -134,7 +135,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(V,xi)
   n = 1 + _maximum(f.orders)
@@ -153,7 +154,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -176,7 +177,7 @@ function return_cache(
   f = fg.fa
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(gradient_type(V,xi),xi)
   n = 1 + _maximum(f.orders)
@@ -196,7 +197,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g, h = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -391,16 +392,17 @@ function _evaluate_nd_mc0!(
 end
 
 @inline function _set_value_mc0!(v::AbstractVector{V},s::T,k,l) where {V,T}
-  m = zero(Mutable(V))
+  ncomp = num_indep_components(V)
+  m = zeros(T,ncomp)
   z = zero(T)
-  js = eachindex(m)
+  js = 1:ncomp
   for j in js
     for i in js
       @inbounds m[i] = z
     end
     @inbounds m[j] = s
     i = k+l*(j-1)
-    @inbounds v[i] = m
+    @inbounds v[i] = V(m...)
   end
   k+1
 end
@@ -461,8 +463,12 @@ end
   k+1
 end
 
+# Indexing and m definition should be fixed if G contains symmetries, that is
+# if the code is  optimized for symmetric tensor V valued FESpaces
+# (if gradient_type(V) returned a symmetric higher order tensor type G)
 @inline function _set_gradient_mc0!(
   v::AbstractVector{G},s,k,l,::Type{V}) where {V,G}
+  @notimplementedif num_indep_components(G) != num_components(G) "Not implemented for symmetric Jacobian or Hessian"
 
   T = eltype(s)
   m = zero(Mutable(G))

diff --git a/src/Polynomials/MonomialBases.jl b/src/Polynomials/MonomialBases.jl
@@ -18,7 +18,7 @@ struct MonomialBasis{D,T} <: AbstractVector{Monomial}
   end
 end
 
-Base.size(a::MonomialBasis{D,T}) where {D,T} = (length(a.terms)*num_components(T),)
+Base.size(a::MonomialBasis{D,T}) where {D,T} = (length(a.terms)*num_indep_components(T),)
 # @santiagobadia : Not sure we want to create the monomial machinery
 Base.getindex(a::MonomialBasis,i::Integer) = Monomial()
 Base.IndexStyle(::MonomialBasis) = IndexLinear()
@@ -122,7 +122,7 @@ function return_cache(f::MonomialBasis{D,T},x::AbstractVector{<:Point}) where {D
   zxi = zero(eltype(eltype(x)))
   Tp = typeof( zT*zxi*zxi + zT*zxi*zxi  )
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(Tp,(np,ndof)))
   v = CachedArray(zeros(Tp,(ndof,)))
@@ -133,7 +133,7 @@ end
 function evaluate!(cache,f::MonomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
   r, v, c = cache
   np = length(x)
-  ndof = length(f.terms)*num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -157,7 +157,7 @@ function _return_cache(
   f = fg.fa
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   r = CachedArray(zeros(T,(np,ndof)))
   v = CachedArray(zeros(T,(ndof,)))
@@ -173,7 +173,7 @@ function _return_cache(
   TisbitsType::Val{false}) where {D,V,T}
 
   cache = _return_cache(fg,x,T,Val{true}())
-  z = CachedArray(zeros(eltype(T),D)) 
+  z = CachedArray(zeros(eltype(T),D))
   (cache...,z)
 end
 
@@ -197,7 +197,7 @@ function _evaluate!(
   r, v, c, g = cache
   z = zero(Mutable(VectorValue{D,eltype(T)}))
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -222,7 +222,7 @@ function _evaluate!(
   f = fg.fa
   r, v, c, g, z = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -255,7 +255,7 @@ function return_cache(
   f = fg.fa
   @check D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
-  ndof = length(f.terms)*num_components(V)
+  ndof = length(f)
   xi = testitem(x)
   T = gradient_type(gradient_type(V,xi),xi)
   n = 1 + _maximum(f.orders)
@@ -275,7 +275,7 @@ function evaluate!(
   f = fg.fa
   r, v, c, g, h = cache
   np = length(x)
-  ndof = length(f.terms) * num_components(T)
+  ndof = length(f)
   n = 1 + _maximum(f.orders)
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
@@ -406,15 +406,16 @@ function _evaluate_nd!(
 end
 
 function _set_value!(v::AbstractVector{V},s::T,k) where {V,T}
-  m = zero(Mutable(V))
+  ncomp = num_indep_components(V)
+  m = zeros(T,ncomp)
   z = zero(T)
-  js = eachindex(m)
+  js = 1:ncomp
   for j in js
     for i in js
       @inbounds m[i] = z
     end
     m[j] = s
-    v[k] = m
+    v[k] = V(m...)
     k += 1
   end
   k
@@ -440,7 +441,7 @@ function _gradient_nd!(
     _evaluate_1d!(c,x,orders[d],d)
     _gradient_1d!(g,x,orders[d],d)
   end
-  
+
   o = one(T)
   k = 1