From 1c033dfdef73fd762544e81eadf15cf5fcf96f8b Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Fri, 12 May 2023 17:47:32 -0500
Subject: [PATCH] fix `approximation_type::Polynomial{Gauss}` constructor
 (#123)

* fix type in constructor

* fix `sparse_low_order_SBP_1D_operators`

* fix low order sparse operators

* add Gauss test

* reduce number of tests
---
 src/RefElemData_SBP.jl        | 26 +++++++++++++++++++-------
 src/RefElemData_polynomial.jl |  5 ++---
 test/reference_elem_tests.jl  | 13 +++++++++++--
 3 files changed, 32 insertions(+), 12 deletions(-)

diff --git a/src/RefElemData_SBP.jl b/src/RefElemData_SBP.jl
index e16da3cb..20c287a1 100644
--- a/src/RefElemData_SBP.jl
+++ b/src/RefElemData_SBP.jl
@@ -312,7 +312,7 @@ function sparse_low_order_SBP_operators(rd::RefElemData{NDIMS}) where {NDIMS}
     return sparse.(Qrst), sparse(E)
 end
 
-function sparse_low_order_SBP_1D_operators(rd::RefElemData{1}) 
+function sparse_low_order_SBP_1D_operators(rd::RefElemData) 
     E = zeros(2, rd.N+1)
     E[1, 1] = 1
     E[2, end] = 1
@@ -329,18 +329,29 @@ end
 sparse_low_order_SBP_operators(rd::RefElemData{1, Line, <:Union{<:SBP, <:Polynomial{Gauss}}}) = 
     sparse_low_order_SBP_1D_operators(rd)
 
+function diagonal_1D_mass_matrix(N, ::SBP)
+    _, w1D = gauss_lobatto_quad(0, 0, N)
+    return Diagonal(w1D)
+end
+
+function diagonal_1D_mass_matrix(N, ::Polynomial{Gauss})
+    _, w1D = gauss_quad(0, 0, N)
+    return Diagonal(w1D)
+end
+
 function sparse_low_order_SBP_operators(rd::RefElemData{2, Quad, <:Union{<:SBP, <:Polynomial{Gauss}}}) 
     (Q1D,), E1D = sparse_low_order_SBP_1D_operators(rd)
 
     # permute face node ordering for the first 2 faces
-    ids = reshape(1:(rd.N+1) * 2, :, 2)    
+    ids = reshape(1:(rd.N+1) * 2, :, 2)
     Er = zeros((rd.N+1) * 2, rd.Np)
     Er[vec(ids'), :] .= kron(I(rd.N+1), E1D)
     Es = kron(E1D, I(rd.N+1))
     E = vcat(Er, Es)
 
-    Qr = kron(I(rd.N+1), Q1D)
-    Qs = kron(Q1D, I(rd.N+1))
+    M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type)
+    Qr = kron(M1D, Q1D)
+    Qs = kron(Q1D, M1D)
 
     return sparse.((Qr, Qs)), sparse(E)
 end
@@ -358,9 +369,10 @@ function sparse_low_order_SBP_operators(rd::RefElemData{3, Hex, <:Union{<:SBP, <
     # create boundary extraction matrix
     E = vcat(Er, Es, Et)
 
-    Qr = kron(I(rd.N+1), Q1D, I(rd.N+1))
-    Qs = kron(I(rd.N+1), I(rd.N+1), Q1D)
-    Qt = kron(Q1D, I(rd.N+1), I(rd.N+1))
+    M1D = diagonal_1D_mass_matrix(rd.N, rd.approximation_type)
+    Qr = kron(M1D, Q1D, M1D)
+    Qs = kron(M1D, M1D, Q1D)
+    Qt = kron(Q1D, M1D, M1D)
 
     return sparse.((Qr, Qs, Qt)), sparse(E)
 end
diff --git a/src/RefElemData_polynomial.jl b/src/RefElemData_polynomial.jl
index 9221c078..16cc9fed 100644
--- a/src/RefElemData_polynomial.jl
+++ b/src/RefElemData_polynomial.jl
@@ -356,9 +356,8 @@ function RefElemData(element_type::Union{Line, Quad, Hex},
                      approximation_type::Polynomial{<:Gauss}, N; kwargs...) 
     # explicitly specify Gauss quadrature rule with (N+1)^d points 
     quad_rule_vol = tensor_product_quadrature(element_type, StartUpDG.gauss_quad(0, 0, N)...)
-    rd = RefElemData(element_type, Polynomial(), N; quad_rule_vol)
-    @set rd.approximation_type = approximation_type
-    return rd
+    rd = RefElemData(element_type, Polynomial(), N; quad_rule_vol)    
+    return @set rd.approximation_type = approximation_type
 end
   
 RefElemData(element_type::Union{Line, Quad, Hex}, ::Gauss, N; kwargs...) = 
diff --git a/test/reference_elem_tests.jl b/test/reference_elem_tests.jl
index c171eeac..062d496f 100644
--- a/test/reference_elem_tests.jl
+++ b/test/reference_elem_tests.jl
@@ -187,8 +187,8 @@ inverse_trace_constant_compare(rd::RefElemData{3, <:Wedge, <:TensorProductWedge}
 
 @testset "Tensor product wedges" begin
     tol = 5e2*eps()
-    @testset "Degree $tri_grad triangle" for tri_grad = [1, 2, 3, 4]
-        @testset "Degree $line_grad line" for line_grad = [1, 2, 3, 4]
+    @testset "Degree $tri_grad triangle" for tri_grad = [2, 3]
+        @testset "Degree $line_grad line" for line_grad = [1, 2]
         line = RefElemData(Line(), line_grad)
         tri  = RefElemData(Tri(), tri_grad)
         tensor = TensorProductWedge(tri, line)
@@ -292,6 +292,15 @@ end
         @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E
     end
 
+    @testset "Quad (SBP)" begin
+        rd = RefElemData(Quad(), SBP(), N)
+        (Qr, Qs), E = sparse_low_order_SBP_operators(rd)
+        @test norm(sum(Qr, dims=2)) < tol
+        @test norm(sum(Qs, dims=2)) < tol
+        @test Qr + Qr' ≈ E' * Diagonal(rd.wf .* rd.nrJ) * E
+        @test Qs + Qs' ≈ E' * Diagonal(rd.wf .* rd.nsJ) * E
+    end
+
     @testset "Quad (Polynomial{Gauss})" begin
         rd = RefElemData(Quad(), Gauss(), N)
         (Qr, Qs), E = sparse_low_order_SBP_operators(rd)