Gridap.TensorValues

This module provide tensor types of supertype MultiValue that have very similar design and API to S/MArray from StaticArrays.jl, but with efficient storage for tensors with dependent components (e.g. symmetry). The main feature of this module is that MultiValue doesn't subtype AbstractArray, but Number!

This allows one to work with them as if they were scalar values in broadcasted operations on arrays of VectorValue objects (also for TensorValue or any <:MultiValue objects). For instance, one can perform the following manipulations:

# Assign a VectorValue to all the entries of an Array of VectorValues
A = zeros(VectorValue{2,Int}, (4,5))
v = VectorValue(12,31)
A .= v # This is possible since  VectorValue <: Number

# Broadcasting of tensor operations in arrays of TensorValues
t = TensorValue(13,41,53,17) # creates a 2x2 TensorValue
g = TensorValue(32,41,3,14) # creates another 2x2 TensorValue
B = fill(t,(1,5))
C = inner.(g,B) # inner product of g against all TensorValues in the array B
@show C
# C = [2494 2494 2494 2494 2494]

Exported names

AbstractSymTensorValue, HighOrderTensorValue, MultiValue, Mutable, QTensorValue, SkewSymTensorValue, SymFourthOrderTensorValue, SymTensorValue, SymTracelessTensorValue, TensorValue, ThirdOrderTensorValue, VectorValue, change_eltype, component_basis, congruent_prod, contracted_product, cross, data_index, det, diagonal_tensor, dot, double_contraction, eigen, indep_comp_getindex, indep_components_names, inner, inv, meas, mutable, n_components, norm, normalize, num_components, num_indep_components, outer, representatives_of_basis_dual, representatives_of_componentbasis_dual, skew_symmetric_part, symmetric_part, tr, ×, ⊗, ⊙, ⋅, ⋅², ⋅¹,

VectorValue{D,T} <: MultiValue{Tuple{D},T,1,D}

Type representing a first-order tensor, that is a vector, of length D.

TensorValue{D1,D2,T,L} <: MultiValue{Tuple{D1,D2},T,2,L}

Type representing a second-order D1×D2 tensor. It must hold L = D1*D2.

If only D1 or no dimension parameter is given to the constructor, D1=D2 is assumed.

SymTensorValue{D,T,L} <: AbstractSymTensorValue{D,T,L}

Type representing a symmetric second-order D×D tensor. It must hold L = D(D+1)/2.

It is constructed by providing the components of index (i,j) for 1 ≤ i ≤ j ≤ D.

SkewSymTensorValue{D,T,L} <: MultiValue{Tuple{D,D},T,2,L}

Type representing a skew symmetric second-order D×D tensor. It must hold L = D(D-1)/2.

It is constructed by providing the components of index (i,j) for 1 ≤ i < j ≤ D.

SymTracelessTensorValue{D,T,L} <: AbstractSymTensorValue{D,T,L}
QTensorValue{D,T,L}

Type representing a symetric second-order D×D tensor with zero trace. It must hold L = D(D+1)/2. This type is used to model the Q-tensor order parameter in nematic liquid cristals.

The constructor determines the value of index (D,D) as minus the sum of the other diagonal values, so it value musn't be provided. The constructor thus expects the L-1 components of indices (i,j) for 1 ≤ i ≤ D-1 and i ≤ j ≤ D.

ThirdOrderTensorValue{D1,D2,D3,T,L} = HighOrderTensorValue{Tuple{D1,D2,D3},T,3,L}

Alias representing a third-order D1×D2×D3 tensor, providing similar contructors than the 1st and 2nd order counterparts. It must hold L = D1*D2*D3.

If only D1 or no dimension parameter is given to the constructor, D1=D2=D3 is assumed.

SymFourthOrderTensorValue{D,T,L} <: MultiValue{Tuple{D,D,D,D},T,4,L}

Type representing a symmetric second-order D×D×D×D tensor, with symmetries ijkl↔jikl and ijkl↔ijlk. It must hold L = (D(D+1)/2)^2.

It is constructed by providing the components of index (i,j,k,l) for 1 ≤ i ≤ j ≤ D and 1 ≤ k ≤ l ≤ D.

struct HighOrderTensorValue{S,T,N,L} <: MultiValue{S,T,N,L}

Type representing a Nth tensor of size S with no dependant components. Thus, it always hold L = prod(S) components.

