Pulled out BridgeBC in to its own branch, since it was causing problems and is not strictly nessecary

Author: xtalax

src/DiffEqOperators.jl

@@ -51,7 +51,7 @@ export MatrixFreeOperator
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
        CenteredDifference, UpwindDifference
-export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MixedBC, MultiDimBC, PeriodicBC, BridgeBC,
+export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MixedBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC,
        compose, decompose, perpsize
 

src/derivative_operators/BC_operators.jl

@@ -153,8 +153,6 @@ struct GeneralBC{T} <:AffineBC{T}
     end
 end
 
-
-
 """
 Quick and dirty way to allow mixed boundary types on each end of an array - may be cleaner and more versatile to split up left and right boundaries going forward
 MixedBC(lowerBC, upperBC) is the interface.
@@ -170,12 +168,6 @@ struct MixedBC{T, R, S} <: AtomicBC{T}
      end
 end
 
-function Base.:*(Q::MixedBC, u::AbstractVector)
-    lower = Q.lower*u
-    upper = Q.upper*u
-    return BoundaryPaddedVector(lower.l, upper.r, u)
-end
-
 #implement Neumann and Dirichlet as special cases of RobinBC
 NeumannBC(α::AbstractVector{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], 1.0, 2.0 )
@@ -186,92 +178,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)
 
-"""
-BridgeBC(u_low::AbstractArray{T,N}, u_up::AbstractArray{T,N}, indslow, indsup) # A different view in to 2 diffferent arrays on each end of the boundary, indslow is an iterable of indicies that index u_low, which extends the lower index end. Analogous for u_up and indsup with the upper boundary.
-
-BridgeBC(u::AbstractArray{T,N}, inds) # The same view in to some array u at the index inds extends the boundary
-
--------------------------------------------------------------------------------------
-
-Allows seperate domains governed by seperate equations to be bridged together with a boundary condition.
-"""
-struct BridgeBC{T,N,I} <: AffineBC{T}
-    a_l::Vector{T} #Dummy vectors so that AffineBC methods still work
-    b_l::SubArray{T,0,Array{T,N},NTuple{N,I},true}
-    a_r::Vector{T}
-    b_r::SubArray{T,0,Array{T,N},NTuple{N,I},true}
-end
-
-BridgeBC(u::AbstractArray, inds) = BridgeBC(u, inds, u, inds)
-
-function BridgeBC(u_low::AbstractArray{T,N}, indslow, u_up::AbstractArray{T,N},  indsup) where {T, N}
-    @assert length(indslow) == N
-    @assert length(indsup) == N
-    @assert mapreduce(x -> typeof(x) <: Integer, (&), indslow)
-    @assert mapreduce(x -> typeof(x) <: Integer, (&), indsup)
-
-    BridgeBC{T, length(indslow), eltype(indslow)}(zeros(T,1), view(u_low, indslow...), zeros(T,1), view(u_up, indsup...))
-end
-
-"""
-    Q1, Q2 = BridgeBC(bc1::AtomicBC, u1::AbstractArray{T,1}, hilo1::String, hilo2::String, u2::AbstractArray{T,1}, bc2::AtomicBC)
--------------------------------------------------------------------------------------
-Creates two BC operators that join array `u1` to `u2` at the `hilo1` end ("high" or "low" index end), and joins `u2` to `u1` with simalar settings given in `hilo2`.
-The ends of `u1` and `u2` that are not connected will use the boundary conditions `bc1` and `bc2` respectively.
-
-Use `Q1` to extend `u1` and `Q2` to extend `u2`.
-
-When using these with a time/space stepping solve, please use elementwise equals on your u1 and u2 to avoid the need to create new BC operators each time, as follows:
-    u_t1 .= L*Q*u_t0
------------------------------------------------------------------------------------
-Connecting two multi dimensional Arrays:
-    Q1, Q2 = BridgeBC(bc1::MultiDimDirectionalBC, u1::AbstractArray{T,N}, dim1::Int, hilo1::String, dim2::Int, hilo2::String, u2::AbstractArray{T,N}, bc2::MultiDimDirectionalBC)
------------------------------------------------------------------------------------
-
-Creates two BC operators that join array `u1` to `u2` at the `hilo1` end ("high" or "low" index end) of dimension `dim1`, and joins `u2` to `u1` with simalar settings given in `hilo2` and `dim2`.
-The ends of `u1` and `u2` that are not connected will use the boundary conditions `bc1` and `bc2` respectively.
-
-Drop `dim1` and `dim2` when your `u1` and `u2` are vectors.
-
-Use `Q1` to extend `u1` and `Q2` to extend `u2`.
-
-When using these with a time/space stepping solve, please use elementwise equals on your u1 and u2 to avoid the need to create new BC operators each time, as follows:
-    u_t1 .= L*Q*u_t0
-"""
-function BridgeBCBridgeBC(bc1::AtomicBC, u1::AbstractArray{T,1}, hilo1::String, hilo2::String, u2::AbstractArray{T,1}, bc2::AtomicBC) where T
-    if hilo1 == "low"
-        view1 = view(u1, 1)
-        if hilo2 == "low"
-            view2 = view(u2, 1)
-            BC1 = MixedBC(BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view2, zeros(T, 1), view2), bc1)
-            BC2 = MixedBC(BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view1, zeros(T, 1), view1), bc2)
-        elseif hilo2 == "high"
-            view2 = view(u2, length(u2))
-            BC1 = MixedBC(BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view2, zeros(T, 1), view2), bc1)
-            BC2 = MixedBC(bc2, BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view1, zeros(T, 1), view1))
-        else
-            throw("hilo2 not recognized, please use \"high\" to connect u1 to u2 along the upper index of dim2 of u2 or \"low\" to connect along the lower index end")
-        end
-    elseif hilo1 == "high"
-        view1 = view(u1, length(u1))
-        if hilo2 == "low"
-            view2 = view(u2, 1)
-            BC1 = MixedBC(bc1, BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view2, zeros(T, 1), view2))
-            BC2 = MixedBC(BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view1, zeros(T, 1), view1), bc2)
-        elseif hilo2 == "high"
-            view2 = view(u2, length(u2))
-            BC1 = MixedBC(bc1, BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view2, zeros(T, 1), view2))
-            BC2 = MixedBC(bc2, BridgeBC{T, 1, eltype(s1)}(zeros(T, 1), view1, zeros(T, 1), view1))
-        else
-            throw("hilo2 not recognized, please use \"high\" to connect u1 to u2 along the upper index of dim2 of u2 or \"low\" to connect along the lower index end")
-        end
-    else
-        throw("hilo1 not recognized, please use \"high\" to connect u1 to u2 along the upper index of dim1 of u1 or \"low\" to connect along the lower index end")
-    end
-    return (BC1, BC2)
-end
-
-Base.:*(Q::BridgeBC{T,I,N}, u::AbstractVector{T}) where {T, I, N} = BoundaryPaddedVector{T, typeof(u)}(Q.b_l[1], Q.b_r[1], u)
+Base.:*(Q::MixedBC, u::AbstractVector) =  BoundaryPaddedVector((Q.lower*u).l, (Q.upper*u).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)
 

src/derivative_operators/multi_dim_bc_operators.jl

@@ -14,13 +14,11 @@ end
 function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
     @assert N != 1
     u_temp = similar(u)
-    pre = axes(u)[1:dim-1]
-    post = axes(u)[dim+1:end]
-    _slice_rmul!(u_temp, A, u, dim, pre, post)
+    _slice_rmul!(u_temp, A, u, dim, axes(u)[1:dim-1], axes(u)[dim+1:end])
     return u_temp
 end
 
-@noinline function _slice_rmul!(lower::Array{T, M}, upper::Array{T, M}, A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,B,N,M}
+@noinline function _slice_rmul!(lower::AbstractArray, upper::AbstractArray, A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,B,N,M}
     for J in CartesianIndices(post)
         for I in CartesianIndices(pre)
             tmp = A[I,J]*u[I, :, J]
@@ -32,11 +30,10 @@ end
 
 function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) where {T, B, N,M}
     @assert N != 1
+    @assert M == N-1
     lower = zeros(T,perpsize(u,dim))
     upper = zeros(T,perpsize(u,dim))
-    pre = axes(u)[1:dim-1]
-    post = axes(u)[dim+1:end]
-    _slice_rmul!(lower, upper, A, u, dim, pre, post)
+    _slice_rmul!(lower, upper, A, u, dim, axes(u)[1:dim-1], axes(u)[dim+1:end])
     return (lower, upper)
 end
 
@@ -100,73 +97,6 @@ GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where
 
 perpsize(A::AbstractArray{T,N}, dim::Integer) where {T,N} = size(A)[setdiff(1:N, dim)] #the size of A perpendicular to dim
 
-# Constructors for Bridge BC to make it easier to join domains together. See docs on BrigeBC in BC_operators.jl for info on usage
-function BridgeBC(bc1::MultiDimDirectionalBC, u1::AbstractArray{T,N}, dim1::Int, hilo1::String, dim2::Int, hilo2::String, u2::AbstractArray{T,N}, bc2::MultiDimDirectionalBC) where {T, N}
-    @assert 1 ≤ dim1 ≤ N "dim1 must be 1≤dim1≤N, got dim1 = $dim1"
-    @assert 1 ≤ dim1 ≤ N "dim2 must be 1≤dim1≤N, got dim1 = $dim1"
-    s1 = perpsize(u1, dim1) #
-    s2 = perpsize(u2, dim2)
-    @assert s1 == s2 "Arrays must be same size along boundary to be joined, got boundary sizes u1 = $s1, u2 = $s2"
-    if hilo1 == "low"
-        view1 = selectdim(u1, dim1, 1)
-        if hilo2 == "low"
-            BC1 = Array{MixedBC{T, BridgeBC{T, N, eltype(s1)}, getboundarytype(bc1)}}(undef, s1...)
-            BC2 = Array{MixedBC{T, BridgeBC{T, N, eltype(s2)}, getboundarytype(bc2)}}(undef, s2...)
-            view2 = selectdim(u2, dim2, 1)
-            R = CartesianIndices(BC1)
-            for I in R
-                BC1[I] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, I), zeros(T, 1), view(view2, I)), bc1.BCs[I])
-                BC2[I] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, I), zeros(T, 1), view(view1, I)), bc2.BCs[I])
-            end
-        elseif hilo2 == "high"
-            BC1 = Array{MixedBC{T, BridgeBC{T, N, eltype(s1)}, getboundarytype(bc1)}}(undef, s1...)
-            BC2 = Array{MixedBC{T, getboundarytype(bc2), BridgeBC{T, N, eltype(s2)}}}(undef, s2...)
-            view2 = selectdim(u2, dim2, size(u2)[dim2])
-            R = CartesianIndices(BC1)
-            for I in R
-                BC1[I] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, I), zeros(T, 1), view(view2, I)), bc1.BCs[I])
-                BC2[I] = MixedBC(bc2.BCs[I], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, I), zeros(T, 1), view(view1, I)))
-            end
-        else
-            throw(ArgumentError("hilo2 not recognized, please use \"high\" to connect u1 to u2 along the upper index of dim2 of u2 or \"low\" to connect along the lower index end"))
-        end
-    elseif hilo1 == "high"
-        view1 = selectdim(u1, dim1, size(u1)[dim1])
-        if hilo2 == "low"
-            BC1 = Array{MixedBC{T, getboundarytype(bc1), BridgeBC{T, N, eltype(s1)}}}(undef, s1...)
-            BC2 = Array{MixedBC{T, BridgeBC{T, N, eltype(s2)}, getboundarytype(bc2)}}(undef, s2...)
-            view2 = selectdim(u2, dim2, 1)
-            R = CartesianIndices(BC1)
-            for I in R
-                BC1[I] = MixedBC(bc1.BCs[I], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, I), zeros(T, 1), view(view2, I)))
-                BC2[I] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, I), zeros(T, 1), view(view1, I)), bc2.BCs[I])
-            end
-        elseif hilo2 == "high"
-            BC1 = Array{MixedBC{T, getboundarytype(bc1), BridgeBC{T, N, eltype(s1)}}}(undef, s1...)
-            BC2 = Array{MixedBC{T, getboundarytype(bc2), BridgeBC{T, N, eltype(s2)}}}(undef, s2...)
-            view2 = selectdim(u2, dim2, size(u2)[dim2])
-            R = CartesianIndices(BC1)
-            for I in R
-                BC1[I] = MixedBC(bc1.BCs[I], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, I), zeros(T, 1), view(view2, I)))
-                BC2[I] = MixedBC(bc2.BCs[I], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, I), zeros(T, 1), view(view1, I)))
-            end
-        else
-            throw(ArgumentError("hilo2 not recognized, please use \"high\" to connect u1 to u2 along the upper index of dim2 of u2 or \"low\" to connect along the lower index end"))
-        end
-    else
-        throw("hilo1 not recognized, please use \"high\" to connect u1 to u2 along the upper index of dim1 of u1 or \"low\" to connect along the lower index end")
-    end
-    return (MultiDimBC(BC1, dim1), MultiDimBC(BC2, dim2))
-end
-
-"""
-    Q1, Q2 = BridgeBC(bc1::MultiDimDirectionalBC, u1::AbstractArray, dim1::Int, dim2::Int, u2::AbstractArray, bc2::MultiDimDirectionalBC)
-Create BC operators that connect `u1` to `u2`, at the high index end of `dim1` for `u1`, and the low index end of `dim2` for `u2`.
-`bc1` and `bc2` are then the boundary conditions at the unconnected ends of `u1` and `u2` respectively.
-"""
-BridgeBC(bc1::MultiDimDirectionalBC, u1::AbstractArray, dim1::Int, dim2::Int, u2::AbstractArray, bc2::MultiDimDirectionalBC) = BridgeBC(u1, dim1, "high", bc1, u2, dim2, "low", bc2)
-
-
 """
 Q = compose(BCs...)
 

test/bridge_bc.jl

@@ -3,7 +3,7 @@ import Base: isapprox
 
 @time @safetestset "Basic Operators Interface" begin include("basic_operators_interface.jl") end
 @time @safetestset "Robin Boundary Condition Operators" begin include("robin.jl") end
-@time @safetestset "Validate BridgeBC and its constructors" begin include("bridge_bc.jl") end
+#@time @safetestset "Validate BridgeBC and its constructors" begin include("bridge_bc.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