Merge pull request #160 from xtalax/multidimbc

Final BC touchups

Author: ChrisRackauckas

src/derivative_operators/BC_operators.jl

@@ -8,6 +8,10 @@ Robin, General, and in general Neumann, Dirichlet and Bridge BCs are all affine
 """
 abstract type AffineBC{T} <: AtomicBC{T} end
 
+struct NeumannBC{N} end
+struct Neumann0BC{N} end
+struct DirichletBC{N} end
+struct Dirichlet0BC{N} end
 """
 q = PeriodicBC{T}()
 
@@ -18,6 +22,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,6 +54,7 @@ struct ComposedMultiDimBC{T, B<:AtomicBC{T}, N,M} <: MultiDimensionalBC{T, N}
     BCs::Vector{Array{B, M}}
 end
 
+
 """
 A multiple dimensional BC, supporting arbitrary BCs at each boundary point.
 To construct an arbitrary BC, pass an Array of BCs with dimension `N-1`, if `N` is the dimensionality of your domain `u`
@@ -82,25 +83,37 @@ 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)
+struct MultiDimBC{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(BC::B, s, dim::Integer) where  {B<:AtomicBC} = MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)]))
+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)]))
 
 #Extra constructor to make a set of BC operators that extend an atomic BC Operator to the whole domain
 #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{dim}(T,s) where dim = MultiDimBC{dim}(PeriodicBC(T), s)
+PeriodicBC(T,s) = MultiDimBC(PeriodicBC(T), s)
 
-PeriodicBC{T}(s) where T = MultiDimBC(PeriodicBC{T}(), s)
-
+NeumannBC{dim}(α::NTuple{2,T}, dx, order, s) where {T,dim} = 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,dim} = RobinBC{dim}((one(T), zero(T), αl), (one(T), zero(T), αr), 1.0, 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) where {dim} = DirichletBC{dim}(zero(T), zero(T), s)
 Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)
+
+Neumann0BC{dim}(T::Type, dx, order, s) where {dim} = 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,dim} = 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,dim} = 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)])
 
 
@@ -118,7 +131,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)
@@ -131,6 +144,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})
 

test/2D_3D_fast_multiplication.jl

@@ -313,7 +313,7 @@ end
     # Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
     A = Ly2+Ly4
     mul!(M_temp, A, M)
-    @test M_temp ≈ ((Ly2*M) + (Ly4*M))
+    @test isapprox(M_temp, ((Ly2*M) + (Ly4*M)), atol = 1e-7)
 
     # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
     A += Ly3

test/MultiDimBC_test.jl

@@ -10,14 +10,14 @@ 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}()
+q2 = PeriodicBC(Float64)
 
 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)
-Qy = MultiDimBC(BCy, 2)
+Qx = MultiDimBC{1}(BCx)
+Qy = MultiDimBC{2}(BCy)
 
 Ax = Qx*A
 Ay = Qy*A
@@ -45,16 +45,22 @@ 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}()
+q2 = PeriodicBC(Float64)
 
 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(Dirichlet0BC(Float64), size(A), 3) #Test the other constructor
+Qx = MultiDimBC{1}(BCx)
+Qy = MultiDimBC{2}(BCy)
+Qz = Dirichlet0BC{3}(Float64, size(A)) #Test the other constructor
 
+Q1 = (Qx+Qy+Qz)
+Q2 = compose(Qx,Qy,Qz) #test addition combinations
+@test_broken Q1 == Q2 #This fails
+@test all(Q1.BCs .== Q2.BCs) # but this passes so it does actually work
+@test_broken Qtmp = Qx + Qz
+@test_skip Qz+Qx+Qy == Qy+Qtmp
 Ax = Qx*A
 Ay = Qy*A
 Az = Qz*A