If Val(N) is given as first argument of the constructor, allequal(S) is assumed.

HighOrderTensorValue is a MultiValue wrapper for SArray{S,T,N,L} with N ≥ 3.

MultiValue{S,T,N,L} <: Number

Abstract type representing a multi-dimensional number value. The parameters are analog to that of StaticArrays.jl:

• S is a Tuple type holding the size of the tensor, e.g. Tuple{3} for a 3d vector or Tuple{2,4} for a 2 rows and 4 columns tensor,
• T is the type of the scalar components, should be subtype of Number,
• N is the order of the tensor, the length of S,
• L is the number of components stored internally.

MultiValues are immutable.

AbstractSymTensorValue{D,T,L} <: MultiValue{Tuple{D,D},T,2,L}

Abstract type representing any symmetric second-order D×D tensor, with symmetry ij↔ji.

See also SymTensorValue, SymTracelessTensorValue.

getindex(arg::MultiValue, inds...)
getindex(arg::MultiValue, i::Integer)
getindex(arg::MultiValue) = arg

MultiValues support a large subset of the indexing methods of AbstractArrays and StaticArrays, including inds argument that is a mixed tuples of Integer, CartesianIndex, slices and common range and array types.

The Number convention is used when no indices are provided: arg[] returns arg.

Examples

julia> t = TensorValue(1:9...)
TensorValue{3, 3, Int64, 9}(1, 2, 3, 4, 5, 6, 7, 8, 9)

julia> t[1,end]
7

julia> t[:,2]
VectorValue{3, Int64}(4, 5, 6)

julia> t[ isodd.(SMatrix(t))]
5-element Vector{Int64}:
 1
 3
 5
 7
 9

Extended help

Similarly to StaticArray.jl, the type of the returned value depend on inferability of the size from the type of the indices in inds.

Indeed, getindex returns

• a scalar (component) if inds contains only "scalars", i.e. Integer and CartesianIndex,
• a <:MultiValue tensor if inds additionally contains statically inferable index ranges/sets such as Colon()/:, SOneTo or StaticArray,
• an Array if inds additionally contains dynamic index ranges/sets such as UnitRange, OneTo or other types of AbstractArray.

warning

Indexing methods that do not return a scalars will loose any symmetry or other known internal constraint, as in SymTensorValue(1:6...)[SOneTo(2),SOneTo(2)]

julia> t[SOneTo(2),SOneTo(3)]
TensorValue{2, 3, Int64, 6}(1, 2, 4, 5, 7, 8)

julia> t[Base.OneTo(2),Base.OneTo(3)]
2×3 Matrix{Int64}:
 1  4  7
 2  5  8

julia> t[MVector(1,2),CartesianIndex(3)]
VectorValue{2, Int64}(7, 8)

julia> t[1:2,CartesianIndex(3)]
2-element Vector{Int64}:
 7
 8


julia> mask = [true false false; false true false; false false true]
3×3 Matrix{Bool}:
 1  0  0
 0  1  0
 0  0  1

julia> t[mask]
3-element Vector{Int64}:
 1
 5
 9

num_components(::Type{<:Number})
num_components(a::Number)

Total number of components of a Number or MultiValue, that is 1 for scalars and the product of the size dimensions for a MultiValue. This is the same as length.

mutable(a::MultiValue)

Converts a into a mutable array of type MArray defined by StaticArrays.jl.

See also Mutable.

Mutable(T::Type{<:MultiValue}) -> ::Type{<:MArray}
Mutable(a::MultiValue)

Return the concrete mutable MArray type (defined by StaticArrays.jl) corresponding to the MultiValue type T or array size and type of a.

See also mutable.

num_indep_components(::Type{<:Number})
num_indep_components(a::Number)

Number of independent components of a Number, that is num_components minus the number of components determined from others by symmetries or constraints.

For example, a TensorValue{3,3} has 9 independent components, a SymTensorValue{3} has 6 and a SymTracelessTensorValue{3} has 5. But they all have 9 (non independent) components.

indep_comp_getindex(a::Number,i)

Get the ith independent component of a. It only differs from getindex(a,i) when the components of a are interdependent, see num_indep_components. i should be in 1:num_indep_components(a).

component_basis(V::Type{<:MultiValue}) -> V[ Vᵢ... ]
component_basis(T::Type{<:Real}) -> [ one(T) ]
component_basis(a::T<:Number)

