wip

Author: charleskawczynski

src/Geometry/Geometry.jl

@@ -22,81 +22,7 @@ export Contravariant1Vector,
     Contravariant23Vector,
     Contravariant123Vector
 
-import LinearAlgebra
-"""
-    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
-Base.similar(::Type{SimpleSymmetric{T, S}}, ax) where {T, S} = SimpleSymmetric{T}(similar(S))
-
-Base.size(A::SimpleSymmetric) = Base.size(A.data)
-Base.eltype(::SimpleSymmetric{T}) where {T} = T
-# Base.iterate(iter, state)
-Base.ndims(::SimpleSymmetric{T, S}) where {T, S} = Base.ndims(S)
-# Base.:*(A::SimpleSymmetric, B::SimpleSymmetric) = A * copyto!(similar(parent(B)), B)
-
-# Base.:*(A::SimpleSymmetric, x::Number) = SimpleSymmetric(A.data*x)
-# Base.:*(x::Number, A::SimpleSymmetric) = SimpleSymmetric(x*A.data)
-# Base.:*(A::SimpleSymmetric, B::AbstractMatrix) = A.data * B
-# Base.:*(A::AbstractMatrix, B::SimpleSymmetric) = A * B.data
-# Base.:/(A::SimpleSymmetric, x::Number) = SimpleSymmetric(A.data / x)
-# Base.:/(A::SimpleSymmetric, B::AbstractMatrix) = A.data * B
-# Base.:/(A::AbstractMatrix, B::SimpleSymmetric) = A * B.data
-
-"""
-    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("simple_symmetric.jl")
 include("coordinates.jl")
 include("axistensors.jl")
 include("localgeometry.jl")

src/Geometry/localgeometry.jl

@@ -18,7 +18,7 @@ end
 
 The necessary local metric information defined at each node.
 """
-struct LocalGeometry{I, C <: AbstractPoint, FT, S}
+struct LocalGeometry{I, C <: AbstractPoint, FT, S, L}
     "Coordinates of the current point"
     coordinates::C
     "Jacobian determinant of the transformation `ξ` to `x`"
@@ -32,12 +32,18 @@ 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}}, SimpleSymmetric{FT, S}}
+    gⁱʲ::Axis2Tensor{
+        FT,
+        Tuple{ContravariantAxis{I}, ContravariantAxis{I}},
+        # SimpleSymmetric{FT, S},
+        SimpleSymmetric{2, FT, L},
+    }
     "Covariant metric tensor (gᵢⱼ), transforms contravariant to covariant vector components"
     gᵢⱼ::Axis2Tensor{
         FT,
         Tuple{CovariantAxis{I}, CovariantAxis{I}},
-        SimpleSymmetric{FT, S},
+        # SimpleSymmetric{FT, S},
+        SimpleSymmetric{2, FT, L},
     }
     @inline function LocalGeometry(
         coordinates,
@@ -54,7 +60,8 @@ struct LocalGeometry{I, C <: AbstractPoint, FT, S}
         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ᵢⱼ)
+        L = triangular_nonzeros(S)
+        return new{I, C, FT, S, L}(coordinates, J, WJ, Jinv, ∂x∂ξ, ∂ξ∂x, gⁱʲ, gᵢⱼ)
     end
 end