From 0db40b7832a19d0b93a000dd0559e8b638893e69 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sat, 8 Aug 2020 23:17:44 +0200
Subject: [PATCH 01/10] Add distinct type for grid vs. function value

---
 src/derivative_operators/BC_operators.jl | 32 ++++++++++++------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index 0f05711e0..ecde7faac 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -22,11 +22,11 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
     b_l::T
     a_r::V
     b_r::T
-    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
+    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::U, order = 1) where {T,U}
         αl, βl, γl = l
         αr, βr, γr = r
 
-        s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order
+        s = calculate_weights(1, one(U), Array(one(U):convert(U,order+1))) #generate derivative coefficients about the boundary of required approximation order
 
         a_l = -s[2:end]./(αl*dx/βl + s[1])
         a_r = s[end:-1:2]./(αr*dx/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*
@@ -36,14 +36,14 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
 
         return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
     end
-    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{T}, order = 1) where {T}
+    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{U}, order = 1) where {T,U}
         αl, βl, γl = l
         αr, βr, γr = r
 
-        s_index = Array(one(T):convert(T,order+1))
+        s_index = Array(one(U):convert(U,order+1))
         dx_l, dx_r = dx[1:length(s_index)], dx[(end-length(s_index)+1):end]
 
-        s = calculate_weights(1, one(T), s_index) #generate derivative coefficients about the boundary of required approximation order
+        s = calculate_weights(1, one(U), s_index) #generate derivative coefficients about the boundary of required approximation order
         denom_l = αl+βl*s[1]/dx_l[1]
         denom_r = αr-βr*s[1]/dx_r[end]
 
@@ -76,24 +76,24 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
     b_l::T
     a_r::R
     b_r::T
-    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::U, order = 1) where {T,U}
         nl = length(αl)
         nr = length(αr)
         S_l = zeros(T, (nl-2, order+nl-2))
         S_r = zeros(T, (nr-2, order+nr-2))
 
         for i in 1:(nl-2)
-            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
+            S_l[i,:] = [transpose(calculate_weights(i, one(U), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
         end
 
         for i in 1:(nr-2)
-            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx^i)
+            S_r[i,:] = [transpose(calculate_weights(i, convert(U, order+i), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nr-2-i)))]./(dx^i)
         end
         s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
         s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
-        denomr = αr[2] .+ αr[3:end] ⋅ s0_r
+        denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
+        denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
         a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -103,7 +103,7 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
     end
 
-    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{U}, order = 1) where {T,U}
 
         nl = length(αl)
         nr = length(αr)
@@ -112,17 +112,17 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         S_r = zeros(T, (nr-2, order+nr-2))
 
         for i in 1:(nl-2)
-            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx_l.^i)
+            S_l[i,:] = [transpose(calculate_weights(i, one(U), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nl-2-i)))]./(dx_l.^i)
         end
 
         for i in 1:(nr-2)
-            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx_r.^i)
+            S_r[i,:] = [transpose(calculate_weights(i, convert(U, order+i), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nr-2-i)))]./(dx_r.^i)
         end
         s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
         s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
-        denomr = αr[2] .+ αr[3:end] ⋅ s0_r
+        denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
+        denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
         a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -145,7 +145,7 @@ Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero
 # other acceptable argument signatures
 #RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)
 
-Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
+Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l'⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r' ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
 Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)
 
 Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?

From b6b39f6b51f340d10004e74196fcce3edf04c480 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sat, 8 Aug 2020 23:19:10 +0200
Subject: [PATCH 02/10] Test Robin BC with complex function values

---
 test/robin_complex.jl | 46 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)
 create mode 100644 test/robin_complex.jl

