Merge pull request #1041 from Antoinemarteau/q_tensor

Q_tensor and various TensorValues improvements.

NEWS.md

@@ -12,6 +12,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added MacroFElements. These are defined as having the basis/dof-basis of a FESpace created on top of a RefinementRule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024).
 - Added Barycentric refinement rule in 2D and 3D. Added Simplexify refinement rule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024).
 - Added names to vector and tensor components in VTK exports, to avoid Paraview's automatic (sometimes wrong) guesses. See `TensorValues.indep_components_names`. Since PR[#1038](https://github.com/gridap/Gridap.jl/pull/1038).
+- Misc improvements of the `TensorValues` module:  See `TensorValues.indep_components_names`. Since PR[#1040](https://github.com/gridap/Gridap.jl/pull/1040).
+    - Documented all symbols exported by the module
+    - Improved and added test for some API function of `MultiValue` (general `diag` of 2nd order tensors, fixed `convert` of 3rd order tensors to SArray, avoid unwanted fallback of `num_components` on `MultiValue` types with undefined dimensions, more autodiff tests, better `double_contraction` API (prevent invalid operation giving indexing errors and enable valid operations)).
+    - Added a clear separation between the physical components access (`getindex`, `num_components`) and the numerical access to the stored independent components (`num_indep_components`, `indep_comp_getindex`) to enable using symmetric tensor types as unknown in FE Spaces.
+    - Implemented automatic differentiation `gradient` and `laplacian` for second order tensor, and `divergence` for third order tensors.
+    - Added `AbstractSymTensorValue`, an abstract type for second order symmetric tensors, and `SymTracelessTensorValue` (aliased to `QTensorValue`), a type for traceless symmetric tensors. `SymTensorValue` is now subtype of `AbstractSymTensorValue`.
+    - A convergence test for Poisson problem of `QTensorValue` unknown field validates the implementation.
 
 ## [0.18.6] - 2024-08-29
 

src/Adaptivity/FineToCoarseReferenceFEs.jl

@@ -44,8 +44,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)
+  @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.
   If you are trying to interpolate a function on a FESpace make sure that
