Kernel split implementation and tests

Project.toml

@@ -4,20 +4,24 @@ version = "0.2.0"
 
 [deps]
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
+ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 ElementaryPDESolutions = "88a69b33-da0f-4502-8c8c-d680cf4d883b"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HAdaptiveIntegration = "eaa5ad34-b243-48e9-b09c-54bc0655cecf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Richardson = "708f8203-808e-40c0-ba2d-98a6953ed40d"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+SpecialMatrices = "928aab9d-ef52-54ac-8ca1-acd7ca42c160"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 
@@ -42,11 +46,13 @@ IntiQPGreenExt = "QPGreen"
 
 [compat]
 Bessels = "0.2"
+ClassicalOrthogonalPolynomials = "0.15.7"
 DataStructures = "0.18 - 0.19"
 ElementaryPDESolutions = "0.2 - 0.3"
 FMM2D = "0.2"
 FMM3D = "1"
 FMMLIB2D = "0.3"
+FastGaussQuadrature = "1.1.0"
 ForwardDiff = "0.10, 1"
 Gmsh = "0.3"
 HAdaptiveIntegration = "0.2"
@@ -56,6 +62,7 @@ LinearMaps = "3"
 Makie = "0.21 - 0.24"
 Meshes = "0.42 - 0.55"
 NearestNeighbors = "0.4"
+OffsetArrays = "1.17.0"
 Pkg = "1"
 Printf = "1"
 QPGreen = "1.0.1"
@@ -64,6 +71,7 @@ Richardson = "1"
 Scratch = "1"
 SparseArrays = "1"
 SpecialFunctions = "2"
+SpecialMatrices = "3.1.0"
 StaticArrays = "1"
 TOML = "1"
 julia = "1.9"

docs/src/refs.bib

@@ -143,3 +143,36 @@ @article{bernardi1989optimal
   year={1989},
   publisher={SIAM}
 }
