From d36232ec77636aced70ce67b05fb20d9ca05b150 Mon Sep 17 00:00:00 2001
From: Kristof Cools <krcools@gmail.com>
Date: Thu, 14 May 2026 13:14:56 +0200
Subject: [PATCH] Proof of concept on how AdaptiveCrossApproximation can be
 used by BEAST.jl

---
 .gitignore        |  1 +
 example/efie2.jl  | 44 ++++++++++++++++++++++++++++++++
 example/mfie.jl   | 64 +++++++++++++++++++++++++++++++++++++++++++++++
 example/pmchwt.jl | 62 +++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 171 insertions(+)
 create mode 100644 example/efie2.jl
 create mode 100644 example/mfie.jl
 create mode 100644 example/pmchwt.jl

diff --git a/.gitignore b/.gitignore
index c3e81f4..591a4cd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@
 /docs/build/
 lcov.info
 *.vscode
+/dev/
diff --git a/example/efie2.jl b/example/efie2.jl
new file mode 100644
index 0000000..eb4d7a0
--- /dev/null
+++ b/example/efie2.jl
@@ -0,0 +1,44 @@
+using LinearAlgebra
+using CompScienceMeshes
+using BEAST
+using ParallelKMeans
+using H2Trees
+using AdaptiveCrossApproximation
+using Krylov
+using PlotlyJS
+
+Γ = meshsphere(1.0, 0.04);
+X = raviartthomas(Γ);
+@show numfunctions(X)
+
+κ, η = 2pi, 1.0;
+t = Maxwell3D.singlelayer(; wavenumber=κ);
+E = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ);
+Eₜ = (n × E) × n;
+
+ttree = KMeansTree(X.pos, 2; minvalues=100)
+tree = BlockTree(ttree, ttree)
+
+function hassemble(op, X, Y; kwargs...)
+    return HMatrix(op, X, Y, tree;
+        tol=1e-4,
+        maxrank=40,
+        isnear=AdaptiveCrossApproximation.isnear(),
+        spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),
+    )
+end
+
+@hilbertspace k
+@hilbertspace j
+
+a = t[k,j]
+A = assemble(a, ∏(X), ∏(X); materialize=hassemble);
+b = assemble(Eₜ[k], ∏(X));
+
+A⁻¹ = BEAST.GMRESSolver(A; reltol=1e-4, maxiter=1000)
+u = A⁻¹ * b
+u = BEAST.FEMFunction(u, ∏(X))
+
+Plot(mesh3d(u[j], colorscale=:Viridis))
+
+
diff --git a/example/mfie.jl b/example/mfie.jl
new file mode 100644
index 0000000..584bf42
--- /dev/null
+++ b/example/mfie.jl
@@ -0,0 +1,64 @@
+using LinearAlgebra
+using CompScienceMeshes
+using BEAST
+using ParallelKMeans
+using H2Trees
+using AdaptiveCrossApproximation
+using Krylov
+using PlotlyJS
+
+Γ = meshsphere(1.0, 0.04);
+X = raviartthomas(Γ);
+Y = buffachristiansen(Γ);
+@show numfunctions(X)
+
+κ, η = 2pi, 1.0;
+ϵ, μ = 1/η, η;
+t = Maxwell3D.singlelayer(; wavenumber=κ);
+E = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ);
+H = -1/(im*μ) * curl(E)
+
+Eₜ = (n × E) × n;
+Hₜ = (n × H) × n;
+
+@hilbertspace k
+@hilbertspace j
+
+function materialize(op, X, Y; kwargs...)
+    if op isa BEAST.IntegralOperator
+        Xtree = KMeansTree(X.pos, 2; minvalues=100)
+        Ytree = KMeansTree(Y.pos, 2; minvalues=100)
+        tree = BlockTree(Xtree, Ytree)
+        return HMatrix(op, X, Y, tree;
+            tol=1e-4,
+            maxrank=40,
+            isnear=AdaptiveCrossApproximation.isnear(),
+            spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),
+        )
+    end
+    return BEAST.assemble(op, X, Y; kwargs...)
+end
+
+K = Maxwell3D.doublelayer(; wavenumber=κ);
+N = BEAST.NCross();
+
+a = K[k,j] + 0.5*N[k,j]
+l = Hₜ[k]
+
+A = assemble(a, ∏(Y), ∏(X); materialize=materialize);
+b = assemble(l, ∏(Y));
+
+A⁻¹ = BEAST.GMRESSolver(A; reltol=1e-4, maxiter=1000)
+u = A⁻¹ * b
+u = BEAST.FEMFunction(u, ∏(X))
+
+Plot(mesh3d(u[j], colorscale=:Viridis))
+
+Xtree = KMeansTree(X.pos, 2; minvalues=100)
+Ytree = KMeansTree(Y.pos, 2; minvalues=100)
+tree = BlockTree(Xtree, Ytree)
+@which HMatrix(K, Y, X, tree;
+    tol=1e-4,
+    maxrank=40,
+    isnear=AdaptiveCrossApproximation.isnear(),
+    spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),)
\ No newline at end of file
diff --git a/example/pmchwt.jl b/example/pmchwt.jl
new file mode 100644
index 0000000..a69266e
--- /dev/null
+++ b/example/pmchwt.jl
@@ -0,0 +1,62 @@
+using LinearAlgebra
+using CompScienceMeshes
+using BEAST
+using ParallelKMeans
+using H2Trees
+using AdaptiveCrossApproximation
+using Krylov
+using PlotlyJS
+
+Γ = meshsphere(1.0, 0.06);
+RT = raviartthomas(Γ);
+X = RT × RT
+
+κ, η = 2pi, 1.0;
+T = Maxwell3D.singlelayer(; wavenumber=κ);
+K = Maxwell3D.doublelayer(; wavenumber=κ);
+N = BEAST.NCross();
+
+κ′, η′ = 2.0 * κ,  η
+T′ = Maxwell3D.singlelayer(; wavenumber=κ′);
+K′ = Maxwell3D.doublelayer(; wavenumber=κ′);
+
+E = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ);
+H = -1/(im*μ) * curl(E)
+
+e = (n × E) × n;
+h = (n × H) × n;
+
+
+function materialize(op, X, Y; kwargs...)
+    ttree = KMeansTree(X.pos, 2; minvalues=100)
+    tree = BlockTree(ttree, ttree)
+    H = HMatrix(op, X, Y, tree;
+        tol=1e-4,
+        maxrank=40,
+        isnear=AdaptiveCrossApproximation.isnear(),
+        spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),
+    )
+    AdaptiveCrossApproximation.storage(H)
+    return H
+end
+
+@hilbertspace p q
+@hilbertspace m j
+
+α, α′ = 1/η, 1/η′
+a = (
+    α*T[p,m]+α′*T′[p,m] + K[p,j]+K′[p,j]
+    -K[q,m]-K′[q,m] + η*T[q,j]+η′*T′[q,j]
+)
+b = -h[p] + e[q] 
+
+𝐀 = assemble(a, X, X; materialize=materialize);
+𝐛 = assemble(b, X);
+
+𝐀⁻¹ = BEAST.GMRESSolver(𝐀; reltol=1e-4, maxiter=1000)
+𝐮 = 𝐀⁻¹ * 𝐛
+𝐮 = BEAST.FEMFunction(𝐮, X)
+
+Plot(mesh3d(𝐮[j], colorscale=:Viridis))
+
+