Extend vdim to support vector-valued kernels

.gitignore

@@ -1,19 +1,5 @@
-*.jl.*.cov
-*.jl.cov
-*.jl.mem
-*.msh
-*.pos
-/Manifest.toml
+Manifest.toml
 /docs/Manifest.toml
-/docs/build/
 /test/Manifest.toml
-.vscode
-/docs/src/pluto-examples/*.md
-/sandbox/*
-/talks/*
+/docs/build/
 /test/vtk_output/*
-
-.markdownlint.json
-Inti.code-workspace
-docs/src/examples/heat.gif
-docs/src/.DS_Store

Project.toml

@@ -55,7 +55,7 @@ HMatrices = "0.2"
 LinearAlgebra = "1"
 LinearMaps = "3"
 Makie = "0.21 - 0.24"
-Meshes = "0.42 - 0.55"
+Meshes = "0.42 - 0.56"
 NearestNeighbors = "0.4"
 OrderedCollections = "1"
 Pkg = "1"

src/Inti.jl

@@ -7,32 +7,42 @@ module Inti
 
 const PROJECT_ROOT = pkgdir(Inti)
 
+import Bessels
+import ElementaryPDESolutions
+import HAdaptiveIntegration
+import SpecialFunctions
+
+import ElementaryPDESolutions: Polynomial
+
 using DataStructures
 using ForwardDiff
 using LinearAlgebra
 using LinearMaps
+using Logging
 using NearestNeighbors
 using Pkg
+using Printf
 using QuadGK
+using Richardson
 using Scratch
 using SparseArrays
 using StaticArrays
-using Printf
 using TOML
-using Richardson
 using OrderedCollections
 
 import ElementaryPDESolutions
 import SpecialFunctions
 import Bessels # faster than SpecialFunctions for Bessel functions with real args
 import HAdaptiveIntegration
 
+
 # helper functions
 include("utils.jl")
 include("blockarray.jl")
 
 # basic interpolation and integration
 include("reference_shapes.jl")
+
 include("polynomials.jl")
 include("reference_interpolation.jl")
 include("quad_rules_tables.jl")

src/api.jl

@@ -176,8 +176,9 @@ function single_double_layer(;
             S = axpy!(true, δS, Smat)
             D = axpy!(true, δD, Dmat)
         else
-            S = LinearMap(Smat) + LinearMap(δS)
-            D = LinearMap(Dmat) + LinearMap(δD)
+            T = default_density_eltype(op)
+            S = LinearMap{T}(Smat) + LinearMap{T}(δS)
+            D = LinearMap{T}(Dmat) + LinearMap{T}(δD)
         end
     elseif compression.method == :fmm
         S = Smat + LinearMap(δS)
@@ -365,7 +366,7 @@ function volume_potential(; op, target, source::Quadrature, compression, correct
     if compression.method ∈ (:hmatrix, :none)
         # TODO: in the hmatrix case, we may want to add the correction directly
         # to the HMatrix so that a direct solver can be later used
-        V = LinearMap(Vmat) + LinearMap(δV)
+        V = LinearMap{default_density_eltype(op)}(Vmat) + LinearMap{default_density_eltype(op)}(δV)
         # if target === source
         #     V = axpy!(true, δV, Vmat)
         # else

src/bdim.jl

@@ -66,6 +66,8 @@ function bdim_correction(
     N = ambient_dimension(source)
     @assert eltype(Dop) == T "eltype of S and D must match"
     m, n = length(target), length(source)
+    # check if we are in debug mode to avoid expensive computations
+    do_debug = debug_mode()
     if isnothing(filter_target_params)
         dict_near = etype_to_nearest_points(target, source; maxdist)
         num_trgs = m
@@ -161,17 +163,21 @@ function bdim_correction(
             # TODO: get ride of all this transposing mumble jumble by assembling
             # the matrix in the correct orientation in the first place
             F = qr!(transpose(Mdata))
-            @debug (imat_cond = max(cond(Mdata), imat_cond)) maxlog = 0
-            @debug (imat_norm = max(norm(Mdata), imat_norm)) maxlog = 0
+            if do_debug
+                imat_cond = max(cond(Mdata), imat_cond)
+                imat_norm = max(norm(Mdata), imat_norm)
+            end
             for i in near_list[n]
                 j = glob_loc_near_trgs[i]
                 Θi .= Θ[j:j, :]
-                @debug (rhs_norm = max(rhs_norm, norm(Θidata))) maxlog = 0
+                if do_debug
+                    rhs_norm = max(rhs_norm, norm(Θidata))
+                end
                 ldiv!(Wdata, F, transpose(Θidata))
-                @debug (
+                if do_debug
                     res_norm = max(norm(Matrix(F) * Wdata - transpose(Θidata)), res_norm)
-                ) maxlog = 0
-                @debug (theta_norm = max(theta_norm, norm(Wdata))) maxlog = 0
+                    theta_norm = max(theta_norm, norm(Wdata))
+                end
                 for k in 1:nq
                     push!(Is, i)
                     push!(Js, jglob[k])

src/entities.jl

@@ -20,16 +20,12 @@ Base.hash(ent::EntityKey, h::UInt) = hash((ent.dim, abs(ent.tag)), h)
 Base.:(==)(e1::EntityKey, e2::EntityKey) = e1.dim == e2.dim && abs(e1.tag) == abs(e2.tag)
 
 # defer some functions on EntityKey to the corresponding GeometricEntity
-for f in (
-        :labels,
-        :boundary,
-        :pushforward,
-        :ambient_dimension,
-        :hasparametrization,
-        :parametrization,
-    )
-    @eval $f(k::EntityKey) = $f(global_get_entity(k))
-end
+labels(k::EntityKey) = labels(global_get_entity(k))
+boundary(k::EntityKey) = boundary(global_get_entity(k))
+pushforward(k::EntityKey) = pushforward(global_get_entity(k))
+ambient_dimension(k::EntityKey) = ambient_dimension(global_get_entity(k))
+hasparametrization(k::EntityKey) = hasparametrization(global_get_entity(k))
+parametrization(k::EntityKey) = parametrization(global_get_entity(k))
 
 function (k::EntityKey)(x)
     hasparametrization(k) || error("$k has no parametrization")

src/kernels.jl

@@ -33,6 +33,17 @@ abstract type AbstractDifferentialOperator{N} end
 
 ambient_dimension(::AbstractDifferentialOperator{N}) where {N} = N
 
+function range_dimension(op::AbstractDifferentialOperator)
+    T = default_density_eltype(op)
+    if T <: Number
+        return 1
+    elseif hasmethod(length, Tuple{Type{T}})
+        return length(T)
+    else
+        error("default_density_eltype($(typeof(op))) = $T does not define length(::Type{$T}); cannot determine range dimension")
+    end
+end
+
 # convenient constructor for e.g. SingleLayerKernel(op,Float64) or DoubleLayerKernel(op,ComplexF64)
 function (::Type{K})(
         op::Op,

src/utils.jl

@@ -243,6 +243,15 @@ function notimplemented()
     return error("not (yet) implemented")
 end
 
+"""
+    debug_mode()
+
+Return `true` if the current logger's minimum enabled level is `Debug` or lower.
+"""
+function debug_mode()
+    return Logging.min_enabled_level(Logging.current_logger()) <= Logging.Debug
+end
+
 """
     MultiIndex{N}