diff --git a/test/robin_complex.jl b/test/robin_complex.jl
new file mode 100644
index 000000000..fab2ab427
--- /dev/null
+++ b/test/robin_complex.jl
@@ -0,0 +1,46 @@
+include("src/DiffEqOperators.jl")
+using .DiffEqOperators
+using LinearAlgebra, Random, Test
+
+# Generate random parameters
+# Generate random parameters
+al = rand(ComplexF64,5)
+bl = rand(ComplexF64,5)
+cl = rand(ComplexF64,5)
+dx = rand(Float64,5)
+ar = rand(ComplexF64,5)
+br = rand(ComplexF64,5)
+cr = rand(ComplexF64,5)
+
+# Construct 5 arbitrary RobinBC operators for each parameter set
+for i in 1:5
+	
+	Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])
+
+	Q_L, Q_b = Array(Q,5i)
+
+	#Check that Q_L is is correctly computed
+	@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
+	@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
+	@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
+
+	#Check that Q_b is computed correctly
+	@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
+
+	# Construct the extended operator and check that it correctly extends u to a (5i+2)
+	# vector, along with encoding boundary condition information.
+	u = rand(ComplexF64,5i)
+
+	Qextended = Q*u
+	CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
+	@test length(Qextended) ≈ 5i+2
+
+	# Check concretization
+	@test Array(Qextended) ≈ CorrectQextended # 	Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, 		Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r
+
+	# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
+	@test Q_L*u + Q_b ≈ CorrectQextended
+
+	@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
+	
+end

From 0d3d564d0e9de4d2db011a86c493386a95174140 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sat, 8 Aug 2020 23:24:22 +0200
Subject: [PATCH 03/10] duplicate comment

---
 test/robin_complex.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/robin_complex.jl b/test/robin_complex.jl
index fab2ab427..e0a737f04 100644
--- a/test/robin_complex.jl
+++ b/test/robin_complex.jl
@@ -2,7 +2,6 @@ include("src/DiffEqOperators.jl")
 using .DiffEqOperators
 using LinearAlgebra, Random, Test
 
-# Generate random parameters
 # Generate random parameters
 al = rand(ComplexF64,5)
 bl = rand(ComplexF64,5)

From 0f635587626a11e6855e641aa83cc53cf2600478 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sun, 9 Aug 2020 00:19:21 +0200
Subject: [PATCH 04/10] Add new test to runtests

---
 test/runtests.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index d0ffa4dfb..8da9cafb0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,6 +12,7 @@ if GROUP == "All" || GROUP == "Interface"
     @time @safetestset "Poisson example" begin include("../examples/poisson.jl") end
     @time @safetestset "Heat equation example" begin include("../examples/heat_equation.jl") end
     @time @safetestset "Robin Boundary Condition Operators" begin include("robin.jl") end
+    @time @safetestset "Robin Boundary Condition Operators, Complex" begin include("robin_complex.jl") end
     @time @safetestset "JacVec Operators Interface" begin include("jacvec_operators.jl") end
     @time @safetestset "Composite Operators Interface" begin include("composite_operators_interface.jl") end
     @time @safetestset "BC and Coefficient Compositions" begin include("bc_coeff_compositions.jl") end

From e8a2d9625fe7b8f75481c7d67de2e63842ee76a7 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sun, 9 Aug 2020 00:19:43 +0200
Subject: [PATCH 05/10] Restrict types to avoid ambiguity

---
 src/derivative_operators/BC_operators.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index ecde7faac..b7ca9888c 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -22,7 +22,7 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
     b_l::T
     a_r::V
     b_r::T
-    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::U, order = 1) where {T,U}
+    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::U, order = 1) where {T<:Number,U<:Real}
         αl, βl, γl = l
         αr, βr, γr = r
 
@@ -36,7 +36,7 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
 
         return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
     end
-    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{U}, order = 1) where {T,U}
+    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{U}, order = 1) where {T<:Number,U<:Real}
         αl, βl, γl = l
         αr, βr, γr = r
 
@@ -76,7 +76,7 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
     b_l::T
     a_r::R
     b_r::T
-    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::U, order = 1) where {T,U}
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::U, order = 1) where {T<:Number,U<:Real}
         nl = length(αl)
         nr = length(αr)
         S_l = zeros(T, (nl-2, order+nl-2))
@@ -103,7 +103,7 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
     end
 
-    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{U}, order = 1) where {T,U}
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{U}, order = 1) where {T<:Number,U<:Real}
 
         nl = length(αl)
         nr = length(αr)

From b558da84cf94443a6dfbee7991097677e7d361e5 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sun, 9 Aug 2020 01:29:18 +0200
Subject: [PATCH 06/10] More tests, type restrictions

---
 src/derivative_operators/BC_operators.jl | 11 ++-
 test/robin.jl                            | 95 ++++++++++++++++++++++++
 test/runtests.jl                         |  1 -
 3 files changed, 102 insertions(+), 5 deletions(-)

diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index b7ca9888c..64fa0921d 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -136,16 +136,19 @@ end
 
 
 #implement Neumann and Dirichlet as special cases of RobinBC
-NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
-DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
+NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{U}, U}, order = 1) where {T<:Number,U<:Real} = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
+DirichletBC(αl::T, αr::T) where {T<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
+DirichletBC(U::Type,αl::T, αr::T) where {T<:Number,U<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(U), 2one(U) )
 #specialized constructors for Neumann0 and Dirichlet0
 Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
-Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)
+Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real} = NeumannBC((zero(T), zero(T)), dx, order)
+Neumann0BC(U::Type,dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,U<:Number} = NeumannBC((zero(U), zero(U)), dx, order)
 
 # other acceptable argument signatures
 #RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)
 
-Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l'⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r' ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
+Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector( Q.a_l'⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r' ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u )
+
 Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)
 
 Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?
diff --git a/test/robin.jl b/test/robin.jl
index d7c8b3882..5237122ee 100644
--- a/test/robin.jl
+++ b/test/robin.jl
@@ -104,3 +104,98 @@ ugeneralextended = G*u
 
 
 #TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary
+
+
+
+
+# Test complex Robin BC, uniform grid
+
+# Generate random parameters
+al = rand(ComplexF64,5)
+bl = rand(ComplexF64,5)
+cl = rand(ComplexF64,5)
+dx = rand(Float64,5)
+ar = rand(ComplexF64,5)
+br = rand(ComplexF64,5)
+cr = rand(ComplexF64,5)
+
+# Construct 5 arbitrary RobinBC operators for each parameter set
+for i in 1:5
+	
+	Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])
+
+	Q_L, Q_b = Array(Q,5i)
+
+	#Check that Q_L is is correctly computed
+	@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
+	@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
+	@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
+
+	#Check that Q_b is computed correctly
+	@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
+
+	# Construct the extended operator and check that it correctly extends u to a (5i+2)
+	# vector, along with encoding boundary condition information.
+	u = rand(ComplexF64,5i)
+
+	Qextended = Q*u
+	CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
+	@test length(Qextended) ≈ 5i+2
+
+	# Check concretization
+	@test Array(Qextended) ≈ CorrectQextended # 	Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, 		Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r
+
+	# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
+	@test Q_L*u + Q_b ≈ CorrectQextended
+
+	@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
+	
+end
+
+# Construct 5 arbitrary RobinBC operators w/non-uniform grid
+al = rand(ComplexF64,5)
+bl = rand(ComplexF64,5)
+cl = rand(ComplexF64,5)
+dx = rand(Float64,5)
+ar = rand(ComplexF64,5)
+br = rand(ComplexF64,5)
+cr = rand(ComplexF64,5)
+for j in 1:2
+    for i in 1:5
+        if j == 1
+            Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]),
+                        dx[i] .* ones(5 * i))
+        else
+            Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]],
+                        dx[i] .* ones(5 * i))
+        end
+        Q_L, Q_b = Array(Q,5i)
+
+        #Check that Q_L is is correctly computed
+        @test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
+        @test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
+        @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
+
+        #Check that Q_b is computed correctly
+        @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
+
+        # Construct the extended operator and check that it correctly extends u to a (5i+2)
+        # vector, along with encoding boundary condition information.
+        u = rand(ComplexF64,5i)
+
+        Qextended = Q*u
+        CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
+        @test length(Qextended) ≈ 5i+2
+
+        # Check concretization
+        @test Array(Qextended) ≈ CorrectQextended
+
+        # Check that Q_L and Q_b correctly compute BoundaryPaddedVector
+        @test Q_L*u + Q_b ≈ CorrectQextended
+
+        @test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
+
+    end
+end
+
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 8da9cafb0..d0ffa4dfb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,7 +12,6 @@ if GROUP == "All" || GROUP == "Interface"
     @time @safetestset "Poisson example" begin include("../examples/poisson.jl") end
     @time @safetestset "Heat equation example" begin include("../examples/heat_equation.jl") end
     @time @safetestset "Robin Boundary Condition Operators" begin include("robin.jl") end
-    @time @safetestset "Robin Boundary Condition Operators, Complex" begin include("robin_complex.jl") end
     @time @safetestset "JacVec Operators Interface" begin include("jacvec_operators.jl") end
     @time @safetestset "Composite Operators Interface" begin include("composite_operators_interface.jl") end
     @time @safetestset "BC and Coefficient Compositions" begin include("bc_coeff_compositions.jl") end

