Add operator matrix (#56)

* small additions to docstrings

* add operator matrix and some refactoring

* remove spurious file

* format

docs/src/interpolation.md

@@ -252,7 +252,7 @@ Kernels can be composed by using [`SumKernel`](@ref) and [`ProductKernel`](@ref)
 which applies a transformation to the input before evaluating the kernel.
 
 However, you can also define your own kernel. A radial-symmetric kernel is a subtype of [`KernelInterpolation.RadialSymmetricKernel`](@ref), which in
-turn is a subtype of [`KernelInterpolation.AbstractKernel`](@ref) and needs to implement the functions [`phi`](@ref) and [`order`](@ref). Let's define a exponential
+turn is a subtype of [`KernelInterpolation.AbstractKernel`](@ref) and needs to implement the functions [`phi`](@ref) and [`order`](@ref). Let's define an exponential
 kernel with ``\phi(r) = \mathrm{e}^{-r^{1.5}}`` and use it for the interpolation problem above.
 
 ```@example interpolation

src/KernelInterpolation.jl

@@ -47,7 +47,8 @@ export NodeSet, empty_nodeset, separation_distance, dim, eachdim, values_along_d
 export interpolation_kernel, nodeset, coefficients, kernel_coefficients,
        polynomial_coefficients, polynomial_basis, polyvars, system_matrix,
        interpolate, solve_stationary, kernel_inner_product, kernel_norm,
-       TemporalInterpolation
+       kernel_matrix, operator_matrix
+export Interpolation, TemporalInterpolation
 export AliveCallback, SaveSolutionCallback, SummaryCallback
 export vtk_save, vtk_read
 export examples_dir, get_examples, default_example

src/discretization.jl

@@ -45,19 +45,16 @@ function solve_stationary(spatial_discretization::SpatialDiscretization{Dim, Rea
                                                                                             RealT
                                                                                             }
     @unpack equations, nodeset_inner, boundary_condition, nodeset_boundary, kernel = spatial_discretization
-    nodeset = merge(nodeset_inner, nodeset_boundary)
 
-    pd_matrix = pde_matrix(equations, nodeset_inner, nodeset, kernel)
-    b_matrix = kernel_matrix(nodeset_boundary, nodeset, kernel)
-    system_matrix = [pd_matrix
-                     b_matrix]
+    system_matrix = pde_boundary_matrix(equations, nodeset_inner, nodeset_boundary, kernel)
     b = [rhs(nodeset_inner, equations); boundary_condition.(nodeset_boundary)]
     c = system_matrix \ b
 
     # Do not support additional polynomial basis for now
     xx = polyvars(Dim)
     ps = monomials(xx, 0:-1)
-    return Interpolation(kernel, nodeset, c, system_matrix, ps, xx)
+    return Interpolation(kernel, merge(nodeset_inner, nodeset_boundary), c, system_matrix,
+                         ps, xx)
 end
 
 """
@@ -84,14 +81,12 @@ function Semidiscretization(spatial_discretization, initial_condition)
     # whole kernel matrix is not needed for rhs, but for initial condition
     k_matrix = [k_matrix_inner
                 k_matrix_boundary]
-    pd_matrix = pde_matrix(equations, nodeset_inner, nodeset, kernel)
-    b_matrix = kernel_matrix(nodeset_boundary, nodeset, kernel)
-    pde_boundary_matrix = [pd_matrix
-                           b_matrix]
+    pdeb_matrix = pde_boundary_matrix(equations, nodeset_inner, nodeset_boundary, kernel)
+
     m_matrix = [k_matrix_inner
                 zeros(eltype(k_matrix_inner), size(k_matrix_boundary)...)]
     cache = (; kernel_matrix = k_matrix, mass_matrix = m_matrix,
-             pde_boundary_matrix = pde_boundary_matrix)
+             pde_boundary_matrix = pdeb_matrix)
     return Semidiscretization{typeof(initial_condition), typeof(cache)}(spatial_discretization,
                                                                         initial_condition,
                                                                         cache)

src/interpolation.jl

@@ -1,19 +1,5 @@
 abstract type AbstractInterpolation{Kernel, Dim, RealT} end
 
-"""
-    interpolation_kernel(itp)
-
-Return the kernel from an interpolation object.
-"""
-interpolation_kernel(itp::AbstractInterpolation) = itp.kernel
-
-"""
-	nodeset(itp)
-
-Return the node set from an interpolation object.
-"""
-nodeset(itp::AbstractInterpolation) = itp.nodeset
-
 @doc raw"""
     Interpolation
 
@@ -51,6 +37,20 @@ Return the dimension of the input variables of the interpolation.
 """
 dim(::Interpolation{Kernel, Dim, RealT, A}) where {Kernel, Dim, RealT, A} = Dim
 
+"""
+    interpolation_kernel(itp)
+
+Return the kernel from an interpolation object.
+"""
+interpolation_kernel(itp::AbstractInterpolation) = itp.kernel
+
+"""
+	nodeset(itp)
+
+Return the node set from an interpolation object.
+"""
+nodeset(itp::AbstractInterpolation) = itp.nodeset
+
 """
     coefficients(itp::Interpolation)
 

src/io.jl

@@ -7,7 +7,7 @@ or vectors to save the values of the functions at the nodes. The functions can a
 [`KernelInterpolation.Interpolation`](@ref) or directly as vectors. The optional keyword argument `keys` is used to
 specify the names of the data arrays in the VTK file.
 
-See [`add_to_pvd`](@ref)
+See also [`add_to_pvd`](@ref), [`vtk_read`](@ref).
 """
 function vtk_save(filename, nodeset::NodeSet, functions_or_vectors...;
                   keys = "value_" .* string.(eachindex(functions_or_vectors)))
@@ -57,6 +57,8 @@ end
 Read a set of nodes from a VTK file and return them as a [`NodeSet`](@ref). Note that the data
 will always be returned as a 3D [`NodeSet`](@ref), even if the data is 1D or 2D. The point data
 will be returned as a dictionary with the keys being the names of the data arrays in the VTK file.
+
+See also [`vtk_save`](@ref).
 """
 function vtk_read(filename)
     vtk = VTKFile(filename)

src/kernel_matrices.jl

@@ -3,7 +3,7 @@
 
 Return the kernel matrix for the nodeset and kernel. The kernel matrix is defined as
 ```math
-A_{ij} = K(x_i, \xi_j),
+    A_{ij} = K(x_i, \xi_j),
 ```
 where ``x_i`` are the nodes in the `nodeset1`, ``\xi_j`` are the nodes in the `nodeset2`,
 and ``K`` the `kernel`.
@@ -25,7 +25,7 @@ end
 
 Return the kernel matrix for the nodeset and kernel. The kernel matrix is defined as
 ```math
-A_{ij} = K(x_i, x_j),
+    A_{ij} = K(x_i, x_j),
 ```
 where ``x_i`` are the nodes in the `nodeset` and ``K`` the `kernel`.
 """
@@ -38,7 +38,7 @@ end
 
 Return the polynomial matrix for the nodeset and polynomials. The polynomial matrix is defined as
 ```math
-A_{ij} = p_j(x_i),
+    A_{ij} = p_j(x_i),
 ```
 where ``x_i`` are the nodes in the `nodeset` and ``p_j`` the polynomials.
 """
@@ -56,23 +56,69 @@ function polynomial_matrix(nodeset, ps)
 end
 
 @doc raw"""
-    pde_matrix(equations, nodeset1, nodeset2, kernel)
+    pde_matrix(diff_op_or_pde, nodeset1, nodeset2, kernel)
 
-Compute the matrix of a partial differential equation with a given kernel. The matrix is defined as
+Compute the matrix of a partial differential equation (or differential operator) with a given kernel. The matrix is defined as
 ```math
-    A_{ij} = \mathcal{L}K(x_i, \xi_j),
+    (\tilde A_\mathcal{L})_{ij} = \mathcal{L}K(x_i, \xi_j),
 ```
-where ``\mathcal{L}`` is the differential operator defined by the `equations`, ``K`` the `kernel`, ``x_i`` are the nodes
+where ``\mathcal{L}`` is the differential operator (defined by the `equations`), ``K`` the `kernel`, ``x_i`` are the nodes
 in `nodeset1` and ``\xi_j`` are the nodes in `nodeset2`.
 """
-function pde_matrix(equations, nodeset1, nodeset2, kernel)
+function pde_matrix(diff_op_or_pde, nodeset1, nodeset2, kernel)
     n = length(nodeset1)
     m = length(nodeset2)
     A = Matrix{eltype(nodeset1)}(undef, n, m)
     for i in 1:n
         for j in 1:m
-            A[i, j] = equations(kernel, nodeset1[i], nodeset2[j])
+            A[i, j] = diff_op_or_pde(kernel, nodeset1[i], nodeset2[j])
         end
     end
     return A
 end
+
+@doc raw"""
+    pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, kernel)
+
+Compute the matrix of a partial differential equation (or differential operator) with a given kernel. The matrix is defined as
+```math
+    A_\mathcal{L} = \begin{pmatrix}\tilde A_\mathcal{L}\\\tilde A}\end{pmatrix},
+```
+where ``\tilde A_\mathcal{L}`` is the matrix of the differential operator (defined by the `equations`) for the inner nodes `x_i`:
+```math
+    (\tilde A_\mathcal{L})_{ij} = \mathcal{L}K(x_i, \xi_j),
+```
+and ``\tilde A`` is the kernel matrix for the boundary nodes:
+```math
+    \tilde A_{ij} = K(x_i, \xi_j),
+```
+where ``\mathcal{L}`` is the differential operator (defined by the `equations`), ``K`` the `kernel`, ``x_i`` are the nodes
+in `nodeset_boundary` and ``\xi_j`` are the nodes in union of `nodeset_inner` and `nodeset_boundary`.
+
+See also [`pde_matrix`](@ref) and [`kernel_matrix`](@ref).
+"""
+function pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, kernel)
+    nodeset = merge(nodeset_inner, nodeset_boundary)
+    pd_matrix = pde_matrix(diff_op_or_pde, nodeset_inner, nodeset, kernel)
+    b_matrix = kernel_matrix(nodeset_boundary, nodeset, kernel)
+    return [pd_matrix
+            b_matrix]
+end
+
+@doc raw"""
+    operator_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, kernel)
+
+Compute the operator matrix `L` discretizing `\mathcal{L}` for a given kernel. The operator matrix is defined as
+```math
+    L = A_\mathcal{L} * A^{-1},
+```
+where ``A_\mathcal{L}`` is the matrix of the differential operator (defined by the `equations`), and ``A`` the kernel matrix.
+
+See also [`pde_boundary_matrix`](@ref) and [`kernel_matrix`](@ref).
+"""
+function operator_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, kernel)
+    nodeset = merge(nodeset_inner, nodeset_boundary)
+    A = kernel_matrix(nodeset, kernel)
+    A_L = pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, kernel)
+    return A_L / A
+end

src/kernels/kernels.jl

@@ -20,5 +20,9 @@ Returns the canonical, human-readable name for the given system of equations.
 get_name(kernel::AbstractKernel) = string(nameof(typeof(kernel))) * "{" *
                                    string(dim(kernel)) * "}"
 
+function (kernel::AbstractKernel)(x)
+    return kernel(x, zero(x))
+end
+
 include("radialsymmetric_kernel.jl")
 include("special_kernel.jl")

src/kernels/radialsymmetric_kernel.jl

@@ -12,6 +12,10 @@ The kernel is then defined by
 ```math
     K(x, y) = \Phi(x - y).
 ```
+
+A `RadialSymmetricKernel` can be evaluated at two points `x` and `y` by calling
+`kernel(x, y)` or at a single point `x` by calling `kernel(x)`, which implicitly
+sets `y` to zero.
 """
 abstract type RadialSymmetricKernel{Dim} <: AbstractKernel{Dim} end
 

src/nodes.jl

@@ -281,7 +281,7 @@ end
 If `n` is integer, create a [`NodeSet`](@ref) with `n` homogeneously distributed nodes in every dimension each of dimension
 `dim` inside a hypercube defined by the bounds `x_min` and `x_max`. If `n` is a `Tuple` of length `dim`,
 then use as many nodes in each dimension as described by `n`. The resulting `NodeSet` will have
-``n^{\textrm{dim}}`` respectively ``\prod_{j = 1}{\textrm{dim}}n_j`` points. If the bounds are given as single values,
+``n^{\textrm{dim}}`` respectively ``\prod_{j = 1}^{\textrm{dim}}n_j`` points. If the bounds are given as single values,
 they are applied for each dimension. If they are `Tuple`s of size `dim`, the hypercube has the according bounds.
 If `dim` is not given explicitly, it is inferred by the lengths of `n`, `x_min` and `x_max` if possible.
 """

src/visualization.jl

@@ -3,12 +3,12 @@
     if Dim == 1
         x = LinRange(x_min, x_max, N)
         title --> get_name(kernel)
-        x, kernel.(Ref(0.0), x)
+        x, kernel.(x)
     elseif Dim == 2
         nodeset = homogeneous_hypercube(N, x_min, x_max; dim = 2)
         x = unique(values_along_dim(nodeset, 1))
         y = unique(values_along_dim(nodeset, 2))
-        z = reshape(kernel.(Ref([0.0, 0.0]), nodeset), (N, N))
+        z = reshape(kernel.(nodeset), (N, N))
         seriestype --> :heatmap # :contourf
         title --> get_name(kernel)
         x, y, z
@@ -20,7 +20,7 @@ end
 @recipe function f(x::AbstractVector, kernel::AbstractKernel)
     xguide --> "r"
     title --> get_name(kernel)
-    x, kernel.(Ref(0.0), x)
+    x, kernel.(x)
 end
 
 @recipe function f(nodeset::NodeSet, kernel::AbstractKernel)
@@ -34,7 +34,7 @@ end
         seriestype --> :scatter
         label --> "nodes"
         title --> get_name(kernel)
-        x, y, kernel.(Ref([0.0, 0.0]), nodeset)
+        x, y, kernel.(nodeset)
     else
         @error("Plotting a kernel is only supported for dimension up to 2, but the set has dimension $(dim(nodeset))")
     end

test/test_unit.jl

@@ -841,6 +841,22 @@ include("test_util.jl")
         for (node, value) in zip(nodeset_boundary, values_boundary)
             @test isapprox(itp(node), value, atol = 1e-14)
         end
+        nodes = nodeset(itp)
+        # Test if u = A * c
+        A = kernel_matrix(nodes, kernel)
+        c = coefficients(itp)
+        u1_values = A * c
+        u2_values = itp.(nodes)
+        for (u1_val, u2_val) in zip(u1_values, u2_values)
+            @test isapprox(u1_val, u2_val, atol = 1e-14)
+        end
+        # Test if L * u = b
+        L = operator_matrix(pde, nodeset_inner, nodeset_boundary, kernel)
+        b = [f1.(nodeset_inner, Ref(pde)); g1.(nodeset_boundary)]
+        b_test = L * u2_values
+        for (b_val, b_test_val) in zip(b, b_test)
+            @test isapprox(b_val, b_test_val, atol = 1e-12)
+        end
         # Test if the solution is close to the analytical solution in other points
         x = [0.1, 0.08]
         @test isapprox(itp(x), u1(x), atol = 0.12)
@@ -875,7 +891,7 @@ include("test_util.jl")
             @test isapprox(titp(t, node), value, atol = 1e-14)
             @test isapprox(titp(t)(node), titp(t, node))
         end
-        t = sol.t[end]
+        # Test if u = A * c
         A = Matrix(semi.cache.kernel_matrix)
         c = sol(t)
         u1_values = A * c