spruced up the extra bridge BC constructors to work to N dimensions

Author: xtalax

src/boundary_padded_arrays.jl

@@ -54,14 +54,14 @@ Composes BoundaryPaddedArrays that extend the same u for each different dimensio
 
 Ax Ay and Az can be passed in any order, as long as there is exactly one BoundaryPaddedArray that extends each dimension.
 """
-function compose(padded_arrays::BoundaryPaddedArray ...)
+function compose(padded_arrays::BoundaryPaddedArray...)
     N = ndims(padded_arrays[1])
     Ds = getaxis.(padded_arrays)
-    (length(padded_arrays) == N) || throw("The padded_arrays must cover every dimension - make sure that the number of padded_arrays is equal to ndims(u).")
+    (length(padded_arrays) == N) || throw(ArgumentError("The padded_arrays must cover every dimension - make sure that the number of padded_arrays is equal to ndims(u)."))
     for D in Ds
-        length(setdiff(Ds, D)) == (N-1) || throw("There are multiple Arrays that extend along $D - make sure every dimension has a unique extension")
+        length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple Arrays that extend along dimension $D - make sure every dimension has a unique extension"))
     end
-    reduce((|), fill(padded_arrays[1].u, (length(padded_arrays),)) .== getfield.(padded_arrays, :u)) || throw("The padded_arrays do not all extend the same u!")
+    reduce((|), fill(padded_arrays[1].u, (length(padded_arrays),)) .== getfield.(padded_arrays, :u)) || throw(ArgumentError("The padded_arrays do not all extend the same u!"))
     padded_arrays = padded_arrays[sortperm([Ds...])]
     lower = [padded_array.lower for padded_array in padded_arrays]
     upper = [padded_array.upper for padded_array in padded_arrays]

src/derivative_operators/BC_operators.jl

@@ -214,7 +214,7 @@ function BridgeBC(u_low::AbstractArray{T,N}, indslow, u_up::AbstractArray{T,N},
 end
 
 """
-    Q1, Q2 = BridgeBC(u1::AbstractVector{T}, hilo1::String, bc1::AtomicBC{T}, u2::AbstractVector{T}, hilo2::AbstractVector{T}, bc2::AtomicBC{T})
+    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.
@@ -225,18 +225,20 @@ When using these with a time/space stepping solve, please use elementwise equals
     u_t1 .= L*Q*u_t0
 -----------------------------------------------------------------------------------
 Connecting two multi dimensional Arrays:
-    Q1, Q2 = BridgeBC(u1::AbstractArray{T,N}, dim1::Int, hilo1::String, bc1, u2::AbstractArray{T,N}, dim2::Int, hilo2::String, bc2)
+    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 BridgeBC(u1::AbstractVector{T}, hilo1::String, bc1::AtomicBC{T}, u2::AbstractVector{T}, hilo2::AbstractVector{T}, bc2::AtomicBC{T}) where T
+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"

src/derivative_operators/multi_dim_bc_operators.jl

@@ -12,6 +12,7 @@ abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end
 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]
@@ -30,6 +31,7 @@ end
 end
 
 function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) where {T, B, N,M}
+    @assert N != 1
     lower = zeros(T,perpsize(u,dim))
     upper = zeros(T,perpsize(u,dim))
     pre = axes(u)[1:dim-1]
