--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -110,3 +110,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
+