From 61066328f721238fea01a3c315a4757ffcb6f62a Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sun, 9 Aug 2020 01:40:44 +0200
Subject: [PATCH 07/10] Fix Dirichlet / Neumann

---
 src/derivative_operators/BC_operators.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index 64fa0921d..a9f8dee88 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -138,11 +138,11 @@ end
 #implement Neumann and Dirichlet as special cases of RobinBC
 NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{U}, U}, order = 1) where {T<:Number,U<:Real} = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
 DirichletBC(αl::T, αr::T) where {T<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
-DirichletBC(U::Type,αl::T, αr::T) where {T<:Number,U<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(U), 2one(U) )
+DirichletBC(::Type{U},αl::T, αr::T) where {T<:Number,U<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(U), 2one(U) )
 #specialized constructors for Neumann0 and Dirichlet0
 Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
 Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real} = NeumannBC((zero(T), zero(T)), dx, order)
-Neumann0BC(U::Type,dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,U<:Number} = NeumannBC((zero(U), zero(U)), dx, order)
+Neumann0BC(::Type{U},dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,U<:Number} = NeumannBC((zero(U), zero(U)), dx, order)
 
 # other acceptable argument signatures
 #RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)

From d477f23d4417f35cbe9d8814491316ee64d973fb Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sun, 9 Aug 2020 12:19:46 +0200
Subject: [PATCH 08/10] Add tests for Neuman/Dirichlet w. Complex numbers

---
 test/robin.jl | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/test/robin.jl b/test/robin.jl
index 5237122ee..5711e3643 100644
--- a/test/robin.jl
+++ b/test/robin.jl
@@ -198,4 +198,19 @@ for j in 1:2
     end
 end
 
+# Test Neumann and Dirichlet as special cases of RobinBC
+dx = [0.121, 0.783, 0.317, 0.518, 0.178]
+αC = (0.539 + 0.653im, 0.842 + 0.47im)
+αR = (0.045, 0.577)
+@test NeumannBC(αC, dx).b_l ≈ -0.065219 - 0.079013im
+@test DirichletBC(αR...).b_r ≈ 0.577
+@test DirichletBC(Float64, αC...) ≈ 0.123 # broken
+
+@test Dirichlet0BC(Float64).a_r ≈ [-0.0,0.0]
+@test Neumann0BC(dx).a_r ≈ [0.3436293436293436]
+@test Neumann0BC(ComplexF64,dx).a_l ≈ [0.15453384418901658 + 0.0im]
+
+@test NeumannBC(αC, first(dx)).b_r ≈ 0.101882 + 0.05687im
+@test Neumann0BC(first(dx)).a_r ≈ [1.0 - 0.0im]
+@test Neumann0BC(ComplexF64,first(dx)).a_l ≈ [1.0 + 0.0im]
 

From 83a9b6c23b125271267dabff227c40a7254579fc Mon Sep 17 00:00:00 2001
From: AndiMD <original-andy@gmx.de>
Date: Sat, 8 Aug 2020 23:17:44 +0200
Subject: [PATCH 09/10] Add distinct type for grid vs. function value Test
 Robin BC, Dirichlet, Neumann with complex function values

---
 src/derivative_operators/BC_operators.jl |  45 ++++-----
 test/robin.jl                            | 112 +++++++++++++++++++++++
 2 files changed, 136 insertions(+), 21 deletions(-)

diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index 0f05711e0..984713e94 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -22,28 +22,28 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
     b_l::T
     a_r::V
     b_r::T
-    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
+    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::U, order = 1) where {T<:Number,U<:Real}
         αl, βl, γl = l
         αr, βr, γr = r
 
-        s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order
+        s = calculate_weights(1, one(U), Array(one(U):convert(U,order+1))) #generate derivative coefficients about the boundary of required approximation order
 
-        a_l = -s[2:end]./(αl*dx/βl + s[1])
-        a_r = s[end:-1:2]./(αr*dx/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*
+        a_l = -βl*s[2:end]./(αl*dx + βl*s[1])
+        a_r = βr*s[end:-1:2]./(αr*dx - βr*s[1]) # for other boundary stencil is flippedlr with *opposite sign*
 
         b_l = γl/(αl+βl*s[1]/dx)
         b_r = γr/(αr-βr*s[1]/dx)
 
         return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
     end
-    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{T}, order = 1) where {T}
+    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{U}, order = 1) where {T<:Number,U<:Real}
         αl, βl, γl = l
         αr, βr, γr = r
 