Given a Number type V with N independent components, return a vector of N values $\{ Vᵢ=V(eᵢ) \}_i$ forming the component basis of $\{ u : u\text{ isa }V\}$, where $\{eᵢ\}_i$ is the Cartesian basis of (eltype(V))ᴺ.

The Vᵢ verify the property that for any u::V,

u = sum( indep_comp_getindex(u,i)*Vᵢ for i ∈ 1:N )

representatives_of_componentbasis_dual(V::Type{<:MultiValue}) -> V[ Vᵢ... ]
representatives_of_componentbasis_dual(T::Type{<:Real}) -> [ one(T) ]
representatives_of_componentbasis_dual(a::V<:Number)

Given a Number type V with N independent components, return a vector of N tensors $\{ Vⁱ \}_i$ that define the form basis $\{ Lⁱ := (u → u ⊙ Vⁱ) \}_i$ that is the dual of the component basis $\{ Vᵢ=V(eᵢ) \}_i$ (where $\{eᵢ\}_i$ is the Cartesian basis of (eltype(V))ᴺ).

The Vᵢ/Vʲ verify the kronecker delta property $Vᵢ⊙Vʲ = δᵢʲ$, and for any u::V,

u = V( [ Lⁱ(u) for i ∈ 1:N ]... )
  = V( [ u⊙Vⁱ  for i ∈ 1:N ]... )

Rq, when V has dependent components, the Vⁱ are not a component basis because $Vʲ ≠ Vᵢ$, $Vᵢ⊙Vⱼ≠δᵢʲ$ and

u ≠ sum( indep_comp_getindex(u,i)*Vⁱ for i ∈ 1:N )

See also representatives_of_basis_dual to compute a dual subspace basis.

representatives_of_basis_dual(basis)

Same as representatives_of_componentbasis_dual, but takes a basis as arguement, i.e. a collection of basis vectors – in the linear algebra sense – of type V<:Number, so basis spans a subpace of Span{component_basis(V)...}.

Computing the dual basis to a subspace basis is usefull for e.g. value types G resulting from applying a derivative to function of value type Vd having dependent components:

Vd = SymTensorValue{2,Float64} # has 3 independent components
P = Point{2,Float64}
G = Fields.gradient_type(Vd, zero(P)) # ThirdOrderTensorValue{2, 2, 2, Float64}

# Span{ component_basis(G)... } is of dimension 8, but the gradient of a Vd
# valued function lives in the 6-dimensional space of basis:
grad_Vd_basis = [ pi⊗vi for pi in component_basis(P) for vi in component_basis(Vd)]

# of dual basis { v -> v ⊙ Gdᵢ}ᵢ for Gdᵢ in:
grad_Vd_dual_basis = representatives_of_basis_dual(grad_Vd_basis)

# kronecker delta property
[ gi ⊙ gdi for gi in grad_Vd_basis, gdi in grad_Vd_dual_basis] ≈ LinearAlgebra.I # true

indep_components_names(::MultiValue)
indep_components_names(::Type{<:MultiValue{S}})

Return an array of strings containing the component labels in the order they are exported in VTK file.

If all dimensions of the tensor shape S are smaller than 3, the components are named with letters "X","Y" and "Z" similarly to the automatic naming of Paraview. Else, if maximum(S)>3, they are labeled by integers starting from "1".

change_eltype(m::Number,::Type{T2})
change_eltype(M::Type{<:Number},::Type{T2})

