Add distinct type for grid vs. function value

Author: AndiMD

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?