|
| 1 | +# Boundary Padded Arrays |
| 2 | +abstract type AbstractBoundaryPaddedArray{T, N} <: AbstractArray{T, N} end |
| 3 | + |
| 4 | +""" |
| 5 | +A vector type that extends a vector u with ghost points at either end |
| 6 | +""" |
| 7 | +struct BoundaryPaddedVector{T,T2 <: AbstractVector{T}} <: AbstractBoundaryPaddedArray{T, 1} |
| 8 | + l::T |
| 9 | + r::T |
| 10 | + u::T2 |
| 11 | +end |
| 12 | +Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2 |
| 13 | +Base.size(Q::BoundaryPaddedVector) = (length(Q.u) + 2,) |
| 14 | +Base.lastindex(Q::BoundaryPaddedVector) = Base.length(Q) |
| 15 | + |
| 16 | +function Base.getindex(Q::BoundaryPaddedVector,i) |
| 17 | + if i == 1 |
| 18 | + return Q.l |
| 19 | + elseif i == length(Q) |
| 20 | + return Q.r |
| 21 | + else |
| 22 | + return Q.u[i-1] |
| 23 | + end |
| 24 | +end |
| 25 | + |
| 26 | +""" |
| 27 | +Higher dimensional generalization of BoundaryPaddedVector, pads an array of dimension N along the dimension D with 2 Arrays of dimension N-1, stored in lower and upper |
| 28 | +
|
| 29 | +""" |
| 30 | +struct BoundaryPaddedArray{T, D, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T,N} |
| 31 | + lower::B #an array of dimension M = N-1, used to extend the lower index boundary |
| 32 | + upper::B #Ditto for the upper index boundary |
| 33 | + u::V |
| 34 | +end |
| 35 | + |
| 36 | +getaxis(Q::BoundaryPaddedArray{T,D,N,M,V,B}) where {T,D,N,M,V,B} = D |
| 37 | + |
| 38 | +function Base.size(Q::BoundaryPaddedArray) |
| 39 | + S = [size(Q.u)...] |
| 40 | + S[getaxis(Q)] += 2 |
| 41 | + return Tuple(S) |
| 42 | +end |
| 43 | + |
| 44 | +""" |
| 45 | +A = compose(padded_arrays::BoundaryPaddedArray...) |
| 46 | +
|
| 47 | +------------------------------------------------------------------------------------- |
| 48 | +
|
| 49 | +Example: |
| 50 | +A = compose(Ax, Ay, Az) # 3D domain |
| 51 | +A = compose(Ax, Ay) # 2D Domain |
| 52 | +
|
| 53 | +Composes BoundaryPaddedArrays that extend the same u for each different dimension that u has in to a ComposedBoundaryPaddedArray |
| 54 | +
|
| 55 | +Ax Ay and Az can be passed in any order, as long as there is exactly one BoundaryPaddedArray that extends each dimension. |
| 56 | +""" |
| 57 | +function compose(padded_arrays::BoundaryPaddedArray...) |
| 58 | + N = ndims(padded_arrays[1]) |
| 59 | + Ds = getaxis.(padded_arrays) |
| 60 | + (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).")) |
| 61 | + for D in Ds |
| 62 | + 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")) |
| 63 | + end |
| 64 | + 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!")) |
| 65 | + padded_arrays = padded_arrays[sortperm([Ds...])] |
| 66 | + lower = [padded_array.lower for padded_array in padded_arrays] |
| 67 | + upper = [padded_array.upper for padded_array in padded_arrays] |
| 68 | + |
| 69 | + ComposedBoundaryPaddedArray{gettype(padded_arrays[1]),N,N-1,typeof(padded_arrays[1].u),typeof(lower[1])}(lower, upper, padded_arrays[1].u) |
| 70 | +end |
| 71 | + |
| 72 | +# Composed BoundaryPaddedArray |
| 73 | + |
| 74 | +struct ComposedBoundaryPaddedArray{T, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T, N} |
| 75 | + lower::Vector{B} |
| 76 | + upper::Vector{B} |
| 77 | + u::V |
| 78 | +end |
| 79 | + |
| 80 | +# Aliases |
| 81 | +AbstractBoundaryPaddedMatrix{T} = AbstractBoundaryPaddedArray{T,2} |
| 82 | +AbstractBoundaryPadded3Tensor{T} = AbstractBoundaryPaddedArray{T,3} |
| 83 | + |
| 84 | +BoundaryPaddedMatrix{T, D, V, B} = BoundaryPaddedArray{T, D, 2, 1, V, B} |
| 85 | +BoundaryPadded3Tensor{T, D, V, B} = BoundaryPaddedArray{T, D, 3, 2, V, B} |
| 86 | + |
| 87 | +ComposedBoundaryPaddedMatrix{T,V,B} = ComposedBoundaryPaddedArray{T,2,1,V,B} |
| 88 | +ComposedBoundaryPadded3Tensor{T,V,B} = ComposedBoundaryPaddedArray{T,3,2,V,B} |
| 89 | + |
| 90 | +Base.size(Q::ComposedBoundaryPaddedArray) = size(Q.u).+2 |
| 91 | + |
| 92 | +""" |
| 93 | +Ax, Ay,... = decompose(A::ComposedBoundaryPaddedArray) |
| 94 | +
|
| 95 | +------------------------------------------------------------------------------------- |
| 96 | +
|
| 97 | +Decomposes a ComposedBoundaryPaddedArray in to components that extend along each dimension individually |
| 98 | +""" |
| 99 | +decompose(A::ComposedBoundaryPaddedArray) = Tuple([BoundaryPaddedArray{gettype(A), ndims(A), ndims(A)-1, typeof(lower[1])}(A.lower[i], A.upper[i], A.u) for i in 1:ndims(A)]) |
| 100 | + |
| 101 | + |
| 102 | +Base.length(Q::AbstractBoundaryPaddedArray) = reduce((*), size(Q)) |
| 103 | +Base.firstindex(Q::AbstractBoundaryPaddedArray, d::Int) = 1 |
| 104 | +Base.lastindex(Q::AbstractBoundaryPaddedArray) = length(Q) |
| 105 | +Base.lastindex(Q::AbstractBoundaryPaddedArray, d::Int) = size(Q)[d] |
| 106 | +gettype(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = T |
| 107 | +Base.ndims(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = N |
| 108 | + |
| 109 | +add_dim(A::AbstractArray, i) = reshape(A, size(A)...,i) |
| 110 | +add_dim(i) = i |
| 111 | + |
| 112 | +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! |
| 113 | + inds = [_inds...] |
| 114 | + S = size(Q) |
| 115 | + dim = D |
| 116 | + otherinds = inds[setdiff(1:N, dim)] |
| 117 | + @assert length(S) == N |
| 118 | + if inds[dim] == 1 |
| 119 | + return Q.lower[otherinds...] |
| 120 | + elseif inds[dim] == S[dim] |
| 121 | + return Q.upper[otherinds...] |
| 122 | + elseif typeof(inds[dim]) <: Integer |
| 123 | + inds[dim] = inds[dim] - 1 |
| 124 | + return Q.u[inds...] |
| 125 | + elseif typeof(inds[dim]) == Colon |
| 126 | + if mapreduce(x -> typeof(x) != Colon, (|), otherinds) |
| 127 | + return vcat(Q.lower[otherinds...], Q.u[inds...], Q.upper[otherinds...]) |
| 128 | + else |
| 129 | + throw("A colon on the extended dim is as yet incompatible with additional colons") |
| 130 | + end |
| 131 | + elseif typeof(inds[dim]) <: AbstractArray |
| 132 | + throw("Range indexing not yet supported!") |
| 133 | + end |
| 134 | +end |
| 135 | + |
| 136 | +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 |
| 137 | + S = size(Q) |
| 138 | + @assert reduce((&), inds .<= S) |
| 139 | + for (dim, index) in enumerate(inds) |
| 140 | + if index == 1 |
| 141 | + _inds = inds[setdiff(1:N, dim)] |
| 142 | + if (1 ∈ _inds) | reduce((|), S[setdiff(1:N, dim)] .== _inds) |
| 143 | + return zero(T) |
| 144 | + else |
| 145 | + return Q.lower[dim][(_inds.-1)...] |
| 146 | + end |
| 147 | + elseif index == S[dim] |
| 148 | + _inds = inds[setdiff(1:N, dim)] |
| 149 | + if (1 ∈ _inds) | reduce((|), S[setdiff(1:N, dim)] .== _inds) |
| 150 | + return zero(T) |
| 151 | + else |
| 152 | + return Q.upper[dim][(_inds.-1)...] |
| 153 | + end |
| 154 | + end |
| 155 | + end |
| 156 | + return Q.u[(inds.-1)...] |
| 157 | +end |
0 commit comments