Add Lavery splines (#23) (#83)

* add lavery splines

* improve performance of 1D Lavery spline. Some restructuring.

* first prototype of 2D implementation

* add contributers

* rename new 1d test file

* simplify the evaluation of the 2D Lavery spline

* small adjusts to 1D files

* apply Horner rule to 2D Lavery spline evaluation

* prototype of the 2D Lavery spline

* adjust new 2D Lavery example

* add 2D test for Lavery spline

* fix formatting of new data file

* add extra test points to exercise all pieces of Lavery evaluation routine

* fix docs refs

* fix broken links in docstrings

* new implementation of the 2D Lavery splines that is much faster

* fix formatting in the new examples

* fix formatting in 1D files

* fix formatting in 2D files

* merge constraints to remove some allocations

* move new Lavery implementation into a package extension

* add NEWS entry with new functionality

* added to docstring for the 1D spline

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>

* adjust names and constructors

* update README

* fix formatting

* try to fix broken docs build

* another attempt to fix docs build

* reorganize and docs now build

* works on generic real type

* adjust reference.md for docs build

* run formatter

* remove old comment

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>

* use lambda and same default everywhere

---------

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>

Author: 3 people

AUTHORS.md

@@ -18,7 +18,10 @@ The following people contributed major additions or modifications to TrixiBottom
 are listed in alphabetical order:
 
 * Maximilian D. Bertrand
+* Logan Eriksson
+* Oscar Jemsson
 * Joshua Lampert
+* Vincent Marks
 * Hendrik Ranocha
 * Michael Schlottke-Lakemper
 * Andrew Winters

NEWS.md

@@ -9,6 +9,7 @@ for human readability.
 
 #### Added
 
+- Added shape preserving splines and interpolation using [JuMP.jl](https://github.com/jump-dev/JuMP.jl) and [HiGHS.jl](https://github.com/jump-dev/HiGHS.jl) [#83]
 - Added scattered radial basis function interpolation using [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl) [#73]
 - Implementation of visualization routines `plot_topography` and `plot_topography_with_interpolation_knots`
 

Project.toml

@@ -10,15 +10,20 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 KernelInterpolation = "95686e15-3c94-4443-8f08-76d06d509f7b"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [extensions]
+JuMPHiGHSExt = ["JuMP", "HiGHS"]
 KernelInterpolationExt = "KernelInterpolation"
 MakieExt = "Makie"
 
 [compat]
 Downloads = "1.6"
+HiGHS = "1.0"
+JuMP = "1.0"
 KernelInterpolation = "0.3.5"
 LinearAlgebra = "1"
 Makie = "0.21, 0.22, 0.24"

README.md

@@ -38,6 +38,7 @@ This package contains the following main functionalities:
 - Setting up a B-spline interpolation structure in one and two dimensions which contains all the relevant information to define a B-spline interpolation function with additional specifications
 - Using the B-spline structure to set up a B-spline interpolation function
 - Setting up a radial basis function interpolation based on scattered data leveraging [KernelInterpolation.jl](https://github.com/JoshuaLampert/KernelInterpolation.jl)
+- Setting up a shape preserving cubic spline structure (Lavery spline) for interpolation that avoids generating new extrema.
 
 A detailed description of the functionalities can be found in the [documentation](https://trixi-framework.github.io/TrixiBottomTopography.jl/stable/) to this package.
 

docs/Project.toml

@@ -2,6 +2,8 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Changelog = "5217a498-cd5d-4ec6-b8c2-9b85a09b6e3e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 KernelInterpolation = "95686e15-3c94-4443-8f08-76d06d509f7b"
 OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
@@ -15,6 +17,8 @@ path = ".."
 CairoMakie = "0.13, 0.15"
 Changelog = "1.1"
 Documenter = "1.3"
+HiGHS = "1.0"
+JuMP = "1.0"
 KernelInterpolation = "0.3.5"
 OrdinaryDiffEqLowStorageRK = "1.3"
 Trixi = "0.12, 0.13, 0.14"

docs/make.jl

@@ -5,6 +5,8 @@ using Changelog: Changelog
 # To load the extensions
 import CairoMakie
 import KernelInterpolation
+import JuMP
+import HiGHS
 
 # Copy list of authors to not need to synchronize it manually.
 # Since the authors header exists twice we create a unique identifier for the docs section.
@@ -94,7 +96,8 @@ rm(joinpath(@__DIR__, "src", "changelog_tmp.md"))
 makedocs(;
          modules = [TrixiBottomTopography,
              Base.get_extension(TrixiBottomTopography, :MakieExt),
-             Base.get_extension(TrixiBottomTopography, :KernelInterpolationExt)],
+             Base.get_extension(TrixiBottomTopography, :KernelInterpolationExt),
+             Base.get_extension(TrixiBottomTopography, :JuMPHiGHSExt)],
          authors = "Andrew R. Winters <andrew.ross.winters@liu.se>, Michael Schlottke-Lakemper <michael@sloede.com>",
          sitename = "TrixiBottomTopography.jl",
          format = Documenter.HTML(;

docs/src/reference.md

@@ -23,3 +23,13 @@ TrixiBottomTopography.RBFInterpolation
 TrixiBottomTopography.RBFInterpolation1D
 TrixiBottomTopography.RBFInterpolation2D
 ```
+
+## JuMPHiGHS.jl extension
+```@docs
+TrixiBottomTopography.LaverySpline1D(::Vector{T}, ::Vector{T}) where {T<:Real}
+TrixiBottomTopography.LaverySpline1D(::String)
+TrixiBottomTopography.LaverySpline2D(::Vector{T}, ::Vector{T}, ::Matrix{T}) where {T<:Real}
+TrixiBottomTopography.LaverySpline2D(::String)
+TrixiBottomTopography.spline_interpolation(::LaverySpline1D, ::Number)
+TrixiBottomTopography.spline_interpolation(::LaverySpline2D, ::Number, ::Number)
+```

examples/lavery_spline_1d.jl

@@ -0,0 +1,40 @@
+##############################################################################
+# Script which uses the functionalities implemented in TrixiBottomTopography #
+# to plot a one dimensional Lavery spline of data with sharp transitions.    #
+##############################################################################
+
+# Include packages
+using TrixiBottomTopography
+using JuMP
+using HiGHS
+
+# Load in data from one of the examples from Eriksson and Jemsson thesis
+# https://liu.diva-portal.org/smash/get/diva2:1918338/FULLTEXT01.pdf
+xData = Vector(0.0:29.0)
+yData = [0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 2.0, 4.0, 3.0, 2.0,
+    1.0, 0.0, 9.0, 9.0, 9.0, 0.0, 0.0, 0.0, 6.0, 7.0,
+    5.0, 6.0, 7.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+# Lavery spline that does not generate new extrema
+spline_struct = LaverySpline1D(xData, yData)
+# Define Lavery spline interpolation function
+spline_func(x) = spline_interpolation(spline_struct, x)
+
+# Evaluate the cubic Lavery spline on a new set of nodes and plot
+if isdefined(Main, :Makie)
+    # Define interpolation points
+    n = 500
+    x_int_pts = Vector(LinRange(spline_struct.x[1], spline_struct.x[end], n))
+
+    # Get interpolated values
+    y_int_pts = spline_func.(x_int_pts)
+
+    # Get the original interpolation knots
+    x_knots = spline_struct.x
+    y_knots = spline_func.(x_knots)
+
+    plot_topography_with_interpolation_knots(x_int_pts, y_int_pts, x_knots, y_knots;
+                                             xlabel = "x",
+                                             ylabel = "y",
+                                             legend_position = :lt)
+end

examples/lavery_spline_2d.jl

@@ -0,0 +1,46 @@
+##############################################################################
+# Script which uses the functionalities implemented in TrixiBottomTopography #
+# to plot a two dimensional Lavery spline of data with sharp transitions.    #
+##############################################################################
+
+# Include packages
+using TrixiBottomTopography
+using JuMP
+using HiGHS
+
+# Load in data from one of the examples from Eriksson and Jemsson thesis
+# https://liu.diva-portal.org/smash/get/diva2:1918338/FULLTEXT01.pdf
+x_knots = Vector(0.0:9.0)
+y_knots = Vector(0.0:9.0)
+z_knots = [1.0 1.0 1.0 0.0 0.0 4.0 4.0 4.0 4.0 1.0;
+           1.0 1.0 1.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0;
+           1.0 1.0 1.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0;
+           2.0 2.0 2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
+           3.0 3.0 3.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
+           4.0 4.0 4.0 1.0 0.0 0.0 2.0 2.0 0.0 0.0;
+           1.0 0.0 3.0 1.0 0.0 0.0 2.0 2.0 0.0 0.0;
+           1.0 0.0 3.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
+           1.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
+           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
+
+# Lavery spline that does not generate new extrema
+spline_struct = LaverySpline2D(x_knots, y_knots, z_knots)
+# Define Lavery spline interpolation function
+spline_func(x, y) = spline_interpolation(spline_struct, x, y)
+
+# Evaluate the bicubic Lavery spline on a new set of nodes and plot
+if isdefined(Main, :Makie)
+    # Define interpolation points
+    n = 150
+    x_int_pts = Vector(LinRange(spline_struct.x[1], spline_struct.x[end], n))
+    y_int_pts = Vector(LinRange(spline_struct.y[1], spline_struct.y[end], n))
+
+    # Get interpolated matrix
+    z_int_pts = evaluate_two_dimensional_interpolant(spline_func, x_int_pts, y_int_pts)
+
+    plot_topography_with_interpolation_knots(x_int_pts, y_int_pts, z_int_pts,
+                                             x_knots, y_knots, z_knots;
+                                             xlabel = "x",
+                                             ylabel = "y",
+                                             zlabel = "z")
+end