+
+@article{helsing2008evaluation,
+  title={On the evaluation of layer potentials close to their sources},
+  author={Helsing, Johan and Ojala, Rikard},
+  journal={Journal of Computational Physics},
+  volume={227},
+  number={5},
+  pages={2899--2921},
+  year={2008},
+  publisher={Elsevier}
+}
+
+@article{helsing2015variants,
+  title={Variants of an explicit kernel-split panel-based Nystr{\"o}m discretization scheme for Helmholtz boundary value problems},
+  author={Helsing, Johan and Holst, Anders},
+  journal={Advances in Computational Mathematics},
+  volume={41},
+  number={3},
+  pages={691--708},
+  year={2015},
+  publisher={Springer}
+}
+
+@article{fryklund2022adaptive,
+  title={An adaptive kernel-split quadrature method for parameter-dependent layer potentials},
+  author={Fryklund, Fredrik and Klinteberg, Ludvig af and Tornberg, Anna-Karin},
+  journal={Advances in Computational Mathematics},
+  volume={48},
+  number={2},
+  pages={12},
+  year={2022},
+  publisher={Springer}
+}

src/Inti.jl

@@ -50,6 +50,9 @@ include("nystrom.jl")
 include("adaptive_correction.jl")
 include("bdim.jl")
 include("vdim.jl")
+# include("kernel_split.jl")
+# include("adaptive_kernel_split.jl")
+include("ksplit.jl")
 
 # some zero-argument methods for the Inti's gmsh extension
 include("gmsh_api.jl")

src/api.jl

@@ -11,7 +11,7 @@ const COMPRESSION_METHODS = [:none, :hmatrix, :fmm]
 Available correction methods for the singular and nearly-singular integrals in
 [`Inti`](@ref).
 """
-const CORRECTION_METHODS = [:none, :dim, :adaptive]
+const CORRECTION_METHODS = [:none, :dim, :adaptive, :ksplit, :ksplit_adaptive]
 
 """
     single_double_layer(; op, target, source::Quadrature, compression,
@@ -163,6 +163,108 @@ function single_double_layer(;
         correction_kw = Base.structdiff(correction, NamedTuple{(:method,)})
         δS = adaptive_correction(Sop; correction_kw...)
         δD = adaptive_correction(Dop; correction_kw...)
+    elseif correction.method == :ksplit
+        # Extract the required geometric args from the correction tuple
+        geo_args = (
+            source_quad_connectivity = correction.connectivity,
+            source_el = correction.elements,
+            velocity_fn = correction.velocity_fn,
+            curvature_fn = correction.curvature_fn,
+            boundary_inv = correction.boundary_inv,
+        )
+
+        # Extract the ksplit-specific optional keywords
+        ksplit_kw_names =
+            (:n_panel_corr, :maxdist, :target_location, :parametric_length)
+        ksplit_kwargs = filter(kv -> kv[1] in ksplit_kw_names, pairs(correction))
+
+        if target === source
+            if !haskey(ksplit_kwargs, :n_panel_corr)
+                error(
+                    "For :ksplit correction with `target === source`, you must provide `n_panel_corr` (e.g., `correction = (method=:ksplit, n_panel_corr=3, ...)`).",
+                )
+            end
+        else # target !== source
+            if !haskey(ksplit_kwargs, :maxdist)
+                error(
+                    "For :ksplit correction with `target !== source`, you must provide `maxdist` (e.g., `correction = (method=:ksplit, maxdist=0.1, ...)`).",
+                )
+            end
+        end
+
+        δS = kernel_split_correction(
+            op,
+            source,
+            geo_args...,
+            Sop,
+            target;
+            layer_type = :single,
+            ksplit_kwargs...,
+        )
+        δD = kernel_split_correction(
+            op,
+            source,
+            geo_args...,
+            Dop,
+            target;
+            layer_type = :double,
+            ksplit_kwargs...,
+        )
+    elseif correction.method == :ksplit_adaptive
+        # Extract the required geometric args
+        geo_args = (
+            source_quad_connectivity = correction.connectivity,
+            source_el = correction.elements,
+            velocity_fn = correction.velocity_fn,
+            curvature_fn = correction.curvature_fn,
+            boundary_inv = correction.boundary_inv,
+        )
+
+        # Extract the adaptive ksplit-specific keywords
+        ksplit_adaptive_kw_names = (
+            :n_panel_corr,
+            :maxdist,
+            :target_location,
+            :Cε,
+            :Rε,
+            :parametric_length,
+            :affine_preimage,
+        )
+        ksplit_adaptive_kwargs = filter(kv -> kv[1] in ksplit_adaptive_kw_names, pairs(correction))
+
+        if target === source
+            if !haskey(ksplit_adaptive_kwargs, :n_panel_corr)
+                error(
+                    "For :ksplit correction with `target === source`, you must provide `n_panel_corr` (e.g., `correction = (method=:ksplit_adaptive, n_panel_corr=3, ...)`).",
+                )
+            end
+        else # target !== source
+            if !haskey(ksplit_adaptive_kwargs, :maxdist)
+                error(
+                    "For :ksplit correction with `target !== source`, you must provide `maxdist` (e.g., `correction = (method=:ksplit_adaptive, maxdist=0.1, ...)`).",
+                )
+            end
+        end
+
+        # Call your adaptive_kernel_split_correction function
+        δS = adaptive_kernel_split_correction(
+            op,
+            source,
+            geo_args...,
+            Sop,
+            target;
+            layer_type = :single,
+            ksplit_adaptive_kwargs...
+        )
+        δD = adaptive_kernel_split_correction(
+            op,
+            source,
+            geo_args...,
+            Dop,
+            target;
+            layer_type = :double,
+            ksplit_adaptive_kwargs...
+        )
     else
         error("Unknown correction method. Available options: $CORRECTION_METHODS")
     end