Use Symmetric Covariant Axis2Tensor

wip

Author: charleskawczynski

src/Geometry/Geometry.jl

@@ -23,6 +23,65 @@ export Contravariant1Vector,
     Contravariant123Vector
 
 
+"""
+    SimpleSymmetric
+
+A simple Symmetric matrix. `LinearAlgebra.Symmetric` has a field
+`uplo::Char`, which results in a failed `check_basetype(T, S)`.
+"""
+struct SimpleSymmetric{T, S <: AbstractMatrix{<:T}} <: AbstractMatrix{T}
+    data::S
+    function SimpleSymmetric{T, S}(data) where {T, S <: AbstractMatrix{<:T}}
+        new{T, S}(data)
+    end
+end
+SimpleSymmetric(A) = cc_symmetric_type(typeof(A))(A)
+
+LinearAlgebra.transpose(A::SimpleSymmetric) = A
+
+"""
+    cc_symmetric_type(T::Type)
+
+The type of the object returned by `symmetric(::T, ::Symbol)`. For matrices, this is an
+appropriately typed `Symmetric`, for `Number`s, it is the original type. If `symmetric` is
+implemented for a custom type, so should be `cc_symmetric_type`, and vice versa.
+"""
+cc_symmetric_type(::Type{T}) where {S, T <: AbstractMatrix{S}} =
+    SimpleSymmetric{
+        Union{S, Base.promote_op(transpose, S), cc_symmetric_type(S)},
+        T,
+    }
+cc_symmetric_type(::Type{T}) where {S <: Number, T <: AbstractMatrix{S}} =
+    SimpleSymmetric{S, T}
+cc_symmetric_type(
+    ::Type{T},
+) where {S <: AbstractMatrix, T <: AbstractMatrix{S}} =
+    SimpleSymmetric{AbstractMatrix, T}
+cc_symmetric_type(::Type{T}) where {T <: Number} = T
+
+Base.@propagate_inbounds function Base.getindex(
+    A::SimpleSymmetric,
+    i::Integer,
+    j::Integer,
+)
+    @boundscheck checkbounds(A, i, j)
+    @inbounds if i == j
+        data_ij = A.data[i, j]
+        if data_ij isa AbstractMatrix
+            return SimpleSymmetric(data_ij)::cc_symmetric_type(eltype(A.data))
+        elseif data_ij isa Number
+            return data_ij::cc_symmetric_type(eltype(A.data))
+        end
+    elseif (i < j)
+        return A.data[i, j]
+    else
+        return LinearAlgebra.transpose(A.data[j, i])
+    end
+end
+
+Base.@propagate_inbounds Base.getindex(A::SimpleSymmetric, i::Integer) =
+    A.data[i]
+
 
 include("coordinates.jl")
 include("axistensors.jl")

src/Geometry/axistensors.jl

@@ -136,7 +136,10 @@ struct AxisTensor{
     T,
     N,
     A <: NTuple{N, AbstractAxis},
-    S <: StaticArray{<:Tuple, T, N},
+    S <: Union{
+        SimpleSymmetric{T, <:StaticArray{<:Tuple, T, N}},
+        StaticArray{<:Tuple, T, N},
+    },
 } <: AbstractArray{T, N}
     axes::A
     components::S
@@ -147,7 +150,10 @@ AxisTensor(
     components::S,
 ) where {
     A <: Tuple{Vararg{AbstractAxis}},
-    S <: StaticArray{<:Tuple, T, N},
+    S <: Union{
+        SimpleSymmetric{T, <:StaticArray{<:Tuple, T, N}},
+        StaticArray{<:Tuple, T, N},
+    },
 } where {T, N} = AxisTensor{T, N, A, S}(axes, components)
 
 AxisTensor(axes::Tuple{Vararg{AbstractAxis}}, components) =

src/Geometry/localgeometry.jl

@@ -1,8 +1,18 @@
+import LinearAlgebra
+import .Geometry: SimpleSymmetric
+
+import InteractiveUtils
 import LinearAlgebra: issymmetric
 
 isapproxsymmetric(A::AbstractMatrix{T}; rtol = 10 * eps(T)) where {T} =
     Base.isapprox(A, A'; rtol)
 
+function SimpleSymmetric(x::AxisTensor{FT}) where {FT}
+    c = SimpleSymmetric(components(x))
+    A = axes(x)
+    return AxisTensor{FT, ndims(x), typeof(A), typeof(c)}(A, c)
+end
+
 """
     LocalGeometry
 
@@ -22,9 +32,13 @@ struct LocalGeometry{I, C <: AbstractPoint, FT, S}
     "Partial derivatives of the map from `x` to `ξ`: `∂ξ∂x[i,j]` is ∂ξⁱ/∂xʲ"
     ∂ξ∂x::Axis2Tensor{FT, Tuple{ContravariantAxis{I}, LocalAxis{I}}, S}
     "Contravariant metric tensor (inverse of gᵢⱼ), transforms covariant to contravariant vector components"
-    gⁱʲ::Axis2Tensor{FT, Tuple{ContravariantAxis{I}, ContravariantAxis{I}}, S}
+    gⁱʲ::Axis2Tensor{FT, Tuple{ContravariantAxis{I}, ContravariantAxis{I}}, SimpleSymmetric{FT, S}}
     "Covariant metric tensor (gᵢⱼ), transforms contravariant to covariant vector components"
-    gᵢⱼ::Axis2Tensor{FT, Tuple{CovariantAxis{I}, CovariantAxis{I}}, S}
+    gᵢⱼ::Axis2Tensor{
+        FT,
+        Tuple{CovariantAxis{I}, CovariantAxis{I}},
+        SimpleSymmetric{FT, S},
+    }
     @inline function LocalGeometry(
         coordinates,
         J,
@@ -34,15 +48,16 @@ struct LocalGeometry{I, C <: AbstractPoint, FT, S}
         ∂ξ∂x = inv(∂x∂ξ)
         C = typeof(coordinates)
         Jinv = inv(J)
-        gⁱʲ = ∂ξ∂x * ∂ξ∂x'
-        gᵢⱼ = ∂x∂ξ' * ∂x∂ξ
-        isapproxsymmetric(components(gⁱʲ)) || error("gⁱʲ is not symmetric.")
-        isapproxsymmetric(components(gᵢⱼ)) || error("gᵢⱼ is not symmetric.")
+        gⁱʲ₀ = ∂ξ∂x * ∂ξ∂x'
+        gᵢⱼ₀ = ∂x∂ξ' * ∂x∂ξ
+        isapproxsymmetric(components(gⁱʲ₀)) || error("gⁱʲ is not symmetric.")
+        isapproxsymmetric(components(gᵢⱼ₀)) || error("gᵢⱼ is not symmetric.")
+        gⁱʲ = SimpleSymmetric(gⁱʲ₀)
+        gᵢⱼ = SimpleSymmetric(gᵢⱼ₀)
         return new{I, C, FT, S}(coordinates, J, WJ, Jinv, ∂x∂ξ, ∂ξ∂x, gⁱʲ, gᵢⱼ)
     end
 end
 
-
 """
     SurfaceGeometry
 

test/Geometry/axistensor_conversion_benchmarks.jl

@@ -1,3 +1,7 @@
+#=
+julia --project
+using Revise; include(joinpath("test", "Geometry", "axistensor_conversion_benchmarks.jl"))
+=#
 using Test, StaticArrays
 #! format: off
 import Random, BenchmarkTools, StatsBase,