@@ -52,8 +54,8 @@ end
 
 """
 A multiple dimensional BC, supporting arbitrary BCs at each boundary point.
-To construct an arbitrary BC, pass an Array of BCs with dimension one less than that of your domain u - denoted N,
-with a size of size(u)[setdiff(1:N, dim)], where dim is the dimension orthogonal to the boundary that you want to extend.
+To construct an arbitrary BC, pass an Array of BCs with dimension `N-1`, if `N` is the dimensionality of your domain `u`
+with a size of `size(u)[setdiff(1:N, dim)]`, where dim is the dimension orthogonal to the boundary that you want to extend.
 
 It is also possible to call
     Q_dim = MultiDimBC(YourBC, size(u), dim)
@@ -69,7 +71,7 @@ In the case where you want to extend the same Robin/GeneralBC to the whole bound
 or
     Qx, Qy, Qz... = GeneralBC(αl, αr, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
 
-There are also constructors for NeumannBC, DirichletBC, Neumann0BC and Dirichlet0BC. Simply replace dx with the tuple as above, and append size(u) to the argument signature.
+There are also constructors for NeumannBC, DirichletBC, Neumann0BC and Dirichlet0BC. Simply replace `dx` in the call with the tuple dxyz... as above, and append `size(u)`` to the argument signature.
 The order is a required argument in this case.
 
 where dx, dy, and dz are vectors of grid steps.
@@ -99,114 +101,72 @@ 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(u1::AbstractArray{T,2}, dim1::Int, hilo1::String, bc1::MultiDimDirectionalBC, u2::AbstractArray{T,2}, dim2::Int, hilo2::String, bc2::MultiDimDirectionalBC) where {T}
-    @assert 1 ≤ dim1 ≤ 2 "dim1 must be 1≤dim1≤2, got dim1 = $dim1"
-    @assert 1 ≤ dim1 ≤ 2 "dim2 must be 1≤dim1≤2, got dim1 = $dim1"
+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, 2, eltype(s1)}, getboundarytype(bc1)}}(undef, s1...)
-            BC2 = Array{MixedBC{T, BridgeBC{T, 2, eltype(s2)}, getboundarytype(bc2)}}(undef, s2...)
+            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)
-            for i in 1:s1[1]
-                BC1[i] = MixedBC(BridgeBC{T, 2, eltype(s1)}(zeros(T, 1), view(view2, i), zeros(T, 1), view(view2, i)), bc1.BCs[i])
-                BC2[i] = MixedBC(BridgeBC{T, 2, eltype(s1)}(zeros(T, 1), view(view1, i), zeros(T, 1), view(view1, i)), bc2.BCs[i])
+            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, 2, eltype(s1)}, getboundarytype(bc1)}}(undef, s1...)
-            BC2 = Array{MixedBC{T, getboundarytype(bc2), BridgeBC{T, 2, eltype(s2)}}}(undef, s2...)
+            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])
-            for i in 1:s1[1]
-                BC1[i] = MixedBC(BridgeBC{T, 2, 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, 2, eltype(s1)}(zeros(T, 1), view(view1, i), zeros(T, 1), view(view1, i)))
+            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("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")
+            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, 2, eltype(s1)}}}(undef, s1...)
-            BC2 = Array{MixedBC{T, BridgeBC{T, 2, eltype(s2)}, getboundarytype(bc2)}}(undef, s2...)
+            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)
-            for i in 1:s1[1]
-                BC1[i] = MixedBC(bc1.BCs[i], BridgeBC{T, 2, eltype(s1)}(zeros(T, 1), view(view2, i), zeros(T, 1), view(view2, i)))
-                BC2[i] = MixedBC(BridgeBC{T, 2, eltype(s1)}(zeros(T, 1), view(view1, i), zeros(T, 1), view(view1, i)), bc2.BCs[i])
+            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, 2, eltype(s1)}}}(undef, s1...)
-            BC2 = Array{MixedBC{T, getboundarytype(bc2), BridgeBC{T, 2, eltype(s2)}}}(undef, s2...)
+            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])
-            for i in 1:s1[1]
-                BC1[i] = MixedBC(bc1.BCs[i], BridgeBC{T, 2, eltype(s1)}(zeros(T, 1), view(view2, i), zeros(T, 1), view(view2, i)))
-                BC2[i] = MixedBC(bc2.BCs[i], BridgeBC{T, 2, eltype(s1)}(zeros(T, 1), view(view1, i), zeros(T, 1), view(view1, i)))
+            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("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")
+            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
 
-function BridgeBC(u1::AbstractArray{T,3}, dim1::Int, hilo1::String, bc1::MultiDimDirectionalBC, u2::AbstractArray{T,3}, dim2::Int, hilo2::String, bc2::MultiDimDirectionalBC) where {T}
-    @assert 1 ≤ dim1 ≤ 3 "dim1 must be 1≤dim1≤3, got dim1 = $dim1"
-    @assert 1 ≤ dim1 ≤ 3 "dim2 must be 1≤dim1≤3, 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, 3, eltype(s1)}, getboundarytype(bc1)}}(undef, s1...)
-            BC2 = Array{MixedBC{T, BridgeBC{T, 3, eltype(s2)}, getboundarytype(bc2)}}(undef, s2...)
-            view2 = selectdim(u2, dim2, 1)
-            for j in 1:s1[2], i in 1:s1[1]
-                BC1[i, j] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, i, j), zeros(T, 1), view(view2, i, j)), bc1.BCs[i,j])
-                BC2[i, j] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, i, j), zeros(T, 1), view(view1, i, j)), bc2.BCs[i,j])
-            end
-        elseif hilo2 == "high"
-            BC1 = Array{MixedBC{T, BridgeBC{T, 3, eltype(s1)}, getboundarytype(bc1)}}(undef, s1...)
-            BC2 = Array{MixedBC{T, getboundarytype(bc2), BridgeBC{T, 3, eltype(s2)}}}(undef, s2...)
-            view2 = selectdim(u2, dim2, size(u2)[dim2])
-            for j in 1:s1[2], i in 1:s1[1]
-                BC1[i, j] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, i, j), zeros(T, 1), view(view2, i, j)), bc1.BCs[i,j])
-                BC2[i, j] = MixedBC(bc2.BCs[i,j], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, i, j), zeros(T, 1), view(view1, i, j)))
-            end
-        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 = selectdim(u1, dim1, size(u1)[dim1])
-        if hilo2 == "low"
-            BC1 = Array{MixedBC{T, getboundarytype(bc1), BridgeBC{T, 3, eltype(s1)}}}(undef, s1...)
-            BC2 = Array{MixedBC{T, BridgeBC{T, 3, eltype(s2)}, getboundarytype(bc2)}}(undef, s2...)
-            view2 = selectdim(u2, dim2, 1)
-            view2 = selectdim(u2, dim2, 1)
-            for j in 1:s1[2], i in 1:s1[1]
-                BC1[i, j] = MixedBC(bc1.BCs[i,j], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, i, j), zeros(T, 1), view(view2, i, j)))
-                BC2[i, j] = MixedBC(BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, i, j), zeros(T, 1), view(view1, i, j)), bc2.BCs[i,j])
-            end
-        elseif hilo2 == "high"
-            BC1 = Array{MixedBC{T, getboundarytype(bc1), BridgeBC{T, 3, eltype(s1)}}}(undef, s1...)
-            BC2 = Array{MixedBC{T, getboundarytype(bc2), BridgeBC{T, 3, eltype(s2)}}}(undef, s2...)
-            view2 = selectdim(u2, dim2, size(u2)[dim2])
-            for j in 1:s1[2], i in 1:s1[1]
-                BC1[i, j] = MixedBC(bc1.BCs[i,j], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view2, i, j), zeros(T, 1), view(view2, i, j)))
-                BC2[i, j] = MixedBC(bc2.BCs[i,j], BridgeBC{T, N, eltype(s1)}(zeros(T, 1), view(view1, i, j), zeros(T, 1), view(view1, i, j)))
-            end
-        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 (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...)
 
@@ -224,9 +184,9 @@ function compose(BCs...)
     T = gettype(BCs[1])
     N = ndims(BCs[1])
     Ds = getaxis.(BCs)
-    (length(BCs) == N) || throw("There must be enough BCs to cover every dimension - check that the number of MultiDimBCs == N")
+    (length(BCs) == N) || throw(ArgumentError("There must be enough BCs to cover every dimension - check that the number of MultiDimBCs == N"))
     for D in Ds
-        length(setdiff(Ds, D)) == (N-1) || throw("There are multiple boundary conditions that extend along $D - make sure every dimension has a unique extension")
+        length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple boundary conditions that extend along $D - make sure every dimension has a unique extension"))
     end
     BCs = BCs[sortperm([Ds...])]
 

test/MultiDimBC_test.jl

@@ -4,8 +4,8 @@ using LinearAlgebra, DiffEqOperators, Random, Test
 ################################################################################
 
 #Create Array
-n = 100
-m = 120
+n = 12
+m = 25
 A = rand(n,m)
 
 #Create atomic BC
@@ -38,9 +38,9 @@ end
 ################################################################################
 
 #Create Array
-n = 100
-m = 120
-o = 78
+n = 31
+m = 25
+o = 12
 A = rand(n,m, o)
 
 #Create atomic BC
@@ -60,19 +60,8 @@ Ay = Qy*A
 Az = Qz*A
 # Test padded array compositions
 
-Aextended = compose(Az,Ax,Ay) #Arrays in wierd order to show that this still works
-@test_broken compose(Ax, Az, Az) #these tests should not work, if they are broken then they are  passing
-
-#test BC compositions
-Q = compose(Qy, Qx,Qz)
-@test_broken compose(Qx, Qx, Qz) #these tests should not work, if they are broken then they are passing
-
 QA = Q*A
 
-for i in 1:(n+2), j in 1:(m+2), k in 1:(o+2)
-    @test Aextended[i,j,k] == QA[i,j,k]
-end
-
 @test size(Ax)[1] == size(A)[1]+2
 @test size(Ay)[2] == size(A)[2]+2
 @test size(Az)[3] == size(A)[3]+2
@@ -85,3 +74,18 @@ end
 for i in 1:n, j in 1:m
     @test Az[i, j, :] == Array(BCz[i, j]*A[i, j, :])
 end
+
+#test compositions to higher dimension
+for N in 2:7
+    sizes = rand(5:10)
+    A = rand(sizes...)
+
+    Q1_N = Neumann0BC(1.0, 3.0, size(A))
+
+    Q = compose(Q1_N...)
+
+    A1_N = Q1_N.*fill(A, N)
+
+    A_extended = Q*A
+    @test Array(A_extended) == Array(compose(A1_N...))
+end