-        s_index = Array(one(T):convert(T,order+1))
+        s_index = Array(one(U):convert(U,order+1))
         dx_l, dx_r = dx[1:length(s_index)], dx[(end-length(s_index)+1):end]
 
-        s = calculate_weights(1, one(T), s_index) #generate derivative coefficients about the boundary of required approximation order
+        s = calculate_weights(1, one(U), s_index) #generate derivative coefficients about the boundary of required approximation order
         denom_l = αl+βl*s[1]/dx_l[1]
         denom_r = αr-βr*s[1]/dx_r[end]
 
@@ -76,24 +76,24 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
     b_l::T
     a_r::R
     b_r::T
-    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::U, order = 1) where {T<:Number,U<:Real}
         nl = length(αl)
         nr = length(αr)
         S_l = zeros(T, (nl-2, order+nl-2))
         S_r = zeros(T, (nr-2, order+nr-2))
 
         for i in 1:(nl-2)
-            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
+            S_l[i,:] = [transpose(calculate_weights(i, one(U), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
         end
 
         for i in 1:(nr-2)
-            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx^i)
+            S_r[i,:] = [transpose(calculate_weights(i, convert(U, order+i), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nr-2-i)))]./(dx^i)
         end
         s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
         s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
-        denomr = αr[2] .+ αr[3:end] ⋅ s0_r
+        denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
+        denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
         a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -103,7 +103,7 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
     end
 
-    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
+    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{U}, order = 1) where {T<:Number,U<:Real}
 
         nl = length(αl)
         nr = length(αr)
@@ -112,17 +112,17 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         S_r = zeros(T, (nr-2, order+nr-2))
 
         for i in 1:(nl-2)
-            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx_l.^i)
+            S_l[i,:] = [transpose(calculate_weights(i, one(U), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nl-2-i)))]./(dx_l.^i)
         end
 
         for i in 1:(nr-2)
-            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx_r.^i)
+            S_r[i,:] = [transpose(calculate_weights(i, convert(U, order+i), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nr-2-i)))]./(dx_r.^i)
         end
         s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
         s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
-        denomr = αr[2] .+ αr[3:end] ⋅ s0_r
+        denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
+        denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
         a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -136,16 +136,19 @@ end
 
 
 #implement Neumann and Dirichlet as special cases of RobinBC
-NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
-DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
+NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{U}, U}, order = 1) where {T<:Number,U<:Real} = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
+DirichletBC(αl::T, αr::T) where {T<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
+DirichletBC(::Type{U},αl::T, αr::T) where {T<:Number,U<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(U), 2one(U) )
 #specialized constructors for Neumann0 and Dirichlet0
 Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
-Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)
+Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real} = NeumannBC((zero(T), zero(T)), dx, order)
+Neumann0BC(::Type{U},dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,U<:Number} = NeumannBC((zero(U), zero(U)), dx, order)
 
 # other acceptable argument signatures
 #RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)
 
-Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
+Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector( Q.a_l'⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r' ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u )
+
 Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)
 
 Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?
diff --git a/test/robin.jl b/test/robin.jl
index d7c8b3882..aaf115036 100644
--- a/test/robin.jl
+++ b/test/robin.jl
@@ -104,3 +104,115 @@ ugeneralextended = G*u
 
 
 #TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary
