From f5fdbe9223539c44248032d730caa28878c130b6 Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Fri, 20 Sep 2024 14:46:14 +0100
Subject: [PATCH 1/2] Minor changes in transform and add some comments

---
 src/continuouspolynomial.jl | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/continuouspolynomial.jl b/src/continuouspolynomial.jl
index f53dc60..d1c20f4 100644
--- a/src/continuouspolynomial.jl
+++ b/src/continuouspolynomial.jl
@@ -112,16 +112,20 @@ function *(Pl::ContinuousPolynomialTransform{T,<:Any,<:Any,Int}, X::AbstractMatr
     dims = Pl.dims
     @assert dims == 1
 
-    M,N = size(X,1), size(X,2)
+    M,N = size(X)
     if size(dat,1) ≥ 1
         cfs[Block(1)[1]] = dat[1,1]
+
+        # check that hat functions between elements are consistent
         for j = 1:N-1
-            isapprox(dat[2,j], dat[1,j+1]; atol=1000*M*eps()) || throw(ArgumentError("Discontinuity in data on order of $(abs(dat[2,j]- dat[1,j+1]))."))
+            # TODO: don't check since its fine to return "bad" approximation
+            isapprox(dat[2,j], dat[1,j+1]; atol=1000*M*eps(T)) || throw(ArgumentError("Discontinuity in data on order of $(abs(dat[2,j]- dat[1,j+1]))."))
         end
         for j = 1:N
-            cfs[Block(1)[j+1]] = dat[2,j]
+            cfs[Block(1)[j+1]] = (dat[2,j] + dat[1,j+1])/2 # average to be more robust
         end
     end
+    # populate coefficients of bubble functions
     cfs[Block.(2:M-1)] .= vec(dat[3:end,:]')
     return cfs
 end

From 7cb838dac6a75d4ae96b3556251f9eb922923eaa Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Fri, 11 Oct 2024 16:23:22 +0100
Subject: [PATCH 2/2] MInor fixes

---
 README.md      |  2 +-
 examples/hp.jl | 10 ++++------
 2 files changed, 5 insertions(+), 7 deletions(-)

diff --git a/README.md b/README.md
index a998304..749d015 100644
--- a/README.md
+++ b/README.md
@@ -16,7 +16,7 @@ plot(C[:,Block.(2:3)])
 
 The mass matrix can be constructed via:
 ```julia
-M = C'C
+M = grammatrix(C)
 ```
 We can also construct the stiffness matrix:
 ```julia
diff --git a/examples/hp.jl b/examples/hp.jl
index 00a6a04..12ef851 100644
--- a/examples/hp.jl
+++ b/examples/hp.jl
@@ -1,4 +1,5 @@
-using PiecewiseOrthogonalPolynomials, Plots
+using PiecewiseOrthogonalPolynomials
+using Plots
 
 ## Goal 1: Do hp solves in 1D where complexity is optimal 
 # regardless of h or p
@@ -7,15 +8,12 @@ r = range(0, 1; length=4)
 
 # C is standard affine FEM combined with mapped (1-x^2) * P_k^(1,1)(x)
 C = ContinuousPolynomial{1}(r)
-
-g = range(0,1; length=1000)
-plot(g, C[g,Block(4)])
+plot(C[:,Block(4)])
 
 x = axes(C,1)
-D = Derivative(x)
 
 M = C'C # Mass matrix
-Δ = -(D*C)'*(D*C) # Weak Laplacian
+Δ = -diff(D)'*diff(C) # Weak Laplacian
 
 
 # For 2D/3D use Fortanto & Townsend have fast spectral Poisson/Helmholtz that