Use a Symmetric matrix in the Covariant Axis2Tensor

.buildkite/pipeline.yml

@@ -180,6 +180,10 @@ steps:
         key: unit_rmul_with_projection
         command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/rmul_with_projection.jl"
 
+      - label: "Unit: simple_symmetric"
+        key: unit_simple_symmetric
+        command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/unit_simple_symmetric.jl"
+
   - group: "Unit: Meshes"
     steps:
 

src/Geometry/Geometry.jl

@@ -22,8 +22,7 @@ export Contravariant1Vector,
     Contravariant23Vector,
     Contravariant123Vector
 
-
-
+include("simple_symmetric.jl")
 include("coordinates.jl")
 include("axistensors.jl")
 include("localgeometry.jl")

src/Geometry/axistensors.jl

@@ -136,7 +136,7 @@ struct AxisTensor{
     T,
     N,
     A <: NTuple{N, AbstractAxis},
-    S <: StaticArray{<:Tuple, T, N},
+    S <: Union{SimpleSymmetric{N, T}, StaticArray{<:Tuple, T, N}},
 } <: AbstractArray{T, N}
     axes::A
     components::S
@@ -147,7 +147,7 @@ AxisTensor(
     components::S,
 ) where {
     A <: Tuple{Vararg{AbstractAxis}},
-    S <: StaticArray{<:Tuple, T, N},
+    S <: Union{SimpleSymmetric{N, T}, StaticArray{<:Tuple, T, N}},
 } where {T, N} = AxisTensor{T, N, A, S}(axes, components)
 
 AxisTensor(axes::Tuple{Vararg{AbstractAxis}}, components) =
@@ -396,8 +396,16 @@ end
 function Base.:*(x::AxisVector, y::AdjointAxisVector)
     AxisTensor((axes(x, 1), axes(y, 2)), components(x) * components(y))
 end
+import InteractiveUtils
 function Base.:*(A::Axis2TensorOrAdj, x::AxisVector)
     check_dual(axes(A, 2), axes(x, 1))
+    c = components(A) * components(x)
+    # @show typeof(A)
+    # @show typeof(x)
+    # s = InteractiveUtils.@which components(A) * components(x)
+    # @show s
+    # # @show typeof(x)
+    # @show c
     return AxisVector(axes(A, 1), components(A) * components(x))
 end
 function Base.:*(A::Axis2TensorOrAdj, B::Axis2TensorOrAdj)

src/Geometry/localgeometry.jl

@@ -1,14 +1,24 @@
+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
 
 The necessary local metric information defined at each node.
 """
-struct LocalGeometry{I, C <: AbstractPoint, FT, S}
+struct LocalGeometry{I, C <: AbstractPoint, FT, S, N, L}
     "Coordinates of the current point"
     coordinates::C
     "Jacobian determinant of the transformation `ξ` to `x`"
@@ -22,9 +32,17 @@ 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{N, FT, L},
+    }
     "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{N, FT, L},
+    }
     @inline function LocalGeometry(
         coordinates,
         J,
@@ -34,15 +52,27 @@ 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.")
-        return new{I, C, FT, S}(coordinates, J, WJ, Jinv, ∂x∂ξ, ∂ξ∂x, gⁱʲ, gᵢⱼ)
+        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ᵢⱼ₀)
+        L = triangular_nonzeros(S)
+        N = size(components(gⁱʲ₀), 1)
+        return new{I, C, FT, S, N, L}(
+            coordinates,
+            J,
+            WJ,
+            Jinv,
+            ∂x∂ξ,
+            ∂ξ∂x,
+            gⁱʲ,
+            gᵢⱼ,
+        )
     end
 end
 
-
 """
     SurfaceGeometry
 
@@ -55,7 +85,7 @@ struct SurfaceGeometry{FT, N}
     normal::N
 end
 
-undertype(::Type{LocalGeometry{I, C, FT, S}}) where {I, C, FT, S} = FT
+undertype(::Type{<:LocalGeometry{I, C, FT}}) where {I, C, FT} = FT
 undertype(::Type{SurfaceGeometry{FT, N}}) where {FT, N} = FT
 
 """