+
+
+
+
+# Test complex Robin BC, uniform grid
+
+# Generate random parameters
+al = rand(ComplexF64,5)
+bl = rand(ComplexF64,5)
+cl = rand(ComplexF64,5)
+dx = rand(Float64,5)
+ar = rand(ComplexF64,5)
+br = rand(ComplexF64,5)
+cr = rand(ComplexF64,5)
+
+# Construct 5 arbitrary RobinBC operators for each parameter set
+for i in 1:5
+	
+	Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])
+
+	Q_L, Q_b = Array(Q,5i)
+
+	#Check that Q_L is is correctly computed
+	@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
+	@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
+	@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
+
+	#Check that Q_b is computed correctly
+	@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
+
+	# Construct the extended operator and check that it correctly extends u to a (5i+2)
+	# vector, along with encoding boundary condition information.
+	u = rand(ComplexF64,5i)
+
+	Qextended = Q*u
+	CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
+	@test length(Qextended) ≈ 5i+2
+
+	# Check concretization
+	@test Array(Qextended) ≈ CorrectQextended # 	Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, 		Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r
+
+	# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
+	@test Q_L*u + Q_b ≈ CorrectQextended
+
+	@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
+	
+end
+
+# Test complex Robin BC, w/non-uniform grid
+al = rand(ComplexF64,5)
+bl = rand(ComplexF64,5)
+cl = rand(ComplexF64,5)
+dx = rand(Float64,5)
+ar = rand(ComplexF64,5)
+br = rand(ComplexF64,5)
+cr = rand(ComplexF64,5)
+for j in 1:2
+    for i in 1:5
+        if j == 1
+            Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]),
+                        dx[i] .* ones(5 * i))
+        else
+            Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]],
+                        dx[i] .* ones(5 * i))
+        end
+        Q_L, Q_b = Array(Q,5i)
+
+        #Check that Q_L is is correctly computed
+        @test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
+        @test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
+        @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
+
+        #Check that Q_b is computed correctly
+        @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
+
+        # Construct the extended operator and check that it correctly extends u to a (5i+2)
+        # vector, along with encoding boundary condition information.
+        u = rand(ComplexF64,5i)
+
+        Qextended = Q*u
+        CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
+        @test length(Qextended) ≈ 5i+2
+
+        # Check concretization
+        @test Array(Qextended) ≈ CorrectQextended
+
+        # Check that Q_L and Q_b correctly compute BoundaryPaddedVector
+        @test Q_L*u + Q_b ≈ CorrectQextended
+
+        @test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
+
+    end
+end
+
+# Test NeumannBC, DirichletBC as special cases of RobinBC
+let
+	dx = [0.121, 0.783, 0.317, 0.518, 0.178]
+	αC = (0.539 + 0.653im, 0.842 + 0.47im)
+	αR = (0.045, 0.577)
+	@test NeumannBC(αC, dx).b_l ≈ -0.065219 - 0.079013im
+	@test DirichletBC(αR...).b_r ≈ 0.577
+	@test DirichletBC(Float64, αC...).b_l ≈ 0.539 + 0.653im
+	@test DirichletBC(Float64, αC...).a_r ≈ [-0.0 + 0.0im, 0.0 + 0.0im]
+
+	@test Dirichlet0BC(Float64).a_r ≈ [-0.0,0.0]
+	@test Neumann0BC(dx).a_r ≈ [0.3436293436293436]
+	@test Neumann0BC(ComplexF64,dx).a_l ≈ [0.15453384418901658 + 0.0im]
+
+	@test NeumannBC(αC, first(dx)).b_r ≈ 0.101882 + 0.05687im
+	@test Neumann0BC(first(dx)).a_r ≈ [1.0 - 0.0im]
+	@test Neumann0BC(ComplexF64,first(dx)).a_l ≈ [1.0 + 0.0im]
+end

From 197c9f284d243580e76d6445162b3a12b5886f65 Mon Sep 17 00:00:00 2001
From: AndiMD <none@gmx.net>
Date: Sat, 5 Sep 2020 13:11:41 +0200
Subject: [PATCH 10/10] Express unconjugated dot product by sum(a.*b)

---
 src/derivative_operators/BC_operators.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl
index 984713e94..debcb6f74 100644
--- a/src/derivative_operators/BC_operators.jl
+++ b/src/derivative_operators/BC_operators.jl
@@ -92,8 +92,8 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
         s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-        denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
-        denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
+        denoml = αl[2] .+ sum(αl[3:end] .* s0_l) # dot product without complex conjugation
+        denomr = αr[2] .+ sum(αr[3:end] .* s0_r)
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
         a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -121,8 +121,8 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
         s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
         s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
 
-        denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
-        denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
+        denoml = αl[2] .+ sum(αl[3:end] .* s0_l)
+        denomr = αr[2] .+ sum(αr[3:end] .* s0_r)
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
         a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -147,7 +147,7 @@ Neumann0BC(::Type{U},dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,
 # other acceptable argument signatures
 #RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)
 
-Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector( Q.a_l'⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r' ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u )
+Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector( sum(Q.a_l .* u[1:length(Q.a_l)]) + Q.b_l, sum(Q.a_r .* u[(end-length(Q.a_r)+1):end]) + Q.b_r, u )
 
 Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)