Added a {dim} version of all BC constructors,

Fleshed out PartiallyComposedBCs

Author: xtalax

src/derivative_operators/BC_operators.jl

@@ -18,6 +18,7 @@ Creates a periodic boundary condition, where the lower index end of some u is ex
 It is not reccomended to concretize this BC type in to a BandedMatrix, since the vast majority of bands will be all 0s. SpatseMatrix concretization is reccomended.
 """
 struct PeriodicBC{T} <: AtomicBC{T}
+    PeriodicBC(T::Type) = new{T}()
 end
 
 """

src/derivative_operators/multi_dim_bc_operators.jl

@@ -54,42 +54,60 @@ struct ComposedMultiDimBC{T, B<:AtomicBC{T}, N,M} <: MultiDimensionalBC{T, N}
     BCs::Vector{Array{B, M}}
 end
 
-struct PartiallyComposedMultiDimBC{T,B<:MultiDimDirectionalBC,N,L} <:MultiDimensionalBC{T,N}
+struct PartiallyComposedMultiDimBC{T,B<:Union{MultiDimDirectionalBC, Nothing},N, is_extended} <:MultiDimensionalBC{T,N}
     BCs::Vector{B}
 end
 
-function Base.:+(Q1::MultiDimDirectionalBC{T,B1,D1,N,M}, Q2::MultiDimDirectionalBC{T,B2,D2,N,M}) where {T,B1, B2,D1,D2,N,M}
-    @assert D1 != D2 "BCs should all extend along a unique axis"
-    return PartiallyComposedMultiDimBC{T,Union{typeof(Q1),typeof(Q2)}, N,2}(Array[Q1,Q2])	
-end
+numextended
 
-function Base.:+(Q1::MultiDimDirectionalBC{T,B1,D1,2,1}, Q2:: MultiDimDirectionalBC{T,B2,D2,2,1}) where {T,B1,B2,D1,D2}
-    return compose(Q1,Q2)
-end
 
-function Base.:+(Q::PartiallyComposedMultiDimBC{T,B,N,L}, Q1::MultiDimDirectionalBC) where {T,B,N,M,L}
-    new_BCs = vcat(Q.BCs,Q1)
-    if L+1 == N
-        return compose(new_BCs...)
-    elseif !(getaxis(Q1) in getaxis.(Q.BCs))
-	return PartiallyComposedMultiDimBC{T,eltype(new_BCs), N,L-1}(new_BCs)
-    else
-	throw("BCs should all extend along a unique axis")
+function Base._:+(BCs::MultiDimDirectionalBC{T}...) where T
+	D = getaxis.(BCs)
+	BCs = BCs[sortperm([D...])]
+	B = typeof.(BCs)
+	dimensionalities = ndims.(BCs)
+	for n in dimensionalities[2:end]
+		@assert dimensionalites[1] == n "There are "
+
+	for d in D
+        @assert length(setdiff(D, d)) == (N-1) "There are multiple boundary conditions that extend along $d - make sure every dimension has a unique extension"
     end
+	is_extended = falses(length(ndims(BCs[1]))
+	is_extended[D] .= true
+	A = Vector{Union{B..., Nothing}}(undef, ndims(BCs[1]))
+	A[!is_extended] .= Nothing()
+	for (i,BC) in enumerate(BCs)
+		A[is_extended][i] = BC
+	end
+	return PartiallyComposedMultiDimBC{T, eltype(A), N, is_extended}(A)
+end
+
+
+function Base.:+(Q::PartiallyComposedMultiDimBC{T,B,N,L}, Q1::MultiDimDirectionalBC{T, B2, D, N, M}) where {T,D, B1, B2,N,M,L}
+    new_BCs = Q.BCs
+	new_BCs[D] = Q1
+	is_extended = L
+	is_extended[D] = true
+    @assert !(getaxis(Q1) in getaxis.(Q.BCs)) "BCs should all extend along a unique axis"
+	return PartiallyComposedMultiDimBC{T,eltype(new_BCs), N,is_extended}(new_BCs)
+
 end
 
 Base.:+(Q1::MultiDimDirectionalBC, Q::PartiallyComposedMultiDimBC) = +(Q,Q1)
 
-function Base.:+(Q1::PartiallyComposedMultiDimBC{T,B1,N,L1} Q2::PartiallyComposedMultiDimBC{T,B2,N,L2) where {T,B1,B2,N,L1,L2}
-    new_BCs = vcat(Q1.BCs, Q2.BCs)
-    if L1 + L2 < N
-	@assert length(setdiff(1:N, getaxis.(new_BCs))) == (N-L1-L2) "BCs should all extend along a unique axis"
-	return PartiallyComposedMultiDimBC{T,Union{B1,B2}, N, L1+ L2}(new_BCs)
-    elseif L1+L2 == N
-        return compose(new_BCs...)
-    else
-        throw("Too many supplied BCs for the number of dimensions")
-    end
+function Base.:+(Q1::PartiallyComposedMultiDimBC{T,B1,N,L1} Q2::PartiallyComposedMultiDimBC{T,B2,N,L2}) where {T,B1,B2,N,L1,L2}
+    new_BCs = deepcopy(Q1.BCs)
+	Lnew = L1 .| L2
+	for (i,q) in enumerate(Q2.BCs)
+		if q <: MultiDimensionalBC
+			if new_BCs[i] <: Nothing
+				new_BCs[i] = q
+			else
+				throw("BCs should all extend along a unique axis")
+			end
+		end
+	end
+	return PartiallyComposedMultiDimBC{T,Union{B1,B2}, N, Lnew}(new_BCs)
 end
 
 """
@@ -121,6 +139,11 @@ For Neumann0BC, please use
 where T is the element type of the domain to be extended
 """
 struct MultiDimBC{N} end
+struct NeumannBC{N} end
+struct Neumann0BC{N} end
+struct DirichletBC{N} end
+struct Dirichlet0BC{N} end
+
 MultiDimBC{dim}(BC::Array{B,N}) where {N, B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC[1]), B, dim, N+1, N}(BC)
 #s should be size of the domain
 MultiDimBC{dim}(BC::B, s) where  {B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)]))
