Add singular quadrature tests and fix

Manifest.toml

@@ -69,6 +69,12 @@ git-tree-sha1 = "8041575f021cba5a099a456b4163c9a08b566a02"
 uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 version = "1.1.0"
 
+[[FastGaussQuadrature]]
+deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"]
+git-tree-sha1 = "5829b25887e53fb6730a9df2ff89ed24baa6abf6"
+uuid = "442a2c76-b920-505d-bb47-c5924d526838"
+version = "0.4.7"
+
 [[IfElse]]
 git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 InverseHankelFunction = "bec4b5d4-d8c4-11e9-2b25-9be465e100f1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

src/OptimalPMLTransformations.jl

@@ -48,6 +48,17 @@ export HankelSeries
 
 export Rip2D, classify_outer_boundary, find_rips
 
+export gausslegendreunit,
+       int_gauss,
+       int_gauss_2d,
+       gausslegendreunittrans,
+       int_gauss_trans,
+       int_gauss_trans_2d,
+       gausslegendretrans_mid,
+       int_gauss_trans_mid,
+       gausslegendretrans_mid,
+       int_gauss_trans_2d_mid
+
 include("coordinates/coordinates.jl")
 
 include("pml_geometries/pml_geometries.jl")

src/integration/integration.jl

@@ -0,0 +1,2 @@
+
+include("singular_quadrature.jl")

src/integration/singular_quadrature.jl

@@ -1,6 +1,5 @@
 
-using FastGaussQuadrature
-using LinearAlgebra
+import FastGaussQuadrature: gausslegendre
 
 "N Gauss-Legendre nodes and weights, but rescaled for integration from 0 to 1"
 function gausslegendreunit(N)

test/runtests.jl

@@ -4,3 +4,5 @@ include("geometry_tests.jl")
 include("planar_wave_tests.jl")
 
 include("single_hankel_tests.jl")
+
+include("singular_quadrature_tests.jl")

test/singular_quadrature_tests.jl

@@ -0,0 +1,94 @@
+
+import Test: @test, @testset
+using OptimalPMLTransformations
+
+import Base: one
+one(x,y) = one(x)
+
+@testset "Singular integration" begin
+
+    @testset "Unit function integration" begin
+
+        for N in 2:7:23
+
+            begin
+                nodes, weights = gausslegendreunit(N)
+                @test all(0 .<= nodes .<= 1)
+                @test sum(weights) ≈ 1
+            end
+
+            @test int_gauss(one, N) ≈ 1
+
+            @test int_gauss_2d(one, N) ≈ 1
+
+            begin
+                nodes, weights = gausslegendreunittrans(N, 0.5)
+                @test all(map(n -> all(0 .<= n .<= 1), nodes))
+                @test sum(weights) ≈ 1
+            end
+
+            @test int_gauss_trans(one, N; a=0.5) ≈ 1
+
+            @test int_gauss_trans_2d(one, N; a=0.5) ≈ 1
+
+            begin
+                a = 0.5
+                x_crit = 0.3
+                nodes, weights = gausslegendreunittrans(N,a)
+                tnodes, tweights = gausslegendretrans_mid(N,nodes,weights,x_crit)
+                @test all(0 .<= nodes .<= 1)
+                @test sum(weights) ≈ 1
+            end
+
+            @test int_gauss_trans_mid(one, N) ≈ 1
+
+            begin
+                a = 0.5
+                x_crit = (0.3, 0.6)
+                nodes, weights = gausslegendreunittrans(N,a)
+                tnodes, tweights = gausslegendretrans_mid(N,nodes,weights,x_crit)
+                @test all(map(n -> all(0 .<= n .<= 1), nodes))
+                @test sum(weights) ≈ 1
+            end
+
+            @test int_gauss_trans_2d_mid(one, N) ≈ 1
+
+        end
+
+    end
+
+    @testset "Integrable singularity integration" begin
+
+        a = 0.5
+        f_1d(x_crit) = x -> complex(x-x_crit)^-a * (1 + 10(x-x_crit) - 8(x-x_crit)^2)
+        antiderivative_of_f_1d(x_crit) = x -> complex(x-x_crit)^-a * ((x-x_crit)/(1-a) + 10(x-x_crit)^2/(2-a) - 8(x-x_crit)^3/(3-a))
+
+        N = 10
+
+        x_crit = 0.0
+        @test int_gauss_trans(f_1d(x_crit), N; a) ≈ antiderivative_of_f_1d(x_crit)(1)
+
+        x_crit = 0.7
+        @test int_gauss_trans_mid(f_1d(x_crit), N; a, x_crit) ≈ antiderivative_of_f_1d(x_crit)(1) - antiderivative_of_f_1d(x_crit)(0)
+
+        f_2d(xy_crit) = (x,y) -> sqrt((x-xy_crit[1])^2+(y-xy_crit[2])^2)^-a * (1 + 10x - 8x^2) * (2 - 5x + 2x^2)
+        f_rip(xy_crit) = (x,y) -> begin
+            z = -(x-xy_crit[1]) + (y-xy_crit[2])*im
+            -1/sqrt(z)
+        end
+
+        # Test for no-throw until we work out closed form
+        xy_crit=(0.0,0.0)
+        int_gauss_trans_2d(f_2d(xy_crit), N; a)
+        int_gauss_trans_2d(f_rip(xy_crit), N; a)
+        @test true
+
+        # Test for no-throw until we work out closed form
+        xy_crit=(0.1,0.9)
+        int_gauss_trans_2d_mid(f_2d(xy_crit), N; a, x_crit=xy_crit)
+        int_gauss_trans_2d_mid(f_rip(xy_crit), N; a, x_crit=xy_crit)
+        @test true
+
+    end
+
+end