--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed May 31 08:59:34 2023 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sat Jan 11 10:19:47 2025 +0100
@@ -1,15 +1,15 @@
 using Test
 
-using Sbplib.SbpOperators
-using Sbplib.Grids
-using Sbplib.LazyTensors
+using Diffinitive.SbpOperators
+using Diffinitive.Grids
+using Diffinitive.LazyTensors
 
 @testset "Laplace" begin
     # Default stencils (4th order)
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
     stencil_set = read_stencil_set(operator_path; order=4)
-    g_1D = equidistant_grid(101, 0.0, 1.)
-    g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+    g_1D = equidistant_grid(0.0, 1., 101)
+    g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
 
     @testset "Constructors" begin
         @testset "1D" begin
@@ -69,8 +69,8 @@
 @testset "laplace" begin
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
     stencil_set = read_stencil_set(operator_path; order=4)
-    g_1D = equidistant_grid(101, 0.0, 1.)
-    g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+    g_1D = equidistant_grid(0.0, 1., 101)
+    g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
 
     @testset "1D" begin
         Δ = laplace(g_1D, stencil_set)
@@ -88,3 +88,66 @@
     end
 end
 
+@testset "sat_tensors" begin
+    # TODO: The following tests should be implemented
+    #       1. Symmetry D'H == H'D (test_broken below)
+    #       2. Test eigenvalues of and/or solution to Poisson
+    #       3. Test tuning of Dirichlet conditions
+    #
+    #       These tests are likely easiest to implement once
+    #       we have support for generating matrices from tensors.
+
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    orders = (2,4)
+    tols = (5e-2,5e-4)
+    sz = (201,401)
+    g = equidistant_grid((0.,0.), (1.,1.), sz...)
+    
+    # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO) 
+    # of the operator D = Δ + SAT, where SAT is the tensor composition of the 
+    # operators from sat_tensor. Note that SAT*u should approximate 0 for the 
+    # conditions chosen.
+
+    @testset "Dirichlet" begin
+        for (o, tol) ∈ zip(orders,tols)
+            stencil_set = read_stencil_set(operator_path; order=o)
+            Δ = Laplace(g, stencil_set)
+            H = inner_product(g, stencil_set)
+            u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
+            Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
+            D = Δ 
+            for id ∈ boundary_identifiers(g)
+                D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
+            end
+            e = D*u .- Δu
+            # Accuracy
+            @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+            # Symmetry
+            r = randn(size(u))
+            @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+        end
+    end
+
+    @testset "Neumann" begin
+        @testset "Dirichlet" begin
+            for (o, tol) ∈ zip(orders,tols)
+                stencil_set = read_stencil_set(operator_path; order=o)
+                Δ = Laplace(g, stencil_set)
+                H = inner_product(g, stencil_set)
+                u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
+                Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
+                D = Δ 
+                for id ∈ boundary_identifiers(g)
+                    D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
+                end
+                e = D*u .- Δu
+                # Accuracy
+                @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+                # Symmetry
+                r = randn(size(u))
+                @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+            end
+        end
+    end
+end
+