@@ -129,17 +152,26 @@ MultiDimBC{dim}(BC::B, s) where  {B<:AtomicBC, dim} = MultiDimDirectionalBC{gett
 #Only valid in the uniform grid case!
 MultiDimBC(BC::B, s) where {B<:AtomicBC} = Tuple([MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])
 
-# Additional constructors for cases when the BC is the same for all boundarties 
-
-PeriodicBC{T}(s) where T = MultiDimBC(PeriodicBC{T}(), s)
+# Additional constructors for cases when the BC is the same for all boundarties
+PeriodicBC{dim}(T,s) where T = MultiDimBC{dim}(PeriodicBC(T), s)
+PeriodicBC(T,s) where T = MultiDimBC(PeriodicBC(T), s)
 
+NeumannBC{dim}(α::NTuple{2,T}, dx, order, s) where T = RobinBC{dim}((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order, s)
 NeumannBC(α::NTuple{2,T}, dxyz, order, s) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dxyz, order, s)
+
+DirichletBC{dim}(αl::T, αr::T, s) where T = RobinBC{dim}((one(T), zero(T), αl), (one(T), zero(T), αr), [ones(T, si) for si in s], 2.0, s)
 DirichletBC(αl::T, αr::T, s) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), [ones(T, si) for si in s], 2.0, s)
 
+Dirichlet0BC{dim}(T::Type, s) = DirichletBC{dim}(zero(T), zero(T), s)
 Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)
+
+Neumann0BC{dim}(T::Type, dx, order, s) = NeumannBC{dim}((zero(T), zero(T)), dx, order, s)
 Neumann0BC(T::Type, dxyz, order, s) = NeumannBC((zero(T), zero(T)), dxyz, order, s)
 
+RobinBC{dim}(l::NTuple{3,T}, r::NTuple{3,T}, dx, order, s) where {T} = MultiDimBC{dim}(RobinBC(l, r, dx, order), s)
 RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), perpindex(s,dim))) for dim in 1:length(s)])
+
+GeneralBC{dim}(αl::AbstractVector{T}, αr::AbstractVector{T}, dx, order, s) where {T} = MultiDimBC{dim}(GeneralBC(αl, αr, dx, order), s)
 GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, GeneralBC{T}, dim, length(s), length(s)-1}(fill(GeneralBC(αl, αr, dxyz[dim], order),perpindex(s,dim))) for dim in 1:length(s)])
 
 
@@ -157,7 +189,7 @@ Creates a ComposedMultiDimBC operator, Q, that extends every boundary when appli
 
 Qx Qy and Qz can be passed in any order, as long as there is exactly one BC operator that extends each dimension.
 """
-function compose(BCs...)
+function compose(BCs::MultiDimDirectionalBC...)
     T = gettype(BCs[1])
     N = ndims(BCs[1])
     Ds = getaxis.(BCs)
@@ -170,6 +202,8 @@ function compose(BCs...)
     ComposedMultiDimBC{T, Union{eltype.(BC.BCs for BC in BCs)...}, N,N-1}([condition.BCs for condition in BCs])
 end
 
+Base.:+(BCs::MultiDimDirectionalBC...) = compose(BCs...)
+
 """
 Qx, Qy,... = decompose(Q::ComposedMultiDimBC{T,N,M})
 
@@ -198,5 +232,18 @@ function Base.:*(Q::ComposedMultiDimBC{T, B, N, K}, u::AbstractArray{T, N}) wher
     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
 
-function Base.:*(Q::PartiallyComposedMultiDimBC, u::AbstractArray{T,N} where {T,N}
-    return Tuple([MultiDimDirectionalBC
+
+
+function Base.:*(Q::PartiallyComposedMultiDimBC{T,B,N,is_extended}, u::AbstractArray{T,N} where {T,N}
+	lower = Vector{Union{Array{T,N-1}, Nothing}}(undef, N)
+	upper = Vector{Union{Array{T,N-1}, Nothing}}(undef, N)
+	for (i,q) in Q.BCs
+		if !(q <: Nothing)
+			lower[i], upper[i] = slice_rmul(q, u)
+		else
+			lower[i] = Nothing()
+			upper[i] = Nothing()
+		end
+	end
+	return PartiallyComposedBoundaryPaddedArray{}(lower,upper,u)
+end