Allow imposing different boundary conditions on each side (#70)

* Allow different BCs in RecombineMatrix * Update .gitignore * Generalise RecombineMatrix constructors * Move natural BC stuff to separate file * Inverse order of arguments in constructors * Fix recombination tests * Add tests * Cleanup * v0.16 * Update ci.yml * Fix new test
jipolanco · May 26, 2023 · e19d73c · e19d73c
1 parent 430783a
commit e19d73c
Show file tree

Hide file tree

Showing 20 changed files with 352 additions and 297 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -21,7 +21,7 @@ jobs:
         experimental: [false]
         version:
           - '1.8'
-          - '~1.9.0-0'
+          - '1.9'
         os:
           - ubuntu-latest
         arch:
@@ -60,7 +60,7 @@ jobs:
       - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '~1.9.0-0'
+          version: '1.9'
       - name: Install dependencies
         run: |
           julia --project=docs -e '

diff --git a/.gitignore b/.gitignore
@@ -1 +1,2 @@
 Manifest.toml
+LocalPreferences.toml
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "BSplineKit"
 uuid = "093aae92-e908-43d7-9660-e50ee39d5a0a"
 authors = ["Juan Ignacio Polanco <[email protected]>"]
-version = "0.15.2"
+version = "0.16.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,7 +18,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
-ArrayLayouts = "0.8, 1"
+ArrayLayouts = "1"
 BandedMatrices = "0.17"
 FastGaussQuadrature = "0.5"
 LazyArrays = "0.22, 1"

diff --git a/README.md b/README.md
@@ -66,7 +66,7 @@ by this package.
 
   ```julia
   bc = Derivative(0) + 3Derivative(1)
-  R = RecombinedBSplineBasis(bc, B)  # satisfies u ∓ 3u' = 0 on the left/right boundary
+  R = RecombinedBSplineBasis(B, bc)  # satisfies u ∓ 3u' = 0 on the left/right boundary
   ```
 
 - Construct [mass matrix](https://en.wikipedia.org/wiki/Mass_matrix) and

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -66,7 +66,7 @@ See the examples in the sidebar for more details.
 
   ```julia
   bc = Derivative(0) + 3Derivative(1)
-  R = RecombinedBSplineBasis(bc, B)  # satisfies u ∓ 3u' = 0 on the left/right boundary
+  R = RecombinedBSplineBasis(B, bc)  # satisfies u ∓ 3u' = 0 on the left/right boundary
   ```
 
 - Construct [mass matrix](https://en.wikipedia.org/wiki/Mass_matrix) and

diff --git a/examples/heat.jl b/examples/heat.jl
@@ -99,7 +99,7 @@ fig
 # recombined bases for many different BCs.
 # In this example we generate a basis satisfying homogeneous Neumann BCs:
 
-R = RecombinedBSplineBasis(Derivative(1), B)
+R = RecombinedBSplineBasis(B, Derivative(1))
 
 fig = Figure()
 ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = L"ϕ_i(x)")
@@ -438,7 +438,7 @@ fig
 hi_res = let
     knots_in = range(-1, 1; length = 101)
     B = BSplineBasis(BSplineOrder(6), knots_in)
-    R = RecombinedBSplineBasis(Derivative(1), B)
+    R = RecombinedBSplineBasis(B, Derivative(1))
     θ₀_spline = approximate(θ₀, R)
     u_init = coefficients(θ₀_spline)
     xcol = collocation_points(R)

diff --git a/notebooks/airy.jl b/notebooks/airy.jl
@@ -38,7 +38,7 @@ breaks =
 B = BSplineBasis(BSplineOrder(6), copy(breaks))
 
 # ╔═╡ e4b55060-6936-11eb-336d-7f6d75edabff
-R = RecombinedBSplineBasis(Derivative(0), B)
+R = RecombinedBSplineBasis(B, Derivative(0))
 
 # ╔═╡ 06d6d560-6937-11eb-03e3-5b2be0c46741
 xs = collocation_points(B);

diff --git a/notebooks/fitting.jl b/notebooks/fitting.jl
@@ -36,7 +36,7 @@ B0 = BSplineBasis(BSplineOrder(4), xbreaks)
 
 # ╔═╡ 823cb81a-79c6-11eb-1e4d-75eb02469abf
 # Impose u' = 0 at the boundaries
-B = RecombinedBSplineBasis(Derivative(1), B0)
+B = RecombinedBSplineBasis(B0, Derivative(1))
 # B = B0
 
 # ╔═╡ 2b263ef0-79af-11eb-3243-eb0ac08d1a46

diff --git a/src/Recombinations/Recombinations.jl b/src/Recombinations/Recombinations.jl
@@ -27,6 +27,7 @@ import ..BSplines:
     boundaries, knots, order, evaluate, nonzero_in_segment, support, summary_basis
 
 include("matrices.jl")
+include("natural.jl")
 include("bases.jl")
 include("splines.jl")
 

diff --git a/src/Recombinations/bases.jl b/src/Recombinations/bases.jl
@@ -97,10 +97,13 @@ For more details on the supported combinations of BCs, see the different
 
 ---
 
-    RecombinedBSplineBasis(op::AbstractDifferentialOp, B::BSplineBasis)
+    RecombinedBSplineBasis(B::BSplineBasis, op::AbstractDifferentialOp)
+    RecombinedBSplineBasis(B::BSplineBasis, op_left, op_right)
 
 Construct `RecombinedBSplineBasis` from B-spline basis `B`, satisfying
-homogeneous boundary conditions associated to the given differential operator.
+homogeneous boundary conditions (BCs) associated to the given differential operator.
+
+The second case allows setting different BCs on each boundary.
 
 For instance, `op = Derivative(0)` and `op = Derivative(1)` correspond to
 homogeneous Dirichlet and Neumann BCs, respectively.
@@ -119,13 +122,13 @@ julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1)
 13-element BSplineBasis of order 4, domain [-1.0, 1.0]
  knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 
-julia> R_neumann = RecombinedBSplineBasis(Derivative(1), B)
+julia> R_neumann = RecombinedBSplineBasis(B, Derivative(1))
 11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
  knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
  BCs left:  (D{1},)
  BCs right: (D{1},)
 
-julia> R_robin = RecombinedBSplineBasis(Derivative(0) + 3Derivative(1), B)
+julia> R_robin = RecombinedBSplineBasis(B, Derivative(0) + 3Derivative(1))
 11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
  knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
  BCs left:  (D{0} + 3 * D{1},)
@@ -134,7 +137,8 @@ julia> R_robin = RecombinedBSplineBasis(Derivative(0) + 3Derivative(1), B)
 
 ---
 
-    RecombinedBSplineBasis(ops, B::BSplineBasis)
+    RecombinedBSplineBasis(B::BSplineBasis, ops)
+    RecombinedBSplineBasis(B::BSplineBasis, ops_left, ops_right)
 
 Construct `RecombinedBSplineBasis` simultaneously satisfying homogeneous BCs
 associated to multiple differential operators.
@@ -184,7 +188,7 @@ For instance, `ops = (Derivative(1), Derivative(0))` fails with an error.
 julia> ops = (Derivative(0), Derivative(1));
 
 
-julia> R1 = RecombinedBSplineBasis(ops, B)
+julia> R1 = RecombinedBSplineBasis(B, ops)
 9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
  knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
  BCs left:  (D{0}, D{1})
@@ -193,13 +197,13 @@ julia> R1 = RecombinedBSplineBasis(ops, B)
 julia> ops = (Derivative(0), Derivative(1) - 4Derivative(2));
 
 
-julia> R2 = RecombinedBSplineBasis(ops, B)
+julia> R2 = RecombinedBSplineBasis(B, ops)
 9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
  knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
  BCs left:  (D{0}, D{1} + -4 * D{2})
  BCs right: (D{0}, D{1} + -4 * D{2})
 
-julia> R3 = RecombinedBSplineBasis(Natural(), B)
+julia> R3 = RecombinedBSplineBasis(B, Natural())
 11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
  knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
  BCs left:  (D{2},)
@@ -213,14 +217,20 @@ struct RecombinedBSplineBasis{
     B :: Parent   # original B-spline basis
     M :: RMatrix  # basis recombination matrix
 
-    function RecombinedBSplineBasis(ops, B::BSplineBasis{k,T}) where {k,T}
+    function RecombinedBSplineBasis(B::BSplineBasis{k,T}, ops_left, ops_right) where {k,T}
         Parent = typeof(B)
-        M = RecombineMatrix(ops, B)
+        M = RecombineMatrix(B, ops_left, ops_right)
         RMatrix = typeof(M)
         new{k,T,Parent,RMatrix}(B, M)
     end
 end
 
+# Same BCs on both boundaries.
+RecombinedBSplineBasis(B::BSplineBasis, ops) = RecombinedBSplineBasis(B, ops, ops)
+
+# This is the old constructor and will be deprecated/removed in the future.
+RecombinedBSplineBasis(ops, B::BSplineBasis) = RecombinedBSplineBasis(B, ops)
+
 BSplines.has_parent_basis(::Type{<:RecombinedBSplineBasis}) = true
 
 Base.:(==)(A::RecombinedBSplineBasis, B::RecombinedBSplineBasis) =
@@ -245,18 +255,6 @@ function Base.show(io::IO, R::RecombinedBSplineBasis)
     end
 end
 
-"""
-    RecombinedBSplineBasis(ops, args...; kwargs...)
-
-Construct [`RecombinedBSplineBasis`](@ref) from B-spline basis, satisfying
-homogeneous boundary conditions associated one or more differential operators.
-
-This variant does not require a previously constructed [`BSplineBasis`](@ref).
-Arguments are passed to the `BSplineBasis` constructor.
-"""
-RecombinedBSplineBasis(ops, args...; kwargs...) =
-    RecombinedBSplineBasis(ops, BSplineBasis(args...; kwargs...))
-
 """
     parent(R::RecombinedBSplineBasis)