Merge branch 'bridge-bc' into SplitBC-dev

Author: xtalax

src/derivative_operators/BC_operators.jl

@@ -165,6 +165,94 @@ 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, L} <: 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},L}
+    a_r::Vector{T}
+    b_r::SubArray{T,0,Array{T,N},NTuple{N,I},L}
+    function BridgeBC{T,N,L}(a_l, b_l::Subarray{T,0,Array{T,N},NTuple{N,I}, L}, a_r, b_l::Subarray{T,0,Array{T,N},NTuple{N,I}, L}) where {T, N, L, I} = new{T,N,L,I}(a_l,b_l,a_r,b_r)
+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::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,16 @@ end
 function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
     @assert N != 1
     u_temp = similar(u)
+
     _slice_rmul!(u_temp, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))
+
     return u_temp
 end
 
 @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 post
         for I in pre
+
             tmp = A[I,J]*u[I, :, J]
             lower[I,J], upper[I,J] = tmp.l, tmp.r
         end
@@ -33,7 +36,9 @@ function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) wher
     @assert M == N-1
     lower = zeros(T,perpsize(u,dim))
     upper = zeros(T,perpsize(u,dim))
+
     _slice_rmul!(lower, upper, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))
+
     return (lower, upper)
 end
 

test/runtests.jl

@@ -4,6 +4,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 "JacVec Operators Interface" begin include("jacvec_operators.jl") end
+@time @safetestset "Validate BridgeBC and its constructors" begin include("bridge_bc.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
 @time @safetestset "Derivative Operators Interface" begin include("derivative_operators_interface.jl") end