Merge pull request #76

Minor improvements when preparing WAVES talk

.gitignore

@@ -10,6 +10,7 @@
 .vscode
 /docs/src/examples/generated/
 /sandbox/*
+/talks/*
 
 .markdownlint.json
 Inti.code-workspace

docs/make.jl

@@ -84,6 +84,7 @@ makedocs(;
         canonical = "https://IntegralEquations.github.io/Inti.jl",
         size_threshold = 2 * 2^20, # 2 MiB
         size_threshold_warn = 1 * 2^20, # 1 MiB
+        sidebar_sitename = false,
     ),
     pages = [
         "Home" => "index.md",

docs/src/examples/graded_mesh.jl

@@ -0,0 +1,103 @@
+using Inti
+using StaticArrays
+using Meshes
+using GLMakie
+using IterativeSolvers
+using FMMLIB2D
+using LinearAlgebra
+
+## Parameters
+k = 10π
+λ = 2π / k
+meshsize = λ / 10
+TEST = false
+PLOT = true # can be quite slow at large k
+
+## Create the geometry as an equilateral
+Γ = let
+    P = 2
+    χ = Inti.kress_change_of_variables_periodic(P)
+    p1 = SVector(0.0, 0.0)
+    p2 = SVector(1.0, 0.0)
+    p3 = SVector(0.5, sqrt(3) / 2)
+
+    l1 = Inti.parametric_curve(
+        t -> p1 + χ(t) * (p2 - p1) + 0.0 * SVector(0, sin(2π * t)),
+        0.0,
+        1.0,
+    )
+    l2 = Inti.parametric_curve(t -> p2 + χ(t) * (p3 - p2), 0.0, 1.0)
+    l3 = Inti.parametric_curve(t -> p3 + χ(t) * (p1 - p3), 0.0, 1.0)
+    Inti.Domain(l1, l2, l3)
+end
+
+msh = Inti.meshgen(Γ; meshsize)
+Q = Inti.Quadrature(msh; qorder = 3)
+
+##
+
+## Create integral operators
+pde = Inti.Helmholtz(; dim = 2, k)
+S, D = Inti.single_double_layer(;
+    pde,
+    target = Q,
+    source = Q,
+    compression = (method = :fmm, tol = 1e-6),
+    correction = (method = :dim,),
+)
+L = I / 2 + D - im * k * S
+
+𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q)
+
+## Test with a source inside for validation
+if TEST
+    G = Inti.SingleLayerKernel(pde)
+    uₑ = x -> G(x, SVector(0.4, 0.4))
+    rhs = map(q -> uₑ(q.coords), Q)
+    σ = gmres(L, rhs; restart = 100, maxiter = 100, abstol = 1e-8, verbose = true)
+    uₕ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x)
+    xₜ = SVector(2, 0)
+    @show uₕ(xₜ), uₑ(xₜ)
+end
+
+## Test with a planewave
+θ = 0 * π / 4
+uᵢ = x -> exp(im * k * dot(x, SVector(cos(θ), sin(θ))))
+rhs = map(q -> -uᵢ(q.coords), Q)
+
+σ = gmres(L, rhs; restart = 1000, maxiter = 1000, abstol = 1e-8, verbose = true)
+uₕ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x)
+@info "Done solving..."
+
+## Visualize a solution
+if PLOT
+    xx = -1:λ/10:2
+    yy = -1:λ/10:2
+    target = [SVector(x, y) for x in xx, y in yy]
+
+    Spot, Dpot = Inti.single_double_layer(;
+        pde,
+        target = target,
+        source = Q,
+        compression = (method = :fmm, tol = 1e-4),
+        correction = (method = :none, target_location = :outside, maxdist = meshsize),
+    )
+
+    vals = Dpot * σ - im * k * Spot * σ
+    vals = map(
+        i -> Inti.isinside(target[i], Q) ? NaN + im * NaN : vals[i] + uᵢ(target[i]),
+        1:length(target),
+    )
+    vals = reshape(vals, length(xx), length(yy))
+
+    heatmap(
+        xx,
+        yy,
+        real.(vals);
+        colorrange = (-1, 1),
+        colormap = :inferno,
+        axis = (aspect = DataAspect(),),
+    )
+    viz!(msh[Γ]; showsegments = true)
+    Makie.current_figure()
+end

src/api.jl

@@ -172,13 +172,8 @@ function single_double_layer(;
         S = axpy!(true, δS, Smat)
         D = axpy!(true, δD, Dmat)
     elseif compression.method == :hmatrix
-        if target === source
-            S = axpy!(true, δS, Smat)
-            D = axpy!(true, δD, Dmat)
-        else
-            S = LinearMap(Smat) + LinearMap(δS)
-            D = LinearMap(Dmat) + LinearMap(δD)
-        end
+        S = LinearMap(Smat) + LinearMap(δS)
+        D = LinearMap(Dmat) + LinearMap(δD)
     elseif compression.method == :fmm
         S = Smat + LinearMap(δS)
         D = Dmat + LinearMap(δD)

src/mesh.jl

@@ -601,3 +601,10 @@ end
     centers = map(center, els)
     return inrange(balltree, centers, tol)
 end
+
+"""
+    Domain(f::Function, msh::AbstractMesh)
+
+Call `Domain(f, ents)` on `ents = entities(msh).`
+"""
+Domain(f::Function, msh::AbstractMesh) = Domain(f, entities(msh))

src/polynomials.jl

@@ -103,7 +103,8 @@ function lagrange_basis(nodes, sp::PolynomialSpace)
     # compute the matrix of coefficients of the lagrange polynomials over the
     # monomomial basis
     V = hcat([basis(x) for x in nodes]...)
-    C = SArray(MArray(V) \ I)
+    # convert to an array to use the backslash operator
+    C = typeof(V)(Matrix(V) \ I)
     lag_basis = x -> C * basis(x)
     return lag_basis
 end

src/utils.jl

@@ -408,3 +408,30 @@ function _rotation_matrix_3d(rot)
     Rz = @SMatrix [cos(rot[3]) sin(rot[3]) 0; -sin(rot[3]) cos(rot[3]) 0; 0 0 1]
     return Rz * Ry * Rx
 end
+
+"""
+    kress_change_of_variables(P)
+
+Return a change of variables mapping `[0,1]` to `[0,1]` with the property that
+the first `P-1` derivatives of the transformation vanish at `x=0`.
+"""
+function kress_change_of_variables(P)
+    v = x -> (1 / P - 1 / 2) * ((1 - x))^3 + 1 / P * ((x - 1)) + 1 / 2
+    return x -> 2v(x)^P / (v(x)^P + v(2 - x)^P)
+end
+
+"""
+    kress_change_of_variables_periodic(P)
+
+Like [`kress_change_of_variables`](@ref), this change of variables maps the interval `[0,1]` onto
+itself, but the first `P` derivatives of the transformation vanish at **both**
+endpoints (thus making it a periodic function).
+
+This change of variables can be used to *periodize* integrals over the interval
+`[0,1]` by mapping the integrand into a new integrand that vanishes (to order P)
+at both endpoints.
+"""
+function kress_change_of_variables_periodic(P)
+    v = (x) -> (1 / P - 1 / 2) * ((1 - 2x))^3 + 1 / P * ((2x - 1)) + 1 / 2
+    return x -> v(x)^P / (v(x)^P + v(1 - x)^P)
+end

src/vdim.jl

@@ -88,13 +88,13 @@ function vdim_correction(
                 iszero(center) || error("SHIFT is not implemented for non-zero center")
                 L̃ = [f((q.coords - c) / r) for f in p, q in view(source, jglob)]
                 S = change_of_basis(multiindices, p, c, r)
-                F = lu(L̃)
+                F = svd(L̃)
                 @debug (vander_cond = max(vander_cond, cond(L̃))) maxlog = 0
                 @debug (shift_norm = max(shift_norm, norm(S))) maxlog = 0
                 @debug (vander_norm = max(vander_norm, norm(L̃))) maxlog = 0
             else
                 L = [f(q.coords) for f in p, q in view(source, jglob)]
-                F = lu(L)
+                F = svd(L)
                 @debug (vander_cond = max(vander_cond, cond(L))) maxlog = 0
                 @debug (shift_norm = max(shift_norm, 1)) maxlog = 0
                 @debug (vander_norm = max(vander_norm, norm(L))) maxlog = 0