Use QR decomposition in DD sharp

src/DiscreteExteriorCalculus.jl

@@ -28,9 +28,9 @@ export DualSimplex, DualV, DualE, DualTri, DualChain, DualForm,
 import Base: ndims
 import Base: *
 import LinearAlgebra: mul!
-using LinearAlgebra: Diagonal, dot, norm, cross, pinv
+using LinearAlgebra: Diagonal, dot, norm, cross, pinv, qr, ColumnNorm
 using SparseArrays
-using StaticArrays: @SVector, SVector, SMatrix
+using StaticArrays: @SVector, SVector, SMatrix, MMatrix
 using GeometryBasics: Point2, Point3
 
 const Point2D = SVector{2,Float64}
@@ -801,7 +801,11 @@ function ♯_mat(s::AbstractDeltaDualComplex2D, ::LLSDDSharp)
       end
     # TODO: Move around ' as appropriate to minimize transposing.
     X = stack(de_vecs)'
-    LLS = pinv(X'*(X))*(X')
+    # See: https://arxiv.org/abs/1102.1845
+    #QRX = qr(X, ColumnNorm())
+    #LLS = (inv(QRX.R) * QRX.Q')[QRX.p,:]
+    #LLS = pinv(X'*(X))*(X')
+    LLS = pinv(X)
     for (i,e) in enumerate(tri_edges)
       ♯_m[t, e] = LLS[:,i]'*weights[i]
     end

test/DiscreteExteriorCalculus.jl

@@ -586,7 +586,7 @@ for (primal_s,s) in flat_meshes
   # This tests that numerically raising indices is equivalent to analytically
   # raising indices.
   # X_continuous = ♯ᵖ_discrete∘♭ᵖ_discrete(X_continuous)
-  X♯ = SVector{3,Float64}(1/√2,1/√2,0)
+  X♯ = SVector{3,Float64}(3/√2,1/√2,0)
   X♭ = eval_constant_dual_form(s, X♯)
   @test all_are_approx(♯_m * X♭,  X♯; atol=1e-13) ||
         all_are_approx(♯_m * X♭, -X♯; atol=1e-13)
@@ -599,8 +599,8 @@ for (primal_s,s) in flat_meshes
   # Note that we check explicitly both cases of signedness, because orient!
   # picks in/outside without regard to the embedding.
   # This is the "right/left-hand-rule trick."
-  @test all(isapprox.(only.(♭_♯_m * X♭),  eval_constant_primal_form(s, X♯), atol=1e-12)) ||
-        all(isapprox.(only.(♭_♯_m * X♭), -eval_constant_primal_form(s, X♯), atol=1e-12))
+  @test all(isapprox.(only.(♭_♯_m * X♭),  eval_constant_primal_form(s, X♯), atol=1e-13)) ||
+        all(isapprox.(only.(♭_♯_m * X♭), -eval_constant_primal_form(s, X♯), atol=1e-13))
 
   # This test shows how the musical isomorphism chaining lets you further
   # define a primal-dual wedge that preserves properties from the continuous
@@ -624,8 +624,8 @@ for (primal_s,s) in flat_meshes
   @test all(Λdp(f̃, g′) .== -1 * Λpd(g′, f̃))
 
   # f∧f = 0 (implied by antisymmetry):
-  @test all(isapprox.(Λpd(f′, f̃), 0, atol=1e-10))
-  @test all(isapprox.(Λpd(g′, g̃), 0, atol=1e-10))
+  @test all(isapprox.(Λpd(f′, f̃), 0, atol=1e-11))
+  @test all(isapprox.(Λpd(g′, g̃), 0, atol=1e-11))
 
   # Test and demonstrate the convenience functions:
   @test all(∧(s, SimplexForm{1}(f′), DualForm{1}(g̃)) .≈ -1 * ∧(s, SimplexForm{1}(g′), DualForm{1}(f̃)))
@@ -654,7 +654,7 @@ for (primal_s,s) in flat_meshes
   @test all(isapprox.(
     dec_hodge_star(2,s) * ∧(s, SimplexForm{1}(u), DualForm{1}(u_star)),
     ff_gg,
-    atol=1e-10))
+    atol=1e-12))
 end
 
 end

test/Operators.jl

@@ -278,7 +278,7 @@ end
       # f∧a = (4dx - dy) ∧ (x + 4y)
       #     = 4(x + 4y)dx -(x + 4y)dy
       #     = (4x + 16y)dx + (-x - 4y)dy
-      @test all(isapprox.(Λ10(f,a)[interior_edges], fΛa_analytic[interior_edges], atol=1e-8))
+      @test all(isapprox.(Λ10(f,a)[interior_edges], fΛa_analytic[interior_edges], atol=1e-10))
     end
 end
 
@@ -393,7 +393,7 @@ end
   # This smaller mesh is proportionally more affected by boundary conditions.
   X♯ = SVector{3,Float64}(1/√2,1/√2,0)
   mag_selfadv, mag_dp, mag_∂ₜu  = euler_equation_test(X♯, tg)
-  @test .60 < (count(mag_selfadv .< 1e-2) / length(mag_selfadv))
+  @test .64 < (count(mag_selfadv .< 1e-2) / length(mag_selfadv))
   @test .64 < (count(mag_dp .< 1e-2) / length(mag_dp))
   @test .60 < (count(mag_∂ₜu .< 1e-2) / length(mag_∂ₜu))