From adba9e7073385ed7c5c6ae633d1a4a03109599b1 Mon Sep 17 00:00:00 2001
From: lvchien <v.chienle@gmail.com>
Date: Tue, 16 Dec 2025 14:27:03 +0100
Subject: [PATCH] Inverse operator struct

---
 src/operator.jl                       | 14 +++++++++
 test/test_assemble_InverseOperator.jl | 43 +++++++++++++++++++++++++++
 2 files changed, 57 insertions(+)
 create mode 100644 test/test_assemble_InverseOperator.jl

diff --git a/src/operator.jl b/src/operator.jl
index 580ef48c..81d64099 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -115,6 +115,20 @@ function assemble(operator::AbstractOperator, test_functions, trial_functions;
     return Z()
 end
 
+struct InverseOperator <: AbstractOperator
+    op::AbstractOperator
+    solver
+end
+
+function InverseOperator(op::AbstractOperator; solver=BEAST.lu)
+    InverseOperator(op, solver)
+end
+
+function assemble(A::InverseOperator, testfns, trialfns; kwargs...)
+    M = assemble(A.op, testfns, trialfns; kwargs...)
+    return A.solver(M)
+end
+
 function assemble(A::LinearMap, testfns, trialfns; kwargs...)
     @assert numfunctions(testfns) == size(A,1)
     @assert numfunctions(trialfns) == size(A,2)
diff --git a/test/test_assemble_InverseOperator.jl b/test/test_assemble_InverseOperator.jl
new file mode 100644
index 00000000..a52f4ea9
--- /dev/null
+++ b/test/test_assemble_InverseOperator.jl
@@ -0,0 +1,43 @@
+@info "Executing test_assemble_InverseOperator.jl"
+
+using CompScienceMeshes
+using LinearAlgebra
+using BEAST
+using Test
+
+fn = joinpath(pkgdir(BEAST), "test", "assets", "sphere45.in")
+Γ1 = CompScienceMeshes.readmesh(fn)
+Γ2 = Mesh([point(x+3.0,y,z) for (x,y,z) in vertices(Γ1)], deepcopy(cells(Γ1)))
+Γ = [Γ1, Γ2]
+numdoms = 2
+
+X = raviartthomas.(Γ)
+Y = buffachristiansen.(Γ)
+
+P = BEAST.DirectProductSpace([Xᵢ×Xᵢ for Xᵢ ∈ X])
+Q = BEAST.DirectProductSpace([Yᵢ×Yᵢ for Yᵢ ∈ Y])
+
+@hilbertspace m j
+@hilbertspace k l
+
+p = BEAST.hilbertspace(:p, numdoms)
+q = BEAST.hilbertspace(:q, numdoms)
+
+N = BEAST.NCross()
+
+iNi = BEAST.InverseOperator(N, solver=BEAST.lu)[k, m] + BEAST.InverseOperator(N, solver=BEAST.GMRES)[l, j]
+iNdiag = sum(iNi[q[i], p[i]] for i in 1:numdoms)
+
+iA = BEAST.assemble(iNdiag, Q, P)
+
+@test iA isa BEAST.LinearMaps.LinearMap
+
+Ni = N[k, m] + N[l, j]
+Ndiag = sum(Ni[q[i], p[i]] for i in 1:numdoms)
+
+A = assemble(Ndiag, Q, P)
+IA = BEAST.lu(Matrix(A))
+
+@test IA isa BEAST.LinearMaps.LinearMap
+
+@test Matrix(iA) ≈ Matrix(IA) atol=1e-6
\ No newline at end of file