@@ -77,8 +77,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
@@ -129,4 +129,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

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 ⊗
 

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
@@ -215,10 +215,10 @@ function _generate_node_to_dof_glue_component_major(
   diri_dof_to_tag = ones(Int8,ndiri_dofs)
   T = change_eltype(z,Int32)
   nnodes = length(node_to_tag)
-  node_and_comp_to_dof = zeros(T,nnodes)
+  node_and_comp_to_dof = Vector{T}(undef,nnodes)
   nfree_dofs = 0
   ndiri_dofs = 0
-  m = zero(Mutable(T))
+  m = zero(MVector{ncomps, Int32})
   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] = Tuple(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

src/FESpaces/ConstantFESpaces.jl

@@ -5,46 +5,42 @@ struct ConstantFESpace <: SingleFieldFESpace
 end
 """
 struct ConstantFESpace{V,T,A,B,C} <: SingleFieldFESpace
-model::DiscreteModel
-cell_basis::A
-cell_dof_basis::B
-cell_dof_ids::C
-function ConstantFESpace(model;
-                       vector_type::Type{V}=Vector{Float64},
-                       field_type::Type{T}=Float64) where {V,T}
-function setup_cell_reffe(model::DiscreteModel,
-  reffe::Tuple{<:ReferenceFEName,Any,Any}; kwargs...)
-  basis, reffe_args,reffe_kwargs = reffe
-  cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
-end
-reffe=ReferenceFE(lagrangian,T,0)
-cell_reffe = setup_cell_reffe(model,reffe)
-cell_basis_array=lazy_map(get_shapefuns,cell_reffe)
-
-cell_basis=SingleFieldFEBasis(
-  cell_basis_array,
-  Triangulation(model),
-  TestBasis(),
-  ReferenceDomain())
-
-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))
-A=typeof(cell_basis)
-B=typeof(cell_dof_basis)
-C=typeof(cell_dof_ids)
-new{V,T,A,B,C}(model,
-               cell_basis,
-               cell_dof_basis,
-               cell_dof_ids)
-end
+  model::DiscreteModel
+  cell_basis::A
+  cell_dof_basis::B
+  cell_dof_ids::C
+
+  function ConstantFESpace(model;
+                         vector_type::Type{V}=Vector{Float64},
+                         field_type::Type{T}=Float64) where {V,T}
+    function setup_cell_reffe(model::DiscreteModel,
+        reffe::Tuple{<:ReferenceFEName,Any,Any}; kwargs...)
+      basis, reffe_args,reffe_kwargs = reffe
+      cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
+    end
+
+    reffe = ReferenceFE(lagrangian,T,0)
+    cell_reffe = setup_cell_reffe(model,reffe)
+    cell_basis_array = lazy_map(get_shapefuns,cell_reffe)
+
+    cell_basis = SingleFieldFEBasis(
+      cell_basis_array,
+      Triangulation(model),
+      TestBasis(),
+      ReferenceDomain())
+
+    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_indep_components(field_type)),num_cells(model))
+    A = typeof(cell_basis)
+    B = typeof(cell_dof_basis)
+    C = typeof(cell_dof_ids)
+    new{V,T,A,B,C}(model, cell_basis, cell_dof_basis, cell_dof_ids)
+  end
 end
 
-# Genuine functions
-function TrialFESpace(f::ConstantFESpace)
-f
-end
+TrialFESpace(f::ConstantFESpace) = f
 
 # Delegated functions
 get_triangulation(f::ConstantFESpace) = Triangulation(f.model)
@@ -70,27 +66,24 @@ num_dirichlet_tags(f::ConstantFESpace) = 0
 get_dirichlet_dof_tag(f::ConstantFESpace) = Int8[]
 
 function scatter_free_and_dirichlet_values(f::ConstantFESpace,fv,dv)
-cell_dof_ids = get_cell_dof_ids(f)
-lazy_map(Broadcasting(PosNegReindex(fv,dv)),cell_dof_ids)
+  cell_dof_ids = get_cell_dof_ids(f)
+  lazy_map(Broadcasting(PosNegReindex(fv,dv)),cell_dof_ids)
 end
 
-function gather_free_and_dirichlet_values!(free_vals,
-                                                      dirichlet_vals,
-                                                      f::ConstantFESpace,
-                                                      cell_vals)
-cell_dofs = get_cell_dof_ids(f)
-cache_vals = array_cache(cell_vals)
-cache_dofs = array_cache(cell_dofs)
-cells = 1:length(cell_vals)
-
-_free_and_dirichlet_values_fill!(
-free_vals,
-dirichlet_vals,
-cache_vals,
-cache_dofs,
-cell_vals,
-cell_dofs,
-cells)
-
-(free_vals,dirichlet_vals)
+function gather_free_and_dirichlet_values!(free_vals, dirichlet_vals, f::ConstantFESpace, cell_vals)
+  cell_dofs = get_cell_dof_ids(f)
+  cache_vals = array_cache(cell_vals)
+  cache_dofs = array_cache(cell_dofs)
+  cells = 1:length(cell_vals)
+
+  _free_and_dirichlet_values_fill!(
+    free_vals,
+    dirichlet_vals,
+    cache_vals,
+    cache_dofs,
+    cell_vals,
+    cell_dofs,
+    cells)
+
+  (free_vals,dirichlet_vals)
 end

src/FESpaces/FESpaces.jl

@@ -10,6 +10,7 @@ using Test
 using FillArrays
 using SparseArrays
 using LinearAlgebra
+using StaticArrays
 using ForwardDiff
 
 using Gridap.Helpers

src/Fields/AutoDiff.jl

@@ -59,6 +59,13 @@ function gradient(f::Function,x::Point,fx::VectorValue)
   TensorValue(transpose(ForwardDiff.jacobian(y->get_array(f(y)),get_array(x))))
 end
 
+# Implementation for all second order tensor values
+# Does not exploit possible symmetries
+function gradient(f::Function,x::Point{A},fx::S) where S<:MultiValue{Tuple{B,C}} where {A,B,C}
+  a = transpose(ForwardDiff.jacobian(y->get_array(f(y)),get_array(x)))
+  ThirdOrderTensorValue(reshape(a, (A,B,C)))
+end
+
 function gradient(f::Function,x::Point,fx::MultiValue)
   @notimplemented
 end
@@ -71,10 +78,20 @@ function divergence(f::Function,x::Point)
   divergence(f,x,return_value(f,x))
 end
 
-function divergence(f::Function,x::Point,fx)
+function divergence(f::Function,x::Point,fx::VectorValue)
   tr(gradient(f,x,fx))
 end
 
+function divergence(f::Function,x::Point{D},fx::S) where S<:MultiValue{Tuple{D,A},T} where {D,A,T}
+  a = ForwardDiff.jacobian(y->get_array(f(y)),get_array(x))
+  VectorValue{A,T}( ntuple(k -> sum(i-> a[(k-1)*D+i,i], 1:D),A) )
+end
+
+function divergence(f::Function,x::Point{D},fx::S) where S<:MultiValue{Tuple{D,A,B},T} where {D,A,B,T}
+  a = ForwardDiff.jacobian(y->get_array(f(y)),get_array(x))
+  TensorValue{A,B,T}( ntuple(k -> sum(i-> a[(k-1)*D+i,i], 1:D),A*B) )
+end
+
 function divergence(f::Function,x::Point,fx::TensorValue{2,2})
   g(x) = SVector(f(x).data)
   a = ForwardDiff.jacobian(g,get_array(x))
@@ -94,7 +111,7 @@ function divergence(f::Function,x::Point,fx::TensorValue{3,3})
    )
 end
 
-function divergence(f::Function,x::Point,fx::SymTensorValue{2})
+function divergence(f::Function,x::Point{2},fx::AbstractSymTensorValue{2})
   g(x) = SVector(f(x).data)
   a = ForwardDiff.jacobian(g,get_array(x))
   VectorValue(
@@ -103,7 +120,7 @@ function divergence(f::Function,x::Point,fx::SymTensorValue{2})
   )
 end
 
-function divergence(f::Function,x::Point,fx::SymTensorValue{3})
+function divergence(f::Function,x::Point{3},fx::AbstractSymTensorValue{3})
   g(x) = SVector(f(x).data)
   a = ForwardDiff.jacobian(g,get_array(x))
   VectorValue(
@@ -113,6 +130,10 @@ function divergence(f::Function,x::Point,fx::SymTensorValue{3})
    )
 end
 
+function divergence(f::Function,x::Point,fx::MultiValue)
+  @notimplemented
+end
+
 function curl(f::Function,x::Point)
   grad2curl(gradient(f,x))
 end
@@ -125,11 +146,25 @@ function laplacian(f::Function,x::Point,fx::Number)
   tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(f,y), get_array(x)))
 end
 
-function laplacian(f::Function,x::Point,fx::VectorValue)
-  A = length(x)
-  B = length(fx)
-  a = ForwardDiff.jacobian(y->transpose(ForwardDiff.jacobian(z->get_array(f(z)),y)), get_array(x))
-  tr(ThirdOrderTensorValue{A,A,B}(Tuple(transpose(a))))
+function laplacian(f::Function,x::Point{A},fx::VectorValue{B,T}) where {A,B,T}
+  a = MMatrix{A*B,A,T}(undef)
+  ForwardDiff.jacobian!(a, y->ForwardDiff.jacobian(z->get_array(f(z)),y), get_array(x))
+  VectorValue{B,T}( ntuple(k -> sum(i-> a[(i-1)*B+k,i], 1:A),B) )
+end
+
+function laplacian(f::Function,x::Point{A},fx::S) where S<:MultiValue{Tuple{B,C},T} where {A,B,C,T}
+  a = MMatrix{A*B*C,A,T}(undef)
+  ForwardDiff.jacobian!(a, y->ForwardDiff.jacobian(z->get_array(f(z)),y), get_array(x))
+  t = ntuple(k -> sum(i-> a[(i-1)*B*C+k,i], 1:A),B*C)
+  S(SMatrix{B,C}(t)) #Necessary cast to build e.g. symmetric tensor values
+end
+
+# Implementation for any third order tensor values
+function laplacian(f::Function,x::Point{A},fx::S) where S<:MultiValue{Tuple{B,C,D},T} where {A,B,C,D,T}
+  a = MMatrix{A*B*C*D,A,T}(undef)
+  ForwardDiff.jacobian!(a, y->ForwardDiff.jacobian(z->get_array(f(z)),y), get_array(x))
+  t = ntuple(k -> sum(i-> a[(i-1)*B*C*D+k,i], 1:A),B*C*D)
+  S(SArray{Tuple{B,C,D}}(t))
 end
 
 function laplacian(f::Function,x::Point,fx::MultiValue)

src/ODEs/TimeDerivatives.jl

@@ -121,14 +121,9 @@ function _time_derivative(T::Type{<:Real}, f, t, x)
   ForwardDiff.derivative(partial, t)
 end
 
-function _time_derivative(T::Type{<:VectorValue}, f, t, x)
+function _time_derivative(T::Type{<:MultiValue}, f, t, x)
   partial(t) = get_array(f(t)(x))
-  VectorValue(ForwardDiff.derivative(partial, t))
-end
-
-function _time_derivative(T::Type{<:TensorValue}, f, t, x)
-  partial(t) = get_array(f(t)(x))
-  TensorValue(ForwardDiff.derivative(partial, t))
+  T(ForwardDiff.derivative(partial, t))
 end
 
 ##########################################