For multivalues, returns M or typeof(m) but with the component type (MultiValue's parametric type T) changed to T2.

For scalars (or any non MultiValue number), change_eltype returns T2.

inner(a::MultiValue{S}, b::MultiValue{S}) -> scalar
a ⊙ b

Inner product of two tensors, that is the full contraction along each indices. The size S of a and b must match.

dot(a::MultiValue{Tuple{...,D}}, b::MultiValue{Tuple{D,...}})
a ⋅¹ b
a ⋅ b

Single contraction of two tensors a and b, of the last index of a with the first index of b. The corresponding dimensions D must match. No symmetry is preserved. On two vectors, this is the same as the inner product.

double_contraction(a::MultiValue{Tuple{...,D,E}}, b::MultiValue{Tuple{D,E,...})
a ⋅² b

Double contraction of two tensors a and b, along the two last indices of a and two first of b. The corresponding dimensions D and E must match, the contraction order is chosen to be consistent with the inner product of second order tensors.

The double_contraction between second- and/or fourth-order symmetric tensors preserves the symmetry (returns a symmetric tensor type).

outer(a,b)
a ⊗ b

Outer product (or tensor-product) of two Numbers and/or MultiValues, that is (a⊗b)[i₁,...,iₙ,j₁,...,jₙ] = a[i₁,...,iₙ]*b[j₁,...,jₙ]. This falls back to standard multiplication if a or b is a scalar.

normalize(a::MultiValue)

Normalization of a tensor value.

contracted_product(Val(n), a::MultiValue, b::MultiValue)

Compute the tensor c that results from the tensor-product of a and b contracted over the n last indices of a and n first indices of b:

c[i₁,...,iₘ₋ₙ,jₙ₊₁,...,jₚ] = Σ_{l₁,...,lₙ} a[i₁,...,iₘ₋ₙ,l₁,...,lₙ]*b[l₁,...,lₙ,jₙ₊₁,...,jₚ]

where m is the order a and p that of b. The sum runs for each cartesian indices (l₁,...,lₙ) of the last n dimension of a, which must match the first n dimensions of b.

This operation generalizes:

• outer / ⊗ for n=0,
• dot / ⋅ for n=1,
• double_contraction / ⋅² for n=2,
• inner / ⊙ for n = length(size(a)) = length(size(b)).

This function is not optimized for tensor types with dependent components (use the specific functions above if possible), but is used as default generic implementation.

det(a::MultiValue{Tuple{D,D},T})

Determinent of square second order tensors.

inv(a::MultiValue{Tuple{D,D}})

Inverse of a second order tensor.

eigen(a::MultiValue{Tuple{D,D}})

Eigenvalue decomposition of a square second order tensor.

symmetric_part(v::MultiValue{Tuple{D,D}})::AbstractSymTensorValue

Return the symmetric part of second order tensor, that is ½(v + vᵀ). Return v if v isa AbstractSymTensorValue, and the zero symmetric tensor if v isa SkewSymTensorValue.

skew_symmetric_part(v::MultiValue{Tuple{D,D}})::SkewSymTensorValue{D}

Return the asymmetric part of v, that is ½(v - vᵀ). Return the zero skew symmetric tensor if v isa AbstractSymTensorValue, and v itself if v isa SkewSymTensorValue.

congruent_prod(a, b)

Given a square second order tensors a and b, return bᵀ⋅a⋅b. The type of the resulting value is (skew) symmetric stable w.r.t. typeof(a).

diagonal_tensor(v::VectorValue{D,T}) -> ::TensorValue{D,D,T}

Return a diagonal D×D tensor with diagonal containing the elements of v.

tr(v::MultiValue{Tuple{D,D}})

Return the trace of a second order square tensor, defined by Σᵢ vᵢᵢ or 0 if D=0.

tr(v::MultiValue{Tuple{D,D,D2}}) -> ::VectorValue{D2}

Return a vector of length D2 of traces computed on the first two indices: resⱼ = Σᵢ vᵢᵢⱼ.

norm(u::MultiValue{Tuple{D}})
norm(u::MultiValue{Tuple{D1,D2}})

Euclidean (2-)norm of u, namely sqrt(inner(u,u)).

meas(a::MultiValue{Tuple{D}})
meas(a::MultiValue{Tuple{1,D2}})

Euclidean norm of a vector.

meas(J::MultiValue{Tuple{D1,D2}})

Returns the absolute D1-dimensional volume of the parallelepiped formed by the rows of J, that is sqrt(det(J⋅Jᵀ)), or abs(det(J)) if D1=D2. This is used to compute the contribution of the Jacobian matrix J of a changes of variables in integrals.

cross(a::VectorValue{3}, b::VectorValue{3}) -> VectorValue{3}
cross(a::VectorValue{2}, b::VectorValue{2}) -> scalar
a × b

Cross product of 2D and 3D vector.

one(::SymFourthOrderTensorValue{D,T}})

Returns the tensor resᵢⱼₖₗ = δᵢₖδⱼₗ(δᵢⱼ + (1-δᵢⱼ)/2).

The scalar type T2 of the result is typeof(one(T)/2).

Previously used to transform Cartesian indices to linear indices that index MultiValue's private internal storage.