Got tests passing

Author: xtalax

src/boundary_padded_arrays.jl

@@ -100,14 +100,16 @@ decompose(A::ComposedBoundaryPaddedArray) = Tuple([BoundaryPaddedArray{gettype(A
 
 
 Base.length(Q::AbstractBoundaryPaddedArray) = reduce((*), size(Q))
-Base.lastindex(Q::AbstractBoundaryPaddedArray) = Base.length(Q)
+Base.firstindex(Q::AbstractBoundaryPaddedArray, d::Int) = 1
+Base.lastindex(Q::AbstractBoundaryPaddedArray) = length(Q)
+Base.lastindex(Q::AbstractBoundaryPaddedArray, d::Int) = size(Q)[d]
 gettype(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = T
 Base.ndims(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = N
 
 add_dim(A::AbstractArray, i) = reshape(A, size(A)...,i)
 add_dim(i) = i
 
-function Base.getindex(Q::BoundaryPaddedArray{T,D,N,M,V,B}, _inds::NTuple{N,Int64}...) where {T,D,N,M,V,B} #supports range and colon indexing!
+function Base.getindex(Q::BoundaryPaddedArray{T,D,N,M,V,B}, _inds::Vararg{Int,N}) where {T,D,N,M,V,B} #supports range and colon indexing!
     inds = [_inds...]
     S = size(Q)
     dim = D
@@ -131,7 +133,7 @@ function Base.getindex(Q::BoundaryPaddedArray{T,D,N,M,V,B}, _inds::NTuple{N,Int6
     end
 end
 
-function Base.getindex(Q::ComposedBoundaryPaddedArray{T, N, M, V, B} , inds::NTuple{N, Int64}...) where {T, N, M, V, B} #as yet no support for range indexing or colon indexing
+function Base.getindex(Q::ComposedBoundaryPaddedArray{T, N, M, V, B} , inds::Vararg{Int, N}) where {T, N, M, V, B} #as yet no support for range indexing or colon indexing
     S = size(Q)
     @assert reduce((&), inds .<= S)
     for (dim, index) in enumerate(inds)
@@ -153,5 +155,3 @@ function Base.getindex(Q::ComposedBoundaryPaddedArray{T, N, M, V, B} , inds::NTu
      end
     return Q.u[(inds.-1)...]
 end
-
-getindex(A::AbstractBoundaryPaddedArray{T,N}, I::CartesianIndex{N}) where {T,N} = A[Tuple(I)...]

src/derivative_operators/multi_dim_bc_operators.jl

@@ -3,8 +3,8 @@ abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end
 
 
 @noinline function _slice_rmul!(u_temp::AbstractArray{T,N}, A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,N}
-    for J in CartesianIndices(post)
-        for I in CartesianIndices(pre)
+    for J in post
+        for I in pre
             u_temp[I, :, J] = A*u[I, :, J]
         end
     end
@@ -14,13 +14,13 @@ 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, axes(u)[1:dim-1], axes(u)[dim+1:end])
+    _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 CartesianIndices(post)
-        for I in CartesianIndices(pre)
+    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 +33,7 @@ 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, axes(u)[1:dim-1], axes(u)[dim+1:end])
+    _slice_rmul!(lower, upper, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))
     return (lower, upper)
 end
 
@@ -68,10 +68,14 @@ 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` in the call with the tuple dxyz... as above, and append `size(u)`` to the argument signature.
+There are also constructors for NeumannBC, DirichletBC 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.
+
+For Neumann0BC, please use
+    Qx, Qy, Qz... = Neumann0BC(T::Type, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
+where T is the element type of the domain to be extended
 """
 MultiDimBC(BC::Array{B,N}, dim::Integer) where {N, B<:AtomicBC} = MultiDimDirectionalBC{gettype(BC[1]), B, dim, N+1, N}(BC)
 #s should be size of the domain
@@ -138,11 +142,15 @@ getboundarytype(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} =
 Base.ndims(Q::MultiDimensionalBC{T,N}) where {T,N} = N
 
 function Base.:*(Q::MultiDimDirectionalBC{T, B, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
+    @assert perpsize(u, D) == size(Q.BCs) "Size of the BCs array in the MultiDimBC is incorrect, needs to be $(perpsize(u,D)) to extend dimension $D, got $(size(Q.BCs))"
     lower, upper = slice_rmul(Q.BCs, u, D)
     return BoundaryPaddedArray{T, D, N, K, typeof(u), typeof(lower)}(lower, upper, u)
 end
 
 function Base.:*(Q::ComposedMultiDimBC{T, B, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
+    for dim in 1:N
+        @assert perpsize(u, dim) == size(Q.BCs[dim]) "Size of the BCs array for dimension $dim in the MultiDimBC is incorrect, needs to be $(perpsize(u,dim)), got $(size(Q.BCs[dim]))"
+    end
     out = slice_rmul.(Q.BCs, fill(u, N), 1:N)
     return ComposedBoundaryPaddedArray{T, N, K, typeof(u), typeof(out[1][1])}([A[1] for A in out], [A[2] for A in out], u)
 end

test/MultiDimBC_test.jl

@@ -4,16 +4,16 @@ using LinearAlgebra, DiffEqOperators, Random, Test
 ################################################################################
 
 #Create Array
-n = 12
-m = 25
+n = 8
+m = 15
 A = rand(n,m)
 
 #Create atomic BC
 q1 = RobinBC([1.0, 2.0, 3.0], [0.0, -1.0, 2.0], 0.1, 4.0)
 q2 = PeriodicBC{Float64}()
 
-BCx = vcat(fill(q1, div(m,2)), fill(q2, div(m,2)))  #The size of BCx has to be all size components *except* for x
-BCy = vcat(fill(q1, div(n,2)), fill(q2, div(n,2)))
+BCx = vcat(fill(q1, div(m,2)), fill(q2, m-div(m,2)))  #The size of BCx has to be all size components *except* for x
+BCy = vcat(fill(q1, div(n,2)), fill(q2, n-div(n,2)))
 
 
 Qx = MultiDimBC(BCx, 1)
@@ -38,29 +38,26 @@ end
 ################################################################################
 
 #Create Array
-n = 31
-m = 25
+n = 8
+m = 11
 o = 12
 A = rand(n,m, o)
 
 #Create atomic BC
 q1 = RobinBC([1.0, 2.0, 3.0], [0.0, -1.0, 2.0], 0.1, 4.0)
 q2 = PeriodicBC{Float64}()
 
-BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (div(m,2), o)))  #The size of BCx has to be all size components *except* for x
-BCy = vcat(fill(q1, (div(n,2), o)), fill(q2, (div(n,2), o)))
-BCz = fill(MixedBC(q1,q2), (n,m))
+BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (m-div(m,2), o)))  #The size of BCx has to be all size components *except* for x
+BCy = vcat(fill(q1, (div(n,2), o)), fill(q2, (n-div(n,2), o)))
+BCz = fill(Dirichlet0BC(Float64), (n,m))
 
 Qx = MultiDimBC(BCx, 1)
 Qy = MultiDimBC(BCy, 2)
-Qz = MultiDimBC(MixedBC(q1, q2), size(A), 3) #Test the other constructor
+Qz = MultiDimBC(Dirichlet0BC(Float64), size(A), 3) #Test the other constructor
 
 Ax = Qx*A
 Ay = Qy*A
 Az = Qz*A
-# Test padded array compositions
-
-QA = Q*A
 
 @test size(Ax)[1] == size(A)[1]+2
 @test size(Ay)[2] == size(A)[2]+2
@@ -77,10 +74,10 @@ end
 
 #test compositions to higher dimension
 for N in 2:7
-    sizes = rand(5:10)
+    sizes = rand(4:7, N)
     A = rand(sizes...)
 
-    Q1_N = Neumann0BC(1.0, 3.0, size(A))
+    Q1_N = Neumann0BC(Float64, Tuple(ones(N)), 3.0, size(A))
 
     Q = compose(Q1_N...)
 

test/boundary_padded_array.jl

@@ -5,10 +5,10 @@ using LinearAlgebra, DiffEqOperators, Random, Test
 
 for dimensionality in 2:5
     for dim in 1:dimensionality
-        sizes = rand(10:20, dimensionality)
+        sizes = rand(4:10, dimensionality)
         A = rand(sizes...)
-        lower = selectdim(A, dim, 1)
-        upper = selectdim(A, dim, size(A)[dim])
+        lower = Array(selectdim(A, dim, 1))
+        upper = Array(selectdim(A, dim, size(A)[dim]))
 
         Apad = DiffEqOperators.BoundaryPaddedArray{Float64, dim, dimensionality, dimensionality-1, typeof(A), typeof(lower)}(lower, upper, selectdim(A, dim, 2:(size(A)[dim]-1)))
 
@@ -24,15 +24,15 @@ end
 # Test ComposedBoundaryPaddedMatrix
 ################################################################################
 
-n = 10
-m = 21
+n = 5
+m = 7
 A = rand(n,m)
 A[1,1] = A[end,1] = A[1,end] = A[end,end] = 0.0
 
 lower = Vector[A[1,2:(end-1)], A[2:(end-1),1]]
 upper = Vector[A[end,2:(end-1)], A[2:(end-1),end]]
 
-Apad = ComposedBoundaryPaddedMatrix{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1)])
+Apad = DiffEqOperators.ComposedBoundaryPaddedMatrix{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1)])
 
 @test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix
 
@@ -44,18 +44,18 @@ end
 # Test ComposedBoundaryPaddedTensor{3}
 ################################################################################
 
-n = 10
-m = 12
+n = 4
+m = 5
 o = 7
 A = rand(n,m,o)
-A[1,1,:] = A[end,1, :] = A[1,end, :] = A[end,end, :] = 0.0
-A[1,:,1] = A[end, :, 1] = A[1,:,end] = A[end,:,end] = 0.0
-A[:,1,1] = A[:,end,1] = A[:,1,end] = A[:,end,end] = 0.0
+A[1,1,:] = A[end,1, :] = A[1,end, :] = A[end,end, :] = zeros(o)
+A[1,:,1] = A[end, :, 1] = A[1,:,end] = A[end,:,end] = zeros(m)
+A[:,1,1] = A[:,end,1] = A[:,1,end] = A[:,end,end] = zeros(n)
 
-lower = Matrix[A[1,2:(end-1), 2:(end-1)], A[2:(end-1),1, 2:(end-1)] A[2:(end-1), 2:(end-1), 1]]
-upper = Matrix[A[end, 2:(end-1), 2:(end-1)], A[2:(end-1), end, 2:(end-1)] A[2:(end-1), 2:(end-1), end]]
+lower = Matrix[A[1,2:(end-1), 2:(end-1)], A[2:(end-1),1, 2:(end-1)], A[2:(end-1), 2:(end-1), 1]]
+upper = Matrix[A[end, 2:(end-1), 2:(end-1)], A[2:(end-1), end, 2:(end-1)], A[2:(end-1), 2:(end-1), end]]
 
-Apad = BoundaryPadded3Tensor{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1), 2:(end-1)])
+Apad = DiffEqOperators.ComposedBoundaryPadded3Tensor{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1), 2:(end-1